]> git.donarmstrong.com Git - lilypond.git/blobdiff - lily/beam-quanting.cc
lilypond-manuals.css: edit color scheme and some spacing
[lilypond.git] / lily / beam-quanting.cc
index 106fffc74134713cbd74fc181ec48b7e44dd2857..f3505772f1a70940c4deb9014d400b9158e485a2 100644 (file)
@@ -1,7 +1,7 @@
 /*
   This file is part of LilyPond, the GNU music typesetter.
 
-  Copyright (C) 1997--2011 Han-Wen Nienhuys <hanwen@xs4all.nl>
+  Copyright (C) 1997--2015 Han-Wen Nienhuys <hanwen@xs4all.nl>
   Jan Nieuwenhuizen <janneke@gnu.org>
 
   LilyPond is free software: you can redistribute it and/or modify
@@ -30,15 +30,22 @@ using namespace std;
 #include "direction.hh"
 #include "directional-element-interface.hh"
 #include "grob.hh"
+#include "grob-array.hh"
+#include "item.hh"
 #include "international.hh"
+#include "interval-minefield.hh"
+#include "least-squares.hh"
 #include "libc-extension.hh"
 #include "main.hh"
+#include "note-head.hh"
 #include "output-def.hh"
 #include "pointer-group-interface.hh"
+#include "spanner.hh"
 #include "staff-symbol-referencer.hh"
 #include "stencil.hh"
 #include "stem.hh"
 #include "warn.hh"
+#include "string-convert.hh"
 
 Real
 get_detail (SCM alist, SCM sym, Real def)
@@ -60,7 +67,10 @@ Beam_quant_parameters::fill (Grob *him)
   REGION_SIZE = get_detail (details, ly_symbol2scm ("region-size"), 2);
 
   // forbidden quants
-  SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0);
+  SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0)
+    // For stems that are non-standard, the forbidden beam quanting
+    // doesn't really work, so decrease their importance.
+    * exp(- 8*fabs (1.0 - robust_scm2double(him->get_property ("length-fraction"), 1.0)));
   STEM_LENGTH_DEMERIT_FACTOR = get_detail (details, ly_symbol2scm ("stem-length-demerit-factor"), 5);
   HORIZONTAL_INTER_QUANT_PENALTY = get_detail (details, ly_symbol2scm ("horizontal-inter-quant"), 500);
 
@@ -73,7 +83,14 @@ Beam_quant_parameters::fill (Grob *him)
 
   // Collisions
   COLLISION_PENALTY = get_detail (details, ly_symbol2scm ("collision-penalty"), 500);
-  COLLISION_PADDING = get_detail (details, ly_symbol2scm ("collision-padding"), 0.5);
+
+  /* For grace notes, beams get scaled down to 80%, but glyphs go down
+     to 63% (magstep -4 for accidentals). To make the padding
+     commensurate with glyph size for grace notes, we take the square
+     of the length fraction, yielding a 64% decrease.
+   */
+  COLLISION_PADDING = get_detail (details, ly_symbol2scm ("collision-padding"), 0.5)
+    * sqr (robust_scm2double(him->get_property ("length-fraction"), 1.0));
   STEM_COLLISION_FACTOR = get_detail (details, ly_symbol2scm ("stem-collision-factor"), 0.1);
 }
 
@@ -127,7 +144,7 @@ Beam_configuration *Beam_configuration::new_config (Interval start,
 Real
 Beam_scoring_problem::y_at (Real x, Beam_configuration const *p) const
 {
-  return p->y[LEFT] + (x - x_span[LEFT]) * p->y.delta () / x_span.delta ();
+  return p->y[LEFT] + x * p->y.delta () / x_span_;
 }
 
 /****************************************************************/
@@ -154,18 +171,7 @@ LY_DEFINE (ly_beam_score_count, "ly:beam-score-count", 0, 0, 0,
 void Beam_scoring_problem::add_collision (Real x, Interval y,
                                           Real score_factor)
 {
-  if (edge_dirs[LEFT] == edge_dirs[RIGHT])
-    {
-      Direction d = edge_dirs[LEFT];
-
-      Real quant_range_y = quant_range[LEFT][-d]
-                           + (x - x_span[LEFT]) * (quant_range[RIGHT][-d] - quant_range[LEFT][-d]) / x_span.delta ();
-
-      if (d * (quant_range_y - minmax (d, y[UP], y[DOWN])) > 0)
-        {
-          return;
-        }
-    }
+  // We used to screen for quant range, but no more.
 
   Beam_collision c;
   c.beam_y_.set_empty ();
@@ -173,183 +179,695 @@ void Beam_scoring_problem::add_collision (Real x, Interval y,
   for (vsize j = 0; j < segments_.size (); j++)
     {
       if (segments_[j].horizontal_.contains (x))
-        c.beam_y_.add_point (segments_[j].vertical_count_ * beam_translation);
+        c.beam_y_.add_point (segments_[j].vertical_count_ * beam_translation_);
       if (segments_[j].horizontal_[LEFT] > x)
         break;
     }
-  c.beam_y_.widen (0.5 * beam_thickness);
+  c.beam_y_.widen (0.5 * beam_thickness_);
 
   c.x_ = x;
 
-  y *= 1 / staff_space;
+  y *= 1 / staff_space_;
   c.y_ = y;
   c.base_penalty_ = score_factor;
   collisions_.push_back (c);
 }
 
-void Beam_scoring_problem::init_collisions (vector<Grob *> grobs)
+void Beam_scoring_problem::init_instance_variables (Grob *me, Drul_array<Real> ys, bool align_broken_intos)
 {
-  Grob *common_x = NULL;
-  segments_ = Beam::get_beam_segments (beam, &common_x);
-  vector_sort (segments_, beam_segment_less);
-  if (common[X_AXIS] != common_x)
+  beam_ = dynamic_cast<Spanner *> (me);
+  unquanted_y_ = ys;
+
+  /*
+    If 'ys' are finite, use them as starting points for y-positions of the
+    ends of the beam, instead of the best-fit through the natural ends of
+    the stems.  Otherwise, we want to do initial slope calculations.
+  */
+  do_initial_slope_calculations_ = false;
+  for (LEFT_and_RIGHT (d))
+    do_initial_slope_calculations_ |= isinf (unquanted_y_[d]) || isnan (unquanted_y_[d]);
+
+  /*
+    Calculations are relative to a unit-scaled staff, i.e. the quants are
+    divided by the current staff_space_.
+  */
+  staff_space_ = Staff_symbol_referencer::staff_space (beam_);
+  beam_thickness_ = Beam::get_beam_thickness (beam_) / staff_space_;
+  line_thickness_ = Staff_symbol_referencer::line_thickness (beam_) / staff_space_;
+
+  // This is the least-squares DY, corrected for concave beams.
+  musical_dy_ = robust_scm2double (beam_->get_property ("least-squares-dy"), 0);
+
+  vector<Spanner *> beams;
+  align_broken_intos_ = align_broken_intos;
+  if (align_broken_intos_)
     {
-      programming_error ("Disagree on common x. Skipping collisions in beam scoring.");
-      return;
+      Spanner *orig = dynamic_cast<Spanner *> (beam_->original ());
+      if (!orig)
+        align_broken_intos_ = false;
+      else if (!orig->broken_intos_.size ())
+        align_broken_intos_ = false;
+      else
+        beams.insert (beams.end (), orig->broken_intos_.begin (), orig->broken_intos_.end ());
     }
+  if (!align_broken_intos_)
+    beams.push_back (beam_);
 
-  set<Grob *> stems;
-  for (vsize i = 0; i < grobs.size (); i++)
+  /*
+    x_span_ is a single scalar, cumulatively summing the length of all the
+    segments the parent beam was broken-into.
+  */
+  x_span_ = 0.0;
+  is_knee_ = false;
+  normal_stem_count_ = 0;
+  for (vsize i = 0; i < beams.size (); i++)
     {
-      Box b;
-      for (Axis a = X_AXIS; a < NO_AXES; incr (a))
-        b[a] = grobs[i]->extent (common[a], a);
+      extract_grob_set (beams[i], "stems", stems);
+      extract_grob_set (beams[i], "covered-grobs", fake_collisions);
+      vector<Grob *> collisions;
+
+      for (vsize j = 0; j < fake_collisions.size (); j++)
+        if (fake_collisions[j]->get_system () == beams[i]->get_system ())
+          collisions.push_back (fake_collisions[j]);
 
-      Real width = b[X_AXIS].length ();
-      Real width_factor = sqrt (width / staff_space);
+      Grob *common[2];
+      for (int a = 2; a--;)
+        common[a] = common_refpoint_of_array (stems, beams[i], Axis (a));
 
-      Direction d = LEFT;
-      do
-        add_collision (b[X_AXIS][d], b[Y_AXIS], width_factor);
-      while (flip (&d) != LEFT);
+      for (LEFT_and_RIGHT (d))
+        common[X_AXIS] = beams[i]->get_bound (d)->common_refpoint (common[X_AXIS], X_AXIS);
 
-      Grob *stem = unsmob_grob (grobs[i]->get_object ("stem"));
-      if (stem && Stem::has_interface (stem) && Stem::is_normal_stem (stem))
+      // positions of the endpoints of this beam segment, including any overhangs
+      const Interval x_pos = robust_scm2interval (beams[i]->get_property ("X-positions"),
+                                                  Interval (0.0, 0.0));
+
+      Drul_array<Grob *> edge_stems (Beam::first_normal_stem (beams[i]),
+                                     Beam::last_normal_stem (beams[i]));
+
+      Drul_array<bool> dirs_found (0, 0);
+
+      Real my_y = beams[i]->relative_coordinate (common[Y_AXIS], Y_AXIS);
+
+      Interval beam_width (-1.0, -1.0);
+      for (vsize j = 0; j < stems.size (); j++)
         {
-          stems.insert (stem);
+          Grob *s = stems[j];
+          beam_multiplicity_.push_back (Stem::beam_multiplicity (stems[j]));
+          head_positions_.push_back (Stem::head_positions (stems[j]));
+          is_normal_.push_back (Stem::is_normal_stem (stems[j]));
+
+          Stem_info si (Stem::get_stem_info (s));
+          si.scale (1 / staff_space_);
+          stem_infos_.push_back (si);
+          chord_start_y_.push_back (Stem::chord_start_y (s));
+          dirs_found[si.dir_] = true;
+
+          bool f = to_boolean (s->get_property ("french-beaming"))
+                   && s != edge_stems[LEFT] && s != edge_stems[RIGHT];
+
+          Real y = Beam::calc_stem_y (beams[i], s, common, x_pos[LEFT], x_pos[RIGHT], CENTER,
+                                      Interval (0, 0), f);
+          base_lengths_.push_back (y / staff_space_);
+          stem_xpositions_.push_back (s->relative_coordinate (common[X_AXIS], X_AXIS) - x_pos[LEFT] + x_span_);
+          stem_ypositions_.push_back (s->relative_coordinate (common[Y_AXIS], Y_AXIS) - my_y);
+
+          if (is_normal_.back ())
+            {
+              if (beam_width[LEFT] == -1.0)
+                beam_width[LEFT] = stem_xpositions_.back ();
+              beam_width[RIGHT] = stem_xpositions_.back ();
+            }
         }
+
+      edge_dirs_ = Drul_array<Direction> (CENTER, CENTER);
+      normal_stem_count_ += Beam::normal_stem_count (beams[i]);
+      if (normal_stem_count_)
+        edge_dirs_ = Drul_array<Direction> (stem_infos_[0].dir_,
+                                            stem_infos_.back ().dir_);
+
+      is_xstaff_ = has_interface<Align_interface> (common[Y_AXIS]);
+      is_knee_ |= dirs_found[DOWN] && dirs_found[UP];
+
+      staff_radius_ = Staff_symbol_referencer::staff_radius (beams[i]);
+      edge_beam_counts_ = Drul_array<int>
+                          (Stem::beam_multiplicity (stems[0]).length () + 1,
+                           Stem::beam_multiplicity (stems.back ()).length () + 1);
+
+      // TODO - why are we dividing by staff_space_?
+      beam_translation_ = Beam::get_beam_translation (beams[i]) / staff_space_;
+
+      for (LEFT_and_RIGHT (d))
+        {
+          quant_range_[d].set_full ();
+          if (!edge_stems[d])
+            continue;
+
+          Real stem_offset = edge_stems[d]->relative_coordinate (common[Y_AXIS], Y_AXIS)
+                             - beams[i]->relative_coordinate (common[Y_AXIS], Y_AXIS);
+          Interval heads = Stem::head_positions (edge_stems[d]) * 0.5 * staff_space_;
+
+          Direction ed = edge_dirs_[d];
+          heads.widen (0.5 * staff_space_
+                       + (edge_beam_counts_[d] - 1) * beam_translation_ + beam_thickness_ * .5);
+          quant_range_[d][-ed] = heads[ed] + stem_offset;
+        }
+
+      segments_ = Beam::get_beam_segments (beams[i]);
+      vector_sort (segments_, beam_segment_less);
+      for (vsize j = 0; j < segments_.size (); j++)
+        segments_[j].horizontal_ += (x_span_ - x_pos[LEFT]);
+
+      set<Grob *> colliding_stems;
+      for (vsize j = 0; j < collisions.size (); j++)
+        {
+          if (!collisions[j]->is_live ())
+            continue;
+
+          if (has_interface<Beam> (collisions[j]) && Beam::is_cross_staff (collisions[j]))
+            continue;
+
+          Box b;
+          for (Axis a = X_AXIS; a < NO_AXES; incr (a))
+            b[a] = collisions[j]->extent (common[a], a);
+
+          if (b[X_AXIS][RIGHT] < x_pos[LEFT] || b[X_AXIS][LEFT] > x_pos[RIGHT])
+            continue;
+          if (b[X_AXIS].is_empty () || b[Y_AXIS].is_empty ())
+            continue;
+
+          b[X_AXIS] += (x_span_ - x_pos[LEFT]);
+          b[Y_AXIS] -= my_y;
+          Real width = b[X_AXIS].length ();
+          Real width_factor = sqrt (width / staff_space_);
+
+          for (LEFT_and_RIGHT (d))
+            add_collision (b[X_AXIS][d], b[Y_AXIS], width_factor);
+
+          Grob *stem = unsmob<Grob> (collisions[j]->get_object ("stem"));
+          if (has_interface<Stem> (stem) && Stem::is_normal_stem (stem))
+            {
+              colliding_stems.insert (stem);
+            }
+        }
+
+      for (set<Grob *>::const_iterator it (colliding_stems.begin ()); it != colliding_stems.end (); it++)
+        {
+          Grob *s = *it;
+          Real x = (robust_relative_extent (s, common[X_AXIS], X_AXIS) - x_pos[LEFT] + x_span_).center ();
+
+          Direction stem_dir = get_grob_direction (*it);
+          Interval y;
+          y.set_full ();
+          y[-stem_dir] = Stem::chord_start_y (*it) + (*it)->relative_coordinate (common[Y_AXIS], Y_AXIS)
+                         - my_y;
+
+          Real factor = parameters_.STEM_COLLISION_FACTOR;
+          if (!unsmob<Grob> (s->get_object ("beam")))
+            factor = 1.0;
+          add_collision (x, y, factor);
+        }
+      x_span_ += beams[i]->spanner_length ();
     }
+}
+
+Beam_scoring_problem::Beam_scoring_problem (Grob *me, Drul_array<Real> ys, bool align_broken_intos)
+{
+  beam_ = dynamic_cast<Spanner *> (me);
+  unquanted_y_ = ys;
+  align_broken_intos_ = align_broken_intos;
+
+  parameters_.fill (me);
+  init_instance_variables (me, ys, align_broken_intos);
+  if (do_initial_slope_calculations_)
+    {
+      least_squares_positions ();
+      slope_damping ();
+      shift_region_to_valid ();
+    }
+}
 
-  for (set<Grob *>::const_iterator it (stems.begin ()); it != stems.end (); it++)
+// Assuming V is not empty, pick a 'reasonable' point inside V.
+static Real
+point_in_interval (Interval v, Real dist)
+{
+  if (isinf (v[DOWN]))
+    return v[UP] - dist;
+  else if (isinf (v[UP]))
+    return v[DOWN] + dist;
+  else
+    return v.center ();
+}
+
+/* Set stem's shorten property if unset.
+
+TODO:
+take some y-position (chord/beam/nearest?) into account
+scmify forced-fraction
+
+This is done in beam because the shorten has to be uniform over the
+entire beam.
+*/
+
+void
+set_minimum_dy (Grob *me, Real *dy)
+{
+  if (*dy)
     {
-      Grob *s = *it;
-      Real x = s->extent (common[X_AXIS], X_AXIS).center ();
+      /*
+        If dy is smaller than the smallest quant, we
+        get absurd direction-sign penalties.
+      */
 
-      Direction stem_dir = get_grob_direction (*it);
-      Interval y;
-      y.set_full ();
-      y[-stem_dir] = Stem::chord_start_y (*it) + (*it)->relative_coordinate (common[Y_AXIS], Y_AXIS)
-                     - beam->relative_coordinate (common[Y_AXIS], Y_AXIS);
+      Real ss = Staff_symbol_referencer::staff_space (me);
+      Real beam_thickness = Beam::get_beam_thickness (me) / ss;
+      Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
+      Real sit = (beam_thickness - slt) / 2;
+      Real inter = 0.5;
+      Real hang = 1.0 - (beam_thickness - slt) / 2;
 
-      Real factor = parameters.STEM_COLLISION_FACTOR;
-      if (!unsmob_grob (s->get_object ("beam")))
-        factor = 1.0;
-      add_collision (x, y, factor);
+      *dy = sign (*dy) * max (fabs (*dy),
+                              min (min (sit, inter), hang));
     }
 }
 
-void Beam_scoring_problem::init_stems ()
+void
+Beam_scoring_problem::no_visible_stem_positions ()
 {
-  extract_grob_set (beam, "covered-grobs", collisions);
-  extract_grob_set (beam, "stems", stems);
-  for (int a = 2; a--;)
+  if (!head_positions_.size ())
     {
-      common[a] = common_refpoint_of_array (stems, beam, Axis (a));
-      common[a] = common_refpoint_of_array (collisions, common[a], Axis (a));
+      unquanted_y_ = Interval (0, 0);
+      return;
+    }
+
+  Interval head_positions;
+  Slice multiplicity;
+  for (vsize i = 0; i < head_positions_.size (); i++)
+    {
+      head_positions.unite (head_positions_[i]);
+      multiplicity.unite (beam_multiplicity_[i]);
     }
 
-  Drul_array<Grob *> edge_stems (Beam::first_normal_stem (beam),
-                                 Beam::last_normal_stem (beam));
-  Direction d = LEFT;
-  do
-    x_span[d] = edge_stems[d] ? edge_stems[d]->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
-  while (flip (&d) != LEFT);
+  Direction dir = get_grob_direction (beam_);
+
+  if (!dir)
+    programming_error ("The beam should have a direction by now.");
+
+  Real y = head_positions.linear_combination (dir)
+           * 0.5 * staff_space_
+           + dir * beam_translation_ * (multiplicity.length () + 1);
+
+  unquanted_y_ = Interval (y, y);
+}
+
+vsize
+Beam_scoring_problem::first_normal_index ()
+{
+  for (vsize i = 0; i < is_normal_.size (); i++)
+    if (is_normal_[i])
+      return i;
+
+  beam_->programming_error ("No normal stems, but asking for first normal stem index.");
+  return 0;
+}
+
+vsize
+Beam_scoring_problem::last_normal_index ()
+{
+  for (vsize i = is_normal_.size (); i--;)
+    if (is_normal_[i])
+      return i;
 
-  Drul_array<bool> dirs_found (0, 0);
-  for (vsize i = 0; i < stems.size (); i++)
+  beam_->programming_error ("No normal stems, but asking for first normal stem index.");
+  return 0;
+}
+
+void
+Beam_scoring_problem::least_squares_positions ()
+{
+  if (!normal_stem_count_)
     {
-      Grob *s = stems[i];
-      if (!Stem::is_normal_stem (s))
-        continue;
+      no_visible_stem_positions ();
+      return;
+    }
+
+  if (stem_infos_.size () < 1)
+    return;
+
+  vsize fnx = first_normal_index ();
+  vsize lnx = last_normal_index ();
+
+  Interval ideal (stem_infos_[fnx].ideal_y_ + stem_ypositions_[fnx],
+                  stem_infos_[lnx].ideal_y_ + stem_ypositions_[lnx]);
+
+  Real y = 0;
+  Real slope = 0;
+  Real dy = 0;
+  Real ldy = 0.0;
+  if (!ideal.delta ())
+    {
+      Interval chord (chord_start_y_[0],
+                      chord_start_y_.back ());
+
+      /* Simple beams (2 stems) on middle line should be allowed to be
+         slightly sloped.
+
+         However, if both stems reach middle line,
+         ideal[LEFT] == ideal[RIGHT] and ideal.delta () == 0.
 
-      Stem_info si (Stem::get_stem_info (s));
-      si.scale (1 / staff_space);
-      stem_infos.push_back (si);
-      dirs_found[si.dir_] = true;
+         For that case, we apply artificial slope */
+      if (!ideal[LEFT] && chord.delta () && stem_infos_.size () == 2)
+        {
+          /* FIXME. -> UP */
+          Direction d = (Direction) (sign (chord.delta ()) * UP);
+          unquanted_y_[d] = Beam::get_beam_thickness (beam_) / 2;
+          unquanted_y_[-d] = -unquanted_y_[d];
+        }
+      else
+        unquanted_y_ = ideal;
+
+      ldy = unquanted_y_[RIGHT] - unquanted_y_[LEFT];
+    }
+  else
+    {
+      vector<Offset> ideals;
+      for (vsize i = 0; i < stem_infos_.size (); i++)
+        if (is_normal_[i])
+          ideals.push_back (Offset (stem_xpositions_[i],
+                                    stem_infos_[i].ideal_y_
+                                    + stem_ypositions_[i]));
 
-      bool f = to_boolean (s->get_property ("french-beaming"))
-               && s != edge_stems[LEFT] && s != edge_stems[RIGHT];
+      minimise_least_squares (&slope, &y, ideals);
 
-      Real y = Beam::calc_stem_y (beam, s, common, x_span[LEFT], x_span[RIGHT], CENTER,
-                                  Interval (0, 0), f);
-      base_lengths.push_back (y / staff_space);
-      stem_xpositions.push_back (s->relative_coordinate (common[X_AXIS], X_AXIS));
+      dy = slope * x_span_;
+
+      set_minimum_dy (beam_, &dy);
+
+      ldy = dy;
+      unquanted_y_ = Interval (y, (y + dy));
     }
 
-  edge_dirs = Drul_array<Direction> (CENTER, CENTER);
-  if (stem_infos.size ())
+  musical_dy_ = ldy;
+  beam_->set_property ("least-squares-dy", scm_from_double (musical_dy_));
+}
+
+/*
+  Determine whether a beam is concave.
+
+  A beam is concave when the middle notes get closer to the
+  beam than the left and right edge notes.
+
+  This is determined in two ways: by looking at the positions of the
+  middle notes, or by looking at the deviation of the inside notes
+  compared to the line connecting first and last.
+
+  The tricky thing is what to do with beams with chords. There are no
+  real guidelines in this case.
+*/
+
+bool
+is_concave_single_notes (vector<int> const &positions, Direction beam_dir)
+{
+  Interval covering;
+  covering.add_point (positions[0]);
+  covering.add_point (positions.back ());
+
+  bool above = false;
+  bool below = false;
+  bool concave = false;
+
+  /*
+    notes above and below the interval covered by 1st and last note.
+  */
+  for (vsize i = 1; i + 1 < positions.size (); i++)
     {
-      edge_dirs = Drul_array<Direction> (stem_infos[0].dir_,
-                                         stem_infos.back ().dir_);
+      above = above || (positions[i] > covering[UP]);
+      below = below || (positions[i] < covering[DOWN]);
     }
 
-  is_xstaff = Align_interface::has_interface (common[Y_AXIS]);
-  is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
+  concave = concave || (above && below);
+  /*
+    A note as close or closer to the beam than begin and end, but the
+    note is reached in the opposite direction as the last-first dy
+  */
+  int dy = positions.back () - positions[0];
+  int closest = max (beam_dir * positions.back (), beam_dir * positions[0]);
+  for (vsize i = 2; !concave && i + 1 < positions.size (); i++)
+    {
+      int inner_dy = positions[i] - positions[i - 1];
+      if (sign (inner_dy) != sign (dy)
+          && (beam_dir * positions[i] >= closest
+              || beam_dir * positions[i - 1] >= closest))
+        concave = true;
+    }
 
-  staff_radius = Staff_symbol_referencer::staff_radius (beam);
-  edge_beam_counts = Drul_array<int>
-                     (Stem::beam_multiplicity (stems[0]).length () + 1,
-                      Stem::beam_multiplicity (stems.back ()).length () + 1);
+  bool all_closer = true;
+  for (vsize i = 1; all_closer && i + 1 < positions.size (); i++)
+    {
+      all_closer = all_closer
+                   && (beam_dir * positions[i] > closest);
+    }
 
-  // TODO - why are we dividing by staff_space?
-  beam_translation = Beam::get_beam_translation (beam) / staff_space;
+  concave = concave || all_closer;
+  return concave;
+}
 
-  d = LEFT;
-  do
+Real
+calc_positions_concaveness (vector<int> const &positions, Direction beam_dir)
+{
+  Real dy = positions.back () - positions[0];
+  Real slope = dy / Real (positions.size () - 1);
+  Real concaveness = 0.0;
+  for (vsize i = 1; i + 1 < positions.size (); i++)
     {
-      quant_range[d].set_full ();
-      if (!edge_stems[d])
-        continue;
+      Real line_y = slope * i + positions[0];
+
+      concaveness += max (beam_dir * (positions[i] - line_y), 0.0);
+    }
+
+  concaveness /= positions.size ();
+
+  /*
+    Normalize. For dy = 0, the slope ends up as 0 anyway, so the
+    scaling of concaveness doesn't matter much.
+  */
+  if (dy)
+    concaveness /= fabs (dy);
+  return concaveness;
+}
+
+Real
+Beam_scoring_problem::calc_concaveness ()
+{
+  SCM conc = beam_->get_property ("concaveness");
+  if (scm_is_number (conc))
+    return scm_to_double (conc);
+
+  if (is_knee_ || is_xstaff_)
+    return 0.0;
 
-      Real stem_offset = edge_stems[d]->relative_coordinate (common[Y_AXIS], Y_AXIS)
-                         - beam->relative_coordinate (common[Y_AXIS], Y_AXIS);
-      Interval heads = Stem::head_positions (edge_stems[d]) * 0.5 * staff_space;
+  Direction beam_dir = CENTER;
+  for (vsize i = is_normal_.size (); i--;)
+    if (is_normal_[i] && stem_infos_[i].dir_)
+      beam_dir = stem_infos_[i].dir_;
 
-      Direction ed = edge_dirs[d];
-      heads.widen (0.5 * staff_space
-                   + (edge_beam_counts[d] - 1) * beam_translation + beam_thickness * .5);
-      quant_range[d][-ed] = heads[ed] + stem_offset;
+  if (normal_stem_count_ <= 2)
+    return 0.0;
+
+  vector<int> close_positions;
+  vector<int> far_positions;
+  for (vsize i = 0; i < is_normal_.size (); i++)
+    if (is_normal_[i])
+      {
+        /*
+          For chords, we take the note head that is closest to the beam.
+
+          Hmmm.. wait, for the beams in the last measure of morgenlied,
+          this doesn't look so good. Let's try the heads farthest from
+          the beam.
+        */
+
+        close_positions.push_back ((int) rint (head_positions_[i][beam_dir]));
+        far_positions.push_back ((int) rint (head_positions_[i][-beam_dir]));
+      }
+
+  Real concaveness = 0.0;
+
+  if (is_concave_single_notes (beam_dir == UP ? close_positions : far_positions, beam_dir))
+    {
+      concaveness = 10000;
+    }
+  else
+    {
+      concaveness = (calc_positions_concaveness (far_positions, beam_dir)
+                     + calc_positions_concaveness (close_positions, beam_dir)) / 2;
+    }
+
+  return concaveness;
+}
+
+void
+Beam_scoring_problem::slope_damping ()
+{
+  if (normal_stem_count_ <= 1)
+    return;
+
+  SCM s = beam_->get_property ("damping");
+  Real damping = scm_to_double (s);
+  Real concaveness = calc_concaveness ();
+  if (concaveness >= 10000)
+    {
+      unquanted_y_[LEFT] = unquanted_y_[RIGHT];
+      musical_dy_ = 0;
+      damping = 0;
     }
-  while (flip (&d) != LEFT);
 
-  init_collisions (collisions);
+  if (damping)
+    {
+      Real dy = unquanted_y_[RIGHT] - unquanted_y_[LEFT];
+
+      Real slope = dy && x_span_ ? dy / x_span_ : 0;
+
+      slope = 0.6 * tanh (slope) / (damping + concaveness);
+
+      Real damped_dy = slope * x_span_;
+
+      set_minimum_dy (beam_, &damped_dy);
+
+      unquanted_y_[LEFT] += (dy - damped_dy) / 2;
+      unquanted_y_[RIGHT] -= (dy - damped_dy) / 2;
+    }
 }
 
-Beam_scoring_problem::Beam_scoring_problem (Grob *me, Drul_array<Real> ys)
+void
+Beam_scoring_problem::shift_region_to_valid ()
 {
-  beam = me;
-  unquanted_y = ys;
+  if (!normal_stem_count_)
+    return;
+
+  Real beam_dy = unquanted_y_[RIGHT] - unquanted_y_[LEFT];
+  Real slope = x_span_ ? beam_dy / x_span_ : 0.0;
 
   /*
-    Calculations are relative to a unit-scaled staff, i.e. the quants are
-    divided by the current staff_space.
+    Shift the positions so that we have a chance of finding good
+    quants (i.e. no short stem failures.)
   */
-  staff_space = Staff_symbol_referencer::staff_space (me);
-  beam_thickness = Beam::get_beam_thickness (me) / staff_space;
-  line_thickness = Staff_symbol_referencer::line_thickness (me) / staff_space;
+  Interval feasible_left_point;
+  feasible_left_point.set_full ();
 
-  // This is the least-squares DY, corrected for concave beams.
-  musical_dy = robust_scm2double (me->get_property ("least-squares-dy"), 0);
+  for (vsize i = 0; i < stem_infos_.size (); i++)
+    {
+      // TODO - check for invisible here...
+      Real left_y
+        = stem_infos_[i].shortest_y_
+          - slope * stem_xpositions_ [i];
+
+      /*
+        left_y is now relative to the stem S. We want relative to
+        ourselves, so translate:
+      */
+      left_y += stem_ypositions_[i];
+      Interval flp;
+      flp.set_full ();
+      flp[-stem_infos_[i].dir_] = left_y;
+
+      feasible_left_point.intersect (flp);
+    }
+
+  vector<Grob *> filtered;
+  /*
+    We only update these for objects that are too large for quanting
+    to find a workaround.  Typically, these are notes with
+    stems, and timesig/keysig/clef, which take out the entire area
+    inside the staff as feasible.
+
+    The code below disregards the thickness and multiplicity of the
+    beam.  This should not be a problem, as the beam quanting will
+    take care of computing the impact those exactly.
+  */
+  Real min_y_size = 2.0;
+
+  // A list of intervals into which beams may not fall
+  vector<Interval> forbidden_intervals;
+
+  for (vsize i = 0; i < collisions_.size (); i++)
+    {
+      if (collisions_[i].x_ < 0 || collisions_[i].x_ > x_span_)
+        continue;
+
+      if (collisions_[i].y_.length () < min_y_size)
+        continue;
+
+      for (LEFT_and_RIGHT (d))
+        {
+          Real dy = slope * collisions_[i].x_;
+
+          Interval disallowed;
+          for (DOWN_and_UP (yd))
+            {
+              Real left_y = collisions_[i].y_[yd] - dy;
+              disallowed[yd] = left_y;
+            }
+
+          forbidden_intervals.push_back (disallowed);
+        }
+    }
+
+  vector_sort (forbidden_intervals, Interval::left_less);
+  Real beam_left_y = unquanted_y_[LEFT];
+  Interval feasible_beam_placements (beam_left_y, beam_left_y);
+
+  Interval_minefield minefield (feasible_beam_placements, 0.0);
+  for (vsize i = 0; i < forbidden_intervals.size (); i++)
+    minefield.add_forbidden_interval (forbidden_intervals[i]);
+  minefield.solve ();
+  feasible_beam_placements = minefield.feasible_placements ();
 
-  parameters.fill (me);
-  init_stems ();
+  // if the beam placement falls out of the feasible region, we push it
+  // to infinity so that it can never be a feasible candidate below
+  for (DOWN_and_UP (d))
+    {
+      if (!feasible_left_point.contains (feasible_beam_placements[d]))
+        feasible_beam_placements[d] = d * infinity_f;
+    }
+
+  if ((feasible_beam_placements[UP] == infinity_f && feasible_beam_placements[DOWN] == -infinity_f) && !feasible_left_point.is_empty ())
+    {
+      // We are somewhat screwed: we have a collision, but at least
+      // there is a way to satisfy stem length constraints.
+      beam_left_y = point_in_interval (feasible_left_point, 2.0);
+    }
+  else if (!feasible_left_point.is_empty ())
+    {
+      // Only one of them offers is feasible solution. Pick that one.
+      if (abs (beam_left_y - feasible_beam_placements[DOWN]) > abs (beam_left_y - feasible_beam_placements[UP]))
+        beam_left_y = feasible_beam_placements[UP];
+      else
+        beam_left_y = feasible_beam_placements[DOWN];
+    }
+  else
+    {
+      // We are completely screwed.
+      beam_->warning (_ ("no viable initial configuration found: may not find good beam slope"));
+    }
+
+  unquanted_y_ = Drul_array<Real> (beam_left_y, (beam_left_y + beam_dy));
 }
 
 void
 Beam_scoring_problem::generate_quants (vector<Beam_configuration *> *scores) const
 {
-  int region_size = (int) parameters.REGION_SIZE;
+  int region_size = (int) parameters_.REGION_SIZE;
 
   // Knees and collisions are harder, lets try some more possibilities
-  if (is_knee)
+  if (is_knee_)
     region_size += 2;
   if (collisions_.size ())
     region_size += 2;
 
   Real straddle = 0.0;
-  Real sit = (beam_thickness - line_thickness) / 2;
+  Real sit = (beam_thickness_ - line_thickness_) / 2;
   Real inter = 0.5;
-  Real hang = 1.0 - (beam_thickness - line_thickness) / 2;
+  Real hang = 1.0 - (beam_thickness_ - line_thickness_) / 2;
   Real base_quants [] = {straddle, sit, inter, hang};
   int num_base_quants = int (sizeof (base_quants) / sizeof (Real));
 
@@ -367,21 +885,19 @@ Beam_scoring_problem::generate_quants (vector<Beam_configuration *> *scores) con
     for (vsize j = 0; j < unshifted_quants.size (); j++)
       {
         Beam_configuration *c
-          = Beam_configuration::new_config (unquanted_y,
+          = Beam_configuration::new_config (unquanted_y_,
                                             Interval (unshifted_quants[i],
                                                       unshifted_quants[j]));
 
-        Direction d = LEFT;
-        do
+        for (LEFT_and_RIGHT (d))
           {
-            if (!quant_range[d].contains (c->y[d]))
+            if (!quant_range_[d].contains (c->y[d]))
               {
                 delete c;
                 c = NULL;
                 break;
               }
           }
-        while (flip (&d) != LEFT);
         if (c)
           scores->push_back (c);
       }
@@ -456,14 +972,17 @@ Beam_scoring_problem::solve () const
   if (configs.empty ())
     {
       programming_error ("No viable beam quanting found.  Using unquanted y value.");
-      return unquanted_y;
+      return unquanted_y_;
     }
 
+  if (to_boolean (beam_->get_property ("skip-quanting")))
+    return unquanted_y_;
+
   Beam_configuration *best = NULL;
 
   bool debug
-    = to_boolean (beam->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring")));
-  SCM inspect_quants = beam->get_property ("inspect-quants");
+    = to_boolean (beam_->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring")));
+  SCM inspect_quants = beam_->get_property ("inspect-quants");
   if (scm_is_pair (inspect_quants))
     {
       debug = true;
@@ -518,32 +1037,44 @@ Beam_scoring_problem::solve () const
         }
 
       string card = best->score_card_ + to_string (" c%d/%d", completed, configs.size ());
-      beam->set_property ("annotation", ly_string2scm (card));
+      beam_->set_property ("annotation", ly_string2scm (card));
     }
 #endif
 
   junk_pointers (configs);
+  if (align_broken_intos_)
+    {
+      Interval normalized_endpoints = robust_scm2interval (beam_->get_property ("normalized-endpoints"), Interval (0, 1));
+      Real y_length = final_positions[RIGHT] - final_positions[LEFT];
+
+      final_positions[LEFT] += normalized_endpoints[LEFT] * y_length;
+      final_positions[RIGHT] -= (1 - normalized_endpoints[RIGHT]) * y_length;
+    }
+
   return final_positions;
 }
 
 void
 Beam_scoring_problem::score_stem_lengths (Beam_configuration *config) const
 {
-  Real limit_penalty = parameters.STEM_LENGTH_LIMIT_PENALTY;
+  Real limit_penalty = parameters_.STEM_LENGTH_LIMIT_PENALTY;
   Drul_array<Real> score (0, 0);
   Drul_array<int> count (0, 0);
 
-  for (vsize i = 0; i < stem_xpositions.size (); i++)
+  for (vsize i = 0; i < stem_xpositions_.size (); i++)
     {
-      Real x = stem_xpositions[i];
-      Real dx = x_span.delta ();
+      if (!is_normal_[i])
+        continue;
+
+      Real x = stem_xpositions_[i];
+      Real dx = x_span_;
       Real beam_y = dx
-                    ? config->y[RIGHT] * (x - x_span[LEFT]) / dx + config->y[LEFT] * (x_span[RIGHT] - x) / dx
+                    ? config->y[RIGHT] * x / dx + config->y[LEFT] * (x_span_ - x) / dx
                     : (config->y[RIGHT] + config->y[LEFT]) / 2;
-      Real current_y = beam_y + base_lengths[i];
-      Real length_pen = parameters.STEM_LENGTH_DEMERIT_FACTOR;
+      Real current_y = beam_y + base_lengths_[i];
+      Real length_pen = parameters_.STEM_LENGTH_DEMERIT_FACTOR;
 
-      Stem_info info = stem_infos[i];
+      Stem_info info = stem_infos_[i];
       Direction d = info.dir_;
 
       score[d] += limit_penalty * max (0.0, (d * (info.shortest_y_ - current_y)));
@@ -554,7 +1085,7 @@ Beam_scoring_problem::score_stem_lengths (Beam_configuration *config) const
       /* We introduce a power, to make the scoring strictly
          convex. Otherwise a symmetric knee beam (up/down/up/down)
          does not have an optimum in the middle. */
-      if (is_knee)
+      if (is_knee_)
         ideal_score = pow (ideal_score, 1.1);
 
       score[d] += length_pen * ideal_score;
@@ -562,10 +1093,18 @@ Beam_scoring_problem::score_stem_lengths (Beam_configuration *config) const
     }
 
   /* Divide by number of stems, to make the measure scale-free. */
-  Direction d = DOWN;
-  do
+  for (DOWN_and_UP (d))
     score[d] /= max (count[d], 1);
-  while (flip (&d) != DOWN);
+
+  /*
+    sometimes, two perfectly symmetric kneed beams will have the same score
+    and can either be quanted up or down.
+
+    we choose the quanting in the direction of the slope so that the first stem
+    always seems longer, reaching to the second, rather than squashed.
+  */
+  if (is_knee_ && count[LEFT] == count[RIGHT] && count[LEFT] == 1 && unquanted_y_.delta ())
+    score[Direction (sign (unquanted_y_.delta ()))] += score[Direction (sign (unquanted_y_.delta ()))] < 1.0 ? 0.01 : 0.0;
 
   config->add (score[LEFT] + score[RIGHT], "L");
 }
@@ -574,7 +1113,7 @@ void
 Beam_scoring_problem::score_slope_direction (Beam_configuration *config) const
 {
   Real dy = config->y.delta ();
-  Real damped_dy = unquanted_y.delta ();
+  Real damped_dy = unquanted_y_.delta ();
   Real dem = 0.0;
   /*
     DAMPING_DIRECTION_PENALTY is a very harsh measure, while for
@@ -587,13 +1126,13 @@ Beam_scoring_problem::score_slope_direction (Beam_configuration *config) const
     {
       if (!dy)
         {
-          if (fabs (damped_dy / x_span.delta ()) > parameters.ROUND_TO_ZERO_SLOPE)
-            dem += parameters.DAMPING_DIRECTION_PENALTY;
+          if (fabs (damped_dy / x_span_) > parameters_.ROUND_TO_ZERO_SLOPE)
+            dem += parameters_.DAMPING_DIRECTION_PENALTY;
           else
-            dem += parameters.HINT_DIRECTION_PENALTY;
+            dem += parameters_.HINT_DIRECTION_PENALTY;
         }
       else
-        dem += parameters.DAMPING_DIRECTION_PENALTY;
+        dem += parameters_.DAMPING_DIRECTION_PENALTY;
     }
 
   config->add (dem, "Sd");
@@ -604,8 +1143,8 @@ void
 Beam_scoring_problem::score_slope_musical (Beam_configuration *config) const
 {
   Real dy = config->y.delta ();
-  Real dem = parameters.MUSICAL_DIRECTION_FACTOR
-             * max (0.0, (fabs (dy) - fabs (musical_dy)));
+  Real dem = parameters_.MUSICAL_DIRECTION_FACTOR
+             * max (0.0, (fabs (dy) - fabs (musical_dy_)));
   config->add (dem, "Sm");
 }
 
@@ -614,14 +1153,14 @@ void
 Beam_scoring_problem::score_slope_ideal (Beam_configuration *config) const
 {
   Real dy = config->y.delta ();
-  Real damped_dy = unquanted_y.delta ();
+  Real damped_dy = unquanted_y_.delta ();
   Real dem = 0.0;
 
-  Real slope_penalty = parameters.IDEAL_SLOPE_FACTOR;
+  Real slope_penalty = parameters_.IDEAL_SLOPE_FACTOR;
 
   /* Xstaff beams tend to use extreme slopes to get short stems. We
      put in a penalty here. */
-  if (is_xstaff)
+  if (is_xstaff_)
     slope_penalty *= 10;
 
   /* Huh, why would a too steep beam be better than a too flat one ? */
@@ -644,11 +1183,11 @@ void
 Beam_scoring_problem::score_horizontal_inter_quants (Beam_configuration *config) const
 {
   if (config->y.delta () == 0.0
-      && abs (config->y[LEFT]) < staff_radius * staff_space)
+      && abs (config->y[LEFT]) < staff_radius_ * staff_space_)
     {
-      Real yshift = config->y[LEFT] - 0.5 * staff_space;
-      if (fabs (my_round (yshift) - yshift) < 0.01 * staff_space)
-        config->add (parameters.HORIZONTAL_INTER_QUANT_PENALTY, "H");
+      Real yshift = config->y[LEFT] - 0.5 * staff_space_;
+      if (fabs (my_round (yshift) - yshift) < 0.01 * staff_space_)
+        config->add (parameters_.HORIZONTAL_INTER_QUANT_PENALTY, "H");
     }
 }
 
@@ -662,42 +1201,51 @@ Beam_scoring_problem::score_forbidden_quants (Beam_configuration *config) const
 {
   Real dy = config->y.delta ();
 
-  Real extra_demerit = parameters.SECONDARY_BEAM_DEMERIT
-                       / max (edge_beam_counts[LEFT], edge_beam_counts[RIGHT]);
-
-  Direction d = LEFT;
+  Real extra_demerit =
+    parameters_.SECONDARY_BEAM_DEMERIT
+    / max (edge_beam_counts_[LEFT], edge_beam_counts_[RIGHT]);
+  
   Real dem = 0.0;
-  Real eps = parameters.BEAM_EPS;
+  Real eps = parameters_.BEAM_EPS;
 
-  do
+  for (LEFT_and_RIGHT (d))
     {
-      for (int j = 1; j <= edge_beam_counts[d]; j++)
+      for (int j = 1; j <= edge_beam_counts_[d]; j++)
         {
-          Direction stem_dir = edge_dirs[d];
+          Direction stem_dir = edge_dirs_[d];
 
           /*
-            The 2.2 factor is to provide a little leniency for
+            The fudge_factor is to provide a little leniency for
             borderline cases. If we do 2.0, then the upper outer line
             will be in the gap of the (2, sit) quant, leading to a
-            false demerit.
+            false demerit. By increasing the fudge factor to 2.2, we
+            fix this case.
           */
-          Real gap1 = config->y[d] - stem_dir * ((j - 1) * beam_translation + beam_thickness / 2 - line_thickness / 2.2);
-          Real gap2 = config->y[d] - stem_dir * (j * beam_translation - beam_thickness / 2 + line_thickness / 2.2);
+          Real fudge_factor = 2.2;
+          Real gap1 = config->y[d] - stem_dir * ((j - 1) * beam_translation_ + beam_thickness_ / 2 - line_thickness_ / fudge_factor);
+          Real gap2 = config->y[d] - stem_dir * (j * beam_translation_ - beam_thickness_ / 2 + line_thickness_ / fudge_factor);
 
           Interval gap;
           gap.add_point (gap1);
           gap.add_point (gap2);
 
-          for (Real k = -staff_radius;
-               k <= staff_radius + eps; k += 1.0)
+          for (Real k = -staff_radius_;
+               k <= staff_radius_ + eps; k += 1.0)
             if (gap.contains (k))
               {
                 Real dist = min (fabs (gap[UP] - k), fabs (gap[DOWN] - k));
 
                 /*
                   this parameter is tuned to grace-stem-length.ly
+                  retuned from 0.40 to 0.39 by MS because of slight increases
+                  in gap.length () resulting from measuring beams at real ends
+                  instead of from the middle of stems.
+
+                  TODO:
+                  This function needs better comments so we know what is forbidden
+                  and why.
                 */
-                Real fixed_demerit = 0.4;
+                Real fixed_demerit = 0.39;
 
                 dem += extra_demerit
                        * (fixed_demerit
@@ -705,48 +1253,47 @@ Beam_scoring_problem::score_forbidden_quants (Beam_configuration *config) const
               }
         }
     }
-  while ((flip (&d)) != LEFT);
 
-  if (max (edge_beam_counts[LEFT], edge_beam_counts[RIGHT]) >= 2)
+  config->add (dem, "Fl");
+  dem = 0.0;
+  if (max (edge_beam_counts_[LEFT], edge_beam_counts_[RIGHT]) >= 2)
     {
       Real straddle = 0.0;
-      Real sit = (beam_thickness - line_thickness) / 2;
+      Real sit = (beam_thickness_ - line_thickness_) / 2;
       Real inter = 0.5;
-      Real hang = 1.0 - (beam_thickness - line_thickness) / 2;
+      Real hang = 1.0 - (beam_thickness_ - line_thickness_) / 2;
 
-      Direction d = LEFT;
-      do
+      for (LEFT_and_RIGHT (d))
         {
-          if (edge_beam_counts[d] >= 2
-              && fabs (config->y[d] - edge_dirs[d] * beam_translation) < staff_radius + inter)
+          if (edge_beam_counts_[d] >= 2
+              && fabs (config->y[d] - edge_dirs_[d] * beam_translation_) < staff_radius_ + inter)
             {
               // TODO up/down symmetry.
-              if (edge_dirs[d] == UP && dy <= eps
+              if (edge_dirs_[d] == UP && dy <= eps
                   && fabs (my_modf (config->y[d]) - sit) < eps)
                 dem += extra_demerit;
 
-              if (edge_dirs[d] == DOWN && dy >= eps
+              if (edge_dirs_[d] == DOWN && dy >= eps
                   && fabs (my_modf (config->y[d]) - hang) < eps)
                 dem += extra_demerit;
             }
 
-          if (edge_beam_counts[d] >= 3
-              && fabs (config->y[d] - 2 * edge_dirs[d] * beam_translation) < staff_radius + inter)
+          if (edge_beam_counts_[d] >= 3
+              && fabs (config->y[d] - 2 * edge_dirs_[d] * beam_translation_) < staff_radius_ + inter)
             {
               // TODO up/down symmetry.
-              if (edge_dirs[d] == UP && dy <= eps
+              if (edge_dirs_[d] == UP && dy <= eps
                   && fabs (my_modf (config->y[d]) - straddle) < eps)
                 dem += extra_demerit;
 
-              if (edge_dirs[d] == DOWN && dy >= eps
+              if (edge_dirs_[d] == DOWN && dy >= eps
                   && fabs (my_modf (config->y[d]) - straddle) < eps)
                 dem += extra_demerit;
             }
         }
-      while (flip (&d) != LEFT);
     }
 
-  config->add (dem, "F");
+  config->add (dem, "Fs");
 }
 
 void
@@ -764,16 +1311,20 @@ Beam_scoring_problem::score_collisions (Beam_configuration *config) const
       Real dist = infinity_f;
       if (!intersection (beam_y, collision_y).is_empty ())
         dist = 0.0;
-      else
+      else 
         dist = min (beam_y.distance (collision_y[DOWN]),
                     beam_y.distance (collision_y[UP]));
 
+      
       Real scale_free
-        = max (parameters.COLLISION_PADDING - dist, 0.0) /
-          parameters.COLLISION_PADDING;
-      demerits
-      += collisions_[i].base_penalty_ *
-         pow (scale_free, 3) * parameters.COLLISION_PENALTY;
+        = max (parameters_.COLLISION_PADDING - dist, 0.0)
+          / parameters_.COLLISION_PADDING;
+      Real collision_demerit = collisions_[i].base_penalty_ *
+         pow (scale_free, 3) * parameters_.COLLISION_PENALTY;
+
+      if (collision_demerit > 0) {
+        demerits += collision_demerit;
+      }
     }
 
   config->add (demerits, "C");