]> git.donarmstrong.com Git - lilypond.git/blobdiff - lily/beam-quanting.cc
Rename collision-distance beam detail property to collision-padding.
[lilypond.git] / lily / beam-quanting.cc
index f358065b8d65331eb8d5bb6a108016b4fec962d3..6b495b379427ac4d7f0f383058162a53a61a8a92 100644 (file)
 /*
-  beam-quanting.cc -- implement Beam quanting functions
-  
-  source file of the GNU LilyPond music typesetter
-  
-  (c) 1997--2005 Han-Wen Nienhuys <hanwen@cs.uu.nl>
+  This file is part of LilyPond, the GNU music typesetter.
+
+  Copyright (C) 1997--2011 Han-Wen Nienhuys <hanwen@xs4all.nl>
   Jan Nieuwenhuizen <janneke@gnu.org>
-  
+
+  LilyPond is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  LilyPond is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with LilyPond.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <math.h>
+#include "beam-scoring-problem.hh"
 
-#include "warn.hh"
-#include "staff-symbol-referencer.hh"
+#include <queue>  
+#include <set>
+#include <algorithm>
+using namespace std;
+
+#include "align-interface.hh"
 #include "beam.hh"
-#include "stem.hh"
+#include "direction.hh"
+#include "directional-element-interface.hh"
+#include "grob.hh"
+#include "international.hh"
+#include "main.hh"
 #include "output-def.hh"
-#include "group-interface.hh"
-#include "align-interface.hh"
+#include "pointer-group-interface.hh"
+#include "staff-symbol-referencer.hh"
+#include "stencil.hh"
+#include "stem.hh"
+#include "warn.hh"
 
-const int INTER_QUANT_PENALTY = 1000;
-const Real SECONDARY_BEAM_DEMERIT  = 10.0;
-const int STEM_LENGTH_DEMERIT_FACTOR = 5;
+Real
+get_detail (SCM alist, SCM sym, Real def)
+{
+  SCM entry = scm_assq (sym, alist);
 
-/*
-  threshold to combat rounding errors.
- */
-const Real BEAM_EPS = 1e-3; 
+  if (scm_is_pair (entry))
+    return robust_scm2double (scm_cdr (entry), def);
+  return def;
+}
 
-// possibly ridiculous, but too short stems just won't do
-const int STEM_LENGTH_LIMIT_PENALTY = 5000;
-const int DAMPING_DIRECTION_PENALTY = 800;
-const int MUSICAL_DIRECTION_FACTOR = 400;
-const int IDEAL_SLOPE_FACTOR = 10;
-const Real ROUND_TO_ZERO_SLOPE = 0.02;
+void
+Beam_quant_parameters::fill (Grob *him)
+{
+  SCM details = him->get_property ("details");
+
+  // General
+  BEAM_EPS = get_detail (details, ly_symbol2scm ("beam-eps"), 1e-3);
+  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);
+  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);
+
+  STEM_LENGTH_LIMIT_PENALTY = get_detail (details, ly_symbol2scm ("stem-length-limit-penalty"), 5000);
+  DAMPING_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("damping-direction-penalty"), 800);
+  HINT_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("hint-direction-penalty"), 20);
+  MUSICAL_DIRECTION_FACTOR = get_detail (details, ly_symbol2scm ("musical-direction-factor"), 400);
+  IDEAL_SLOPE_FACTOR = get_detail (details, ly_symbol2scm ("ideal-slope-factor"), 10);
+  ROUND_TO_ZERO_SLOPE = get_detail (details, ly_symbol2scm ("round-to-zero-slope"), 0.02);
+
+  // Collisions
+  COLLISION_PENALTY = get_detail (details, ly_symbol2scm ("collision-penalty"), 500);
+  COLLISION_PADDING = get_detail (details, ly_symbol2scm ("collision-padding"), 0.5);
+  STEM_COLLISION_FACTOR = get_detail (details, ly_symbol2scm ("stem-collision-factor"), 0.1);
+}
 
+// Add x if x is positive, add |x|*fac if x is negative.
 static Real
 shrink_extra_weight (Real x, Real fac)
 {
   return fabs (x) * ((x < 0) ? fac : 1.0);
 }
 
+/****************************************************************/
 
-struct Quant_score
+Beam_configuration::Beam_configuration ()
 {
-  Real yl;
-  Real yr;
-  Real demerits;
+  y = Interval (0.0, 0.0);
+  demerits = 0.0;
+  next_scorer_todo = ORIGINAL_DISTANCE;
+}
 
-#if DEBUG_QUANTING
-  String score_card_;
+bool Beam_configuration::done () const
+{
+  return next_scorer_todo >= NUM_SCORERS;
+}
+
+void Beam_configuration::add (Real demerit, const string &reason)
+{
+  demerits += demerit;
+
+#if DEBUG_BEAM_SCORING
+  if (demerit) 
+    score_card_ += to_string (" %s %.2f", reason.c_str (), demerit);
 #endif
-};
+}
+  
+Beam_configuration* Beam_configuration::new_config (Interval start,
+                                                    Interval offset)
+{
+  Beam_configuration* qs = new Beam_configuration;
+  qs->y = Interval (int (start[LEFT]) + offset[LEFT],
+                    int (start[RIGHT]) + offset[RIGHT]);
+
+  // This orders the sequence so we try combinations closest to the
+  // the ideal offset first.
+  Real start_score = abs (offset[RIGHT]) + abs (offset[LEFT]);
+  qs->demerits = start_score / 1000.0;
+  qs->next_scorer_todo = ORIGINAL_DISTANCE + 1;
+  
+  return qs;
+}
+
+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();
+}
 
+/****************************************************************/
 
 /*
   TODO:
-  
-   - Make all demerits customisable
-
-   - One sensible check per demerit (what's this --hwn)
 
-   - Add demerits for quants per se, as to forbid a specific quant
-     entirely
+  - Make all demerits customisable
 
+  - Add demerits for quants per se, as to forbid a specific quant
+  entirely
 */
 
-int
-best_quant_score_idx (Array<Quant_score>  const & qscores)
+// This is a temporary hack to see how much we can gain by using a
+// priority queue on the beams to score.
+static int score_count = 0;
+LY_DEFINE (ly_beam_score_count, "ly:beam-score-count", 0, 0, 0,
+          (),
+          "count number of beam scores.") {
+  return scm_from_int (score_count);
+}
+
+void Beam_scoring_problem::add_collision (Real x, Interval y,
+                                          Real score_factor)
 {
-  Real best = 1e6;
-  int best_idx = -1;
-  for (int i = qscores.size (); i--;)
-    {
-      if (qscores[i].demerits < best)
-       {
-         best = qscores [i].demerits ;
-         best_idx = i;
-       }
+  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;
     }
+  }
 
-  if (best_idx < 0)
+  Beam_collision c;
+  c.beam_y_.set_empty ();
+
+  for (vsize j = 0; j < segments_.size (); j++)
     {
-      programming_error ("Huh? No best beam quant score?");
-      best_idx = 0;
+      if (segments_[j].horizontal_.contains(x))
+        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);
   
-    
-  return best_idx;
+  c.x_ = x;
+  c.y_ = y;
+  c.base_penalty_ = score_factor;
+  collisions_.push_back (c);
 }
-  
-MAKE_SCHEME_CALLBACK (Beam, quanting, 1);
-SCM
-Beam::quanting (SCM smob)
+
+void Beam_scoring_problem::init_collisions (vector<Grob*> grobs)
 {
-  Grob *me = unsmob_grob (smob);
+  Grob* common_x = NULL;
+  segments_ = Beam::get_beam_segments (beam, &common_x);
+  vector_sort (segments_, beam_segment_less);
+  if (common[X_AXIS] != common_x)
+    {
+      programming_error ("Disagree on common x. Skipping collisions in beam scoring.");
+      return;
+    }
 
-  SCM s = me->get_property ("positions");
-  Real yl = scm_to_double (scm_car (s));
-  Real yr = scm_to_double (scm_cdr (s));
+  set<Grob*> stems;
+  for (vsize i = 0; i < grobs.size (); i++) {
+    Box b;
 
+    for (Axis a = X_AXIS; a < NO_AXES; incr (a))
+      b[a] = grobs[i]->extent(common[a], a);
 
-  /*
-    Calculations are relative to a unit-scaled staff, i.e. the quants are
-    divided by the current staff_space.
-    
-   */
-  Real ss = Staff_symbol_referencer::staff_space (me);
-  Real thickness = Beam::get_thickness (me) / ss ;
-  Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
+    Real width = b[X_AXIS].length();
 
-  Real dy_mus = robust_scm2double (me->get_property ("least-squares-dy"), 0);
-  Real straddle = 0.0;
-  Real sit = (thickness - slt) / 2;
-  Real inter = 0.5;
-  Real hang = 1.0 - (thickness - slt) / 2;
-  Real quants [] = {straddle, sit, inter, hang };
+    Real width_factor = sqrt(width / staff_space);
 
+    Direction d = LEFT;
+    do
+      add_collision (b[X_AXIS][d], b[Y_AXIS], width_factor);
+    while (flip (&d) != LEFT);
 
+    Grob* stem = unsmob_grob (grobs[i]->get_object ("stem"));
+    if (stem && Stem::has_interface (stem) && Stem::is_normal_stem (stem))
+      {
+        stems.insert (stem);
+      }
+  }
   
-  int num_quants = int (sizeof (quants)/sizeof (Real));
-  Array<Real> quantsl;
-  Array<Real> quantsr;
-
-  /*
-    going to REGION_SIZE == 2, yields another 0.6 second with
-    wtk1-fugue2.
-
-
-    (result indexes between 70 and 575)  ? --hwn. 
-
-  */
-
-
+  for (set<Grob*>::const_iterator it(stems.begin ()); it != stems.end (); it++)
+    {
+      Grob *s = *it;
+      Real x = s->extent (common[X_AXIS], X_AXIS).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)
+        - beam->relative_coordinate (common[Y_AXIS], Y_AXIS);
+
+      Real factor = parameters.STEM_COLLISION_FACTOR;
+      if (!unsmob_grob (s->get_object ("beam"))
+          && !Stem::flag (s).is_empty ())
+        factor = 1.0; 
+      add_collision (x, y, factor);
+    }
+}
   
-  /*
-    Do stem computations.  These depend on YL and YR linearly, so we can
-    precompute for every stem 2 factors.
-   */
-  Link_array<Grob> stems =
-    Pointer_group_interface__extract_grobs (me, (Grob*)0, "stems");
-  Array<Stem_info> stem_infos;
-  Array<Real> base_lengths;
-  Array<Real> stem_xposns;  
-
-  Drul_array<bool> dirs_found (0,0);
-  Grob *common[2];
+void Beam_scoring_problem::init_stems ()
+{
+  extract_grob_set (beam, "covered-grobs", collisions);
+  extract_grob_set (beam, "stems", stems);
   for (int a = 2; a--;)
-    common[a] = common_refpoint_of_array (stems, me, Axis (a));
-
-  Grob * fvs = first_visible_stem (me);
-  Grob *lvs = last_visible_stem (me);
-  Real xl = fvs ? fvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
-  Real xr = fvs ? lvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
-
-  /*
-    We store some info to quickly interpolate.
-
-    Sometimes my head is screwed on backwards.  The stemlength are
-    AFFINE linear in YL and YR. If YL == YR == 0, then we might have
-    stem_y != 0.0, when we're cross staff.
-    
-   */
-  for (int i = 0; i < stems.size (); i++)
     {
-      Grob*s = stems[i];
-
+      common[a] = common_refpoint_of_array (stems, beam, Axis (a));
+      common[a] = common_refpoint_of_array (collisions, common[a], Axis (a));
+    }
+  
+  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);
+  
+  Drul_array<bool> dirs_found (0, 0);
+  for (vsize i = 0; i < stems.size (); i++)
+    {
+      Grob *s = stems[i];
+      if (!Stem::is_normal_stem (s))
+        continue;
+      
       Stem_info si (Stem::get_stem_info (s));
-      si.scale (1 / ss);
-      stem_infos.push (si);
-      dirs_found[stem_infos.top ().dir_] = true;
+      si.scale (1 / staff_space);
+      stem_infos.push_back (si);
+      dirs_found[si.dir_] = true;
 
       bool f = to_boolean (s->get_property ("french-beaming"))
-        && s != lvs && s != fvs;
+        && s != edge_stems[LEFT] && s != edge_stems[RIGHT];
 
-      base_lengths.push (calc_stem_y (me, s, common, xl, xr,
-                                     Interval (0,0), f) / ss);
-      stem_xposns.push (s->relative_coordinate (common[X_AXIS], X_AXIS));
+      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));
     }
-
-  bool xstaff = false;
-  if (lvs && fvs)
+  
+  edge_dirs = Drul_array<Direction> (CENTER, CENTER);
+  if (stem_infos.size ())
     {
-      Grob *commony = fvs->common_refpoint (lvs, Y_AXIS);
-      xstaff = Align_interface::has_interface (commony);
+      edge_dirs = Drul_array<Direction> (stem_infos[0].dir_,
+                                         stem_infos.back().dir_);
     }
+
+  is_xstaff = Align_interface::has_interface (common[Y_AXIS]);
+  is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
   
-  Direction ldir = Direction (stem_infos[0].dir_);
-  Direction rdir = Direction (stem_infos.top ().dir_);
-  bool is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
+  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);
+
+  // TODO - why are we dividing by staff_space?
+  beam_translation = Beam::get_beam_translation (beam) / staff_space;
+
+  d = LEFT;
+  do
+    {
+      quant_range[d].set_full ();
+      if (!edge_stems[d])
+        continue;
+      
+      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 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;
+    }
+  while (flip (&d) != LEFT);
 
+  init_collisions (collisions);
+}
 
-  int region_size = REGION_SIZE;
+Beam_scoring_problem::Beam_scoring_problem (Grob *me, Drul_array<Real> ys)
+{
+  beam = me;
+  unquanted_y = ys;
+    
   /*
-    Knees are harder, lets try some more possibilities for knees. 
-   */
+    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 (me);
+  beam_thickness = Beam::get_beam_thickness (me) / staff_space;
+  line_thickness = Staff_symbol_referencer::line_thickness (me) / staff_space;
+
+  // This is the least-squares DY, corrected for concave beams.
+  musical_dy = robust_scm2double (me->get_property ("least-squares-dy"), 0);
+
+  parameters.fill (me);
+  init_stems ();
+}
+
+void
+Beam_scoring_problem::generate_quants (vector<Beam_configuration*> *scores) const
+{
+  int region_size = (int) parameters.REGION_SIZE;
+
+  // Knees and collisions are harder, lets try some more possibilities
   if (is_knee)
     region_size += 2;
+  if (collisions_.size ())
+    region_size += 2;
+  
+  Real straddle = 0.0;
+  Real sit = (beam_thickness - line_thickness) / 2;
+  Real inter = 0.5;
+  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));
 
   /*
     Asymetry ? should run to <= region_size ?
-   */
-  for (int i = -region_size ; i < region_size; i++)
-    for (int j = 0; j < num_quants; j++)
+  */
+  vector<Real> unshifted_quants;
+  for (int i = -region_size; i < region_size; i++)
+    for (int j = 0; j < num_base_quants; j++)
       {
-       quantsl.push (i + quants[j] + int (yl));
-       quantsr.push (i + quants[j] + int (yr));
+        unshifted_quants.push_back (i + base_quants[j]);
       }
 
-  Array<Quant_score> qscores;
-  
-  for (int l = 0; l < quantsl.size (); l++)  
-    for (int r = 0; r < quantsr.size (); r++)
+  for (vsize i = 0; i < unshifted_quants.size (); i++)
+    for (vsize j = 0; j < unshifted_quants.size (); j++)
       {
-       Quant_score qs;
-       qs.yl = quantsl[l];
-       qs.yr = quantsr[r];
-       qs.demerits = 0.0;
-       
-       qscores.push (qs);
+        Beam_configuration *c =
+          Beam_configuration::new_config (unquanted_y,
+                                          Interval (unshifted_quants[i],
+                                                    unshifted_quants[j]));
+        
+        Direction d = LEFT;
+        do
+          {
+            if (!quant_range[d].contains (c->y[d]))
+              {
+                delete c;
+                c = NULL;
+                break;
+              }
+          }
+        while (flip (&d) != LEFT);
+        if (c)   
+          scores->push_back (c);
       }
+    
+}
 
-  /* This is a longish function, but we don't separate this out into
-     neat modular separate subfunctions, as the subfunctions would be
-     called for many values of YL, YR. By precomputing various
-     parameters outside of the loop, we can save a lot of time. */
-  for (int i = qscores.size (); i--;)
-    {
-      Real d =  score_slopes_dy (qscores[i].yl, qscores[i].yr,
-                                dy_mus, yr- yl, 
-                                xr - xl,
-                                xstaff);
-      qscores[i].demerits += d;
-
-#if DEBUG_QUANTING
-      qscores[i].score_card_ += to_string ("S%.2f",d);
-#endif
-    }
 
-  Real rad = Staff_symbol_referencer::staff_radius (me);
-  Drul_array<int> edge_beam_counts
-    (Stem::beam_multiplicity (stems[0]).length  () + 1,
-     Stem::beam_multiplicity (stems.top ()).length  () + 1);
-  
-  Real beam_translation = get_beam_translation (me) / ss;
+void Beam_scoring_problem::one_scorer (Beam_configuration* config) const
+{
+  score_count ++;
+  switch (config->next_scorer_todo) {
+  case SLOPE_IDEAL:
+    score_slope_ideal (config);
+    break;
+  case SLOPE_DIRECTION:
+    score_slope_direction (config);
+    break;
+  case SLOPE_MUSICAL:
+    score_slope_musical (config);
+    break;
+  case FORBIDDEN:
+    score_forbidden_quants (config);
+    break;
+  case STEM_LENGTHS:
+    score_stem_lengths (config);
+    break;
+  case COLLISIONS:
+    score_collisions (config);
+    break;
+  case HORIZONTAL_INTER:
+    score_horizontal_inter_quants (config);
+    break;
+    
+  case NUM_SCORERS:
+  case ORIGINAL_DISTANCE:
+  default:
+    assert (false);
+  }
+  config->next_scorer_todo++;
+}                                  
 
-  Real reasonable_score = (is_knee) ? 200000 : 100;
-  for (int i = qscores.size (); i--;)
-    if (qscores[i].demerits < reasonable_score)
-      {
-       Real d = score_forbidden_quants (qscores[i].yl, qscores[i].yr,
-                                        rad, slt, thickness, beam_translation,
-                                        edge_beam_counts, ldir, rdir); 
-       qscores[i].demerits += d;
 
-#if DEBUG_QUANTING
-       qscores[i].score_card_ += to_string (" F %.2f", d);
-#endif
-      }
+Beam_configuration *
+Beam_scoring_problem::force_score (SCM inspect_quants, const vector<Beam_configuration*> &configs) const
+{
+  Drul_array<Real> ins = ly_scm2interval (inspect_quants);
+  Real mindist = 1e6;
+  Beam_configuration *best = NULL; 
+  for (vsize i = 0; i < configs.size (); i++)
+    {
+      Real d = fabs (configs[i]->y[LEFT]- ins[LEFT]) + fabs (configs[i]->y[RIGHT] - ins[RIGHT]);
+      if (d < mindist)
+        {
+          best = configs[i];
+          mindist = d;
+        }
+    }
+  if (mindist > 1e5)
+    programming_error ("cannot find quant");
 
-  for (int i = qscores.size (); i--;)
-    if (qscores[i].demerits < reasonable_score)
-      {
-       Real d = score_stem_lengths (stems, stem_infos,
-                                base_lengths, stem_xposns,
-                                xl, xr,
-                                is_knee,
-                                qscores[i].yl, qscores[i].yr);
-       qscores[i].demerits +=  d;
-
-#if DEBUG_QUANTING
-       qscores[i].score_card_ += to_string (" L %.2f", d);
-#endif
-      }
+  while (!best->done ())
+    one_scorer (best);
+      
+  return best;
+}
+
+Drul_array<Real>
+Beam_scoring_problem::solve () const {
+  vector<Beam_configuration*> configs;
+  generate_quants (&configs);
 
-  int best_idx = best_quant_score_idx (qscores);
+  Beam_configuration *best = NULL;  
 
-#if DEBUG_QUANTING
-  SCM inspect_quants = me->get_property ("inspect-quants");
-  if (to_boolean (me->get_layout ()->lookup_variable (ly_symbol2scm ("debug-beam-quanting")))
-      && scm_is_pair (inspect_quants))
+  bool debug =
+    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;
+      best = force_score (inspect_quants, configs);
+    }
+  else
     {
-      Drul_array<Real> ins = ly_scm2interval (inspect_quants);
+      std::priority_queue<Beam_configuration*, std::vector<Beam_configuration*>,
+                          Beam_configuration_less> queue;
+      for (vsize i = 0; i < configs.size(); i++)
+        queue.push(configs[i]);
+
+      /*
+        TODO
+
+        It would be neat if we generated new configurations on the
+        fly, depending on the best complete score so far, eg.
+
+        if (best->done()) {
+          if (best->demerits < sqrt(queue.size())
+            break;
+          while (best->demerits > sqrt(queue.size()) {
+            generate and insert new configuration
+          }
+        }
+
+        that would allow us to do away with region_size altogether.
+      */
+      while (true) {
+        best = queue.top ();
+        if (best->done ())
+          break;
+
+        queue.pop ();
+        one_scorer (best);
+        queue.push (best);
+      }
+    }
 
-      int i = 0;
+  Interval final_positions = best->y;
 
-      Real mindist = 1e6;
-      for (; i < qscores.size (); i ++)
-       {
-         Real d = fabs (qscores[i].yl- ins[LEFT]) + fabs (qscores[i].yr - ins[RIGHT]);
-         if (d < mindist)
-           {
-             best_idx = i;
-             mindist = d;
-           }
-       }
-      if (mindist > 1e5)
-       programming_error ("Could not find quant.");
-    }
-#endif
-  
-  me->set_property ("positions",
-                   ly_interval2scm (Drul_array<Real> (qscores[best_idx].yl,
-                                                      qscores[best_idx].yr)));
-#if DEBUG_QUANTING
-  if (to_boolean (me->get_layout ()->lookup_variable (ly_symbol2scm ("debug-beam-quanting"))))
+#if DEBUG_BEAM_SCORING
+  if (debug)
     {
-      qscores[best_idx].score_card_ += to_string ("i%d", best_idx);
-      
       // debug quanting
-      me->set_property ("quant-score",
-                       scm_makfrom0str (qscores[best_idx].score_card_.to_str0 ()));
+      int completed = 0;
+      for (vsize i = 0; i < configs.size (); i++)
+        {
+          if (configs[i]->done ())
+            completed++;
+        }
+
+      string card = best->score_card_ + to_string (" c%d/%d", completed, configs.size());
+      beam->set_property ("annotation", ly_string2scm (card));
     }
 #endif
 
-  return SCM_UNSPECIFIED;
+  junk_pointers (configs);
+  return final_positions;
 }
 
-Real
-Beam::score_stem_lengths (Link_array<Grob> const &stems,
-                         Array<Stem_info> const &stem_infos,
-                         Array<Real> const &base_stem_ys,
-                         Array<Real> const &stem_xs,
-                         Real xl, Real xr, 
-                         bool knee, 
-                         Real yl, Real yr)
+void
+Beam_scoring_problem::score_stem_lengths (Beam_configuration* config) const
 {
-  Real limit_penalty = 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 (int i = 0; i < stems.size (); i++)
+
+  for (vsize i = 0; i < stem_xpositions.size (); i++)
     {
-      Grob* s = stems[i];
-      if (Stem::is_invisible (s))
-       continue;
-
-      Real x = stem_xs[i];
-      Real dx = xr-xl;
-      Real beam_y = dx ? yr *(x - xl)/dx + yl * ( xr - x)/dx : (yr + yl)/2;
-      Real current_y = beam_y + base_stem_ys[i];
-      Real length_pen = STEM_LENGTH_DEMERIT_FACTOR;
-      
+      Real x = stem_xpositions[i];
+      Real dx = x_span.delta ();
+      Real beam_y = dx
+        ? config->y[RIGHT] * (x - x_span[LEFT]) / dx + config->y[LEFT] * (x_span[RIGHT] - 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;
+
       Stem_info info = stem_infos[i];
       Direction d = info.dir_;
 
-      score[d] += limit_penalty * (0 >? (d * (info.shortest_y_ - current_y)));
-      
+      score[d] += limit_penalty * max (0.0, (d * (info.shortest_y_ - current_y)));
+
       Real ideal_diff = d * (current_y - info.ideal_y_);
       Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
-      
+
       /* 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 (knee)
+      if (is_knee)
        ideal_score = pow (ideal_score, 1.1);
-      
-      score[d] += length_pen * ideal_score;
 
-      count[d] ++;
+      score[d] += length_pen * ideal_score;
+      count[d]++;
     }
 
+  /* Divide by number of stems, to make the measure scale-free. */
   Direction d = DOWN;
   do
-    { 
-      score[d] /= (count[d] >? 1);
-    }
+    score[d] /= max (count[d], 1);
   while (flip (&d) != DOWN);
 
-  return score[LEFT]+score[RIGHT];
+  config->add (score[LEFT] + score[RIGHT], "L");
 }
 
-Real
-Beam::score_slopes_dy (Real yl, Real yr,
-                      Real dy_mus, Real dy_damp,
-                      Real dx,
-                      bool xstaff)
+void
+Beam_scoring_problem::score_slope_direction (Beam_configuration *config) const
 {
-  Real dy = yr - yl;
+  Real dy = config->y.delta ();
+  Real damped_dy = unquanted_y.delta();
   Real dem = 0.0;
-
   /*
     DAMPING_DIRECTION_PENALTY is a very harsh measure, while for
     complex beaming patterns, horizontal is often a good choice.
 
     TODO: find a way to incorporate the complexity of the beam in this
     penalty.
-   */
-  if (fabs (dy/dx) > ROUND_TO_ZERO_SLOPE
-      && sign (dy_damp) != sign (dy))
+  */
+  if (sign (damped_dy) != sign (dy))
     {
-      dem += DAMPING_DIRECTION_PENALTY;
+      if (!dy)
+       {
+         if (fabs (damped_dy / x_span.delta ()) > parameters.ROUND_TO_ZERO_SLOPE)
+           dem += parameters.DAMPING_DIRECTION_PENALTY;
+         else
+           dem += parameters.HINT_DIRECTION_PENALTY;
+       }
+      else
+       dem += parameters.DAMPING_DIRECTION_PENALTY;
     }
 
-   dem += MUSICAL_DIRECTION_FACTOR * (0 >? (fabs (dy) - fabs (dy_mus)));
+  config->add (dem, "Sd");
+}
 
+// Score for going against the direction of the musical pattern 
+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)));
+  config->add (dem, "Sm");
+}
 
-   Real slope_penalty = IDEAL_SLOPE_FACTOR;
+// Score deviation from calculated ideal slope.
+void
+Beam_scoring_problem::score_slope_ideal (Beam_configuration *config) const
+{
+  Real dy = config->y.delta ();
+  Real damped_dy = unquanted_y.delta();
+  Real dem = 0.0;
+  
+  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 (xstaff)
-     slope_penalty *= 10;
+  /* Xstaff beams tend to use extreme slopes to get short stems. We
+     put in a penalty here. */
+  if (is_xstaff)
+    slope_penalty *= 10;
 
-   /* Huh, why would a too steep beam be better than a too flat one ? */
-   dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
-     * slope_penalty;
+  /* Huh, why would a too steep beam be better than a too flat one ? */
+  dem += shrink_extra_weight (fabs (damped_dy) - fabs (dy), 1.5)
+    * slope_penalty;
 
-   return dem;
+  config->add (dem, "Si");
 }
 
-
 static Real
 my_modf (Real x)
 {
   return x - floor (x);
 }
 
+// TODO - there is some overlap with forbidden quants, but for
+// horizontal beams, it is much more serious to have stafflines
+// appearing in the wrong place, so we have a separate scorer.
+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)
+    {
+      Real yshift = config->y[LEFT] - 0.5 * staff_space;
+      if (abs (round(yshift) - yshift) < 0.01 * staff_space)
+        config->add (parameters.HORIZONTAL_INTER_QUANT_PENALTY, "H");
+    }
+}
 
 /*
   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
-   because for 32nd and 64th beams the forbidden quants are relatively
-   more important than stem lengths.
+  because for 32nd and 64th beams the forbidden quants are relatively
+  more important than stem lengths.
 */
-Real
-Beam::score_forbidden_quants (Real yl, Real yr,
-                             Real radius,
-                             Real slt,
-                             Real thickness, Real beam_translation,
-                             Drul_array<int> beam_counts,
-                             Direction ldir, Direction rdir)
+void
+Beam_scoring_problem::score_forbidden_quants (Beam_configuration *config) const
 {
-  Real dy = yr - yl;
-  Drul_array<Real> y (yl,yr);
-  Drul_array<Direction> dirs (ldir,rdir);
-  
-  Real extra_demerit = SECONDARY_BEAM_DEMERIT / (beam_counts[LEFT] >? beam_counts[RIGHT]);
+  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 dem = 0.0;
+  Real eps = parameters.BEAM_EPS;
   
-
   do
     {
-      for (int j = 1; j <= beam_counts[d]; j++)
+      for (int j = 1; j <= edge_beam_counts[d]; j++)
        {
-         Direction stem_dir = dirs[d];
+         Direction stem_dir = edge_dirs[d];
 
          /*
            The 2.2 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
+           will be in the gap of the (2, sit) quant, leading to a
            false demerit.
          */
-         Real gap1 = y[d] - stem_dir * ((j-1) * beam_translation + thickness / 2 - slt/2.2 );
-         Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt/2.2);
+         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);
 
          Interval gap;
          gap.add_point (gap1);
          gap.add_point (gap2);
 
-         for (Real k = - radius ;
-              k <= radius + BEAM_EPS; k += 1.0) 
+         for (Real k = -staff_radius;
+              k <= staff_radius + eps; k += 1.0)
            if (gap.contains (k))
              {
-               Real dist = fabs (gap[UP]-k) <? fabs (gap[DOWN] - k);
+               Real dist = min (fabs (gap[UP] - k), fabs (gap[DOWN] - k));
 
                /*
                  this parameter is tuned to grace-stem-length.ly
                */
                Real fixed_demerit = 0.4;
-               
+
                dem += extra_demerit
-                 * (fixed_demerit +
-                    (1-fixed_demerit) * (dist / gap.length())* 2);
+                 * (fixed_demerit
+                    + (1 - fixed_demerit) * (dist / gap.length ()) * 2);
              }
        }
     }
-  while ((flip (&d))!= LEFT); 
-
+  while ((flip (&d)) != LEFT);
 
-  if ((beam_counts[LEFT] >? beam_counts[RIGHT]) >= 2)
+  if (max (edge_beam_counts[LEFT], edge_beam_counts[RIGHT]) >= 2)
     {
       Real straddle = 0.0;
-      Real sit = (thickness - slt) / 2;
+      Real sit = (beam_thickness - line_thickness) / 2;
       Real inter = 0.5;
-      Real hang = 1.0 - (thickness - slt) / 2;
-
+      Real hang = 1.0 - (beam_thickness - line_thickness) / 2;
 
       Direction d = LEFT;
       do
        {
-         if (beam_counts[d] >= 2
-             && fabs (y[d] - dirs[d] * beam_translation) < radius + inter)
+         if (edge_beam_counts[d] >= 2
+             && fabs (config->y[d] - edge_dirs[d] * beam_translation) < staff_radius + inter)
            {
-             if (dirs[d] == UP && dy <= BEAM_EPS
-                 && fabs (my_modf (y[d]) - sit) < BEAM_EPS)
+              // TODO up/down symmetry.
+             if (edge_dirs[d] == UP && dy <= eps
+                 && fabs (my_modf (config->y[d]) - sit) < eps)
                dem += extra_demerit;
-         
-             if (dirs[d] == DOWN && dy >= BEAM_EPS
-                 && fabs (my_modf (y[d]) - hang) < BEAM_EPS)
+
+             if (edge_dirs[d] == DOWN && dy >= eps
+                 && fabs (my_modf (config->y[d]) - hang) < eps)
                dem += extra_demerit;
            }
 
-         if (beam_counts[d] >= 3
-             && fabs (y[d] - 2 * dirs[d] * beam_translation) < radius + inter)
+         if (edge_beam_counts[d] >= 3
+             && fabs (config->y[d] - 2 * edge_dirs[d] * beam_translation) < staff_radius + inter)
            {
-             if (dirs[d] == UP && dy <= BEAM_EPS
-                 && fabs (my_modf (y[d]) - straddle) < BEAM_EPS)
+              // TODO up/down symmetry.
+             if (edge_dirs[d] == UP && dy <= eps
+                 && fabs (my_modf (config->y[d]) - straddle) < eps)
                dem += extra_demerit;
-             
-             if (dirs[d] == DOWN && dy >= BEAM_EPS
-                 && fabs (my_modf (y[d]) - straddle) < BEAM_EPS)
+
+             if (edge_dirs[d] == DOWN && dy >= eps
+                 && fabs (my_modf (config->y[d]) - straddle) < eps)
                dem += extra_demerit;
            }
        }
       while (flip (&d) != LEFT);
     }
-  
-  return dem;
+
+  config->add (dem, "F");
 }
 
-  
+void
+Beam_scoring_problem::score_collisions (Beam_configuration *config) const
+{  
+  Real demerits = 0.0;
+  for (vsize i = 0; i < collisions_.size (); i++)
+    {
+      Interval collision_y = collisions_[i].y_;
+      Real x = collisions_[i].x_;
+
+      Real center_beam_y = y_at (x, config);
+      Interval beam_y = center_beam_y + collisions_[i].beam_y_;
+
+      Real dist = infinity_f;
+      if (!intersection (beam_y, collision_y).is_empty ())
+        dist = 0.0;
+      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;
+    }
+
+  config->add (demerits, "C");
+}