lily/beam-quanting.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 1997--2011 Han-Wen Nienhuys <hanwen@xs4all.nl>
   5   Jan Nieuwenhuizen <janneke@gnu.org>
   6
   7   LilyPond is free software: you can redistribute it and/or modify
   8   it under the terms of the GNU General Public License as published by
   9   the Free Software Foundation, either version 3 of the License, or
  10   (at your option) any later version.
  11
  12   LilyPond is distributed in the hope that it will be useful,
  13   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15   GNU General Public License for more details.
  16
  17   You should have received a copy of the GNU General Public License
  18   along with LilyPond.  If not, see <http://www.gnu.org/licenses/>.
  19 */
  20
  21 #include "beam-scoring-problem.hh"
  22
  23 #include <algorithm>
  24 #include <queue>
  25 #include <set>
  26 using namespace std;
  27
  28 #include "align-interface.hh"
  29 #include "beam.hh"
  30 #include "direction.hh"
  31 #include "directional-element-interface.hh"
  32 #include "grob.hh"
  33 #include "international.hh"
  34 #include "libc-extension.hh"
  35 #include "main.hh"
  36 #include "output-def.hh"
  37 #include "pointer-group-interface.hh"
  38 #include "staff-symbol-referencer.hh"
  39 #include "stencil.hh"
  40 #include "stem.hh"
  41 #include "warn.hh"
  42
  43 Real
  44 get_detail (SCM alist, SCM sym, Real def)
  45 {
  46   SCM entry = scm_assq (sym, alist);
  47
  48   if (scm_is_pair (entry))
  49     return robust_scm2double (scm_cdr (entry), def);
  50   return def;
  51 }
  52
  53 void
  54 Beam_quant_parameters::fill (Grob *him)
  55 {
  56   SCM details = him->get_property ("details");
  57
  58   // General
  59   BEAM_EPS = get_detail (details, ly_symbol2scm ("beam-eps"), 1e-3);
  60   REGION_SIZE = get_detail (details, ly_symbol2scm ("region-size"), 2);
  61
  62   // forbidden quants
  63   SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0);
  64   STEM_LENGTH_DEMERIT_FACTOR = get_detail (details, ly_symbol2scm ("stem-length-demerit-factor"), 5);
  65   HORIZONTAL_INTER_QUANT_PENALTY = get_detail (details, ly_symbol2scm ("horizontal-inter-quant"), 500);
  66
  67   STEM_LENGTH_LIMIT_PENALTY = get_detail (details, ly_symbol2scm ("stem-length-limit-penalty"), 5000);
  68   DAMPING_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("damping-direction-penalty"), 800);
  69   HINT_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("hint-direction-penalty"), 20);
  70   MUSICAL_DIRECTION_FACTOR = get_detail (details, ly_symbol2scm ("musical-direction-factor"), 400);
  71   IDEAL_SLOPE_FACTOR = get_detail (details, ly_symbol2scm ("ideal-slope-factor"), 10);
  72   ROUND_TO_ZERO_SLOPE = get_detail (details, ly_symbol2scm ("round-to-zero-slope"), 0.02);
  73
  74   // Collisions
  75   COLLISION_PENALTY = get_detail (details, ly_symbol2scm ("collision-penalty"), 500);
  76   COLLISION_PADDING = get_detail (details, ly_symbol2scm ("collision-padding"), 0.5);
  77   STEM_COLLISION_FACTOR = get_detail (details, ly_symbol2scm ("stem-collision-factor"), 0.1);
  78 }
  79
  80 // Add x if x is positive, add |x|*fac if x is negative.
  81 static Real
  82 shrink_extra_weight (Real x, Real fac)
  83 {
  84   return fabs (x) * ((x < 0) ? fac : 1.0);
  85 }
  86
  87 /****************************************************************/
  88
  89 Beam_configuration::Beam_configuration ()
  90 {
  91   y = Interval (0.0, 0.0);
  92   demerits = 0.0;
  93   next_scorer_todo = ORIGINAL_DISTANCE;
  94 }
  95
  96 bool Beam_configuration::done () const
  97 {
  98   return next_scorer_todo >= NUM_SCORERS;
  99 }
 100
 101 void Beam_configuration::add (Real demerit, const string &reason)
 102 {
 103   demerits += demerit;
 104
 105 #if DEBUG_BEAM_SCORING
 106   if (demerit)
 107     score_card_ += to_string (" %s %.2f", reason.c_str (), demerit);
 108 #endif
 109 }
 110
 111 Beam_configuration *Beam_configuration::new_config (Interval start,
 112                                                     Interval offset)
 113 {
 114   Beam_configuration *qs = new Beam_configuration;
 115   qs->y = Interval (int (start[LEFT]) + offset[LEFT],
 116                     int (start[RIGHT]) + offset[RIGHT]);
 117
 118   // This orders the sequence so we try combinations closest to the
 119   // the ideal offset first.
 120   Real start_score = abs (offset[RIGHT]) + abs (offset[LEFT]);
 121   qs->demerits = start_score / 1000.0;
 122   qs->next_scorer_todo = ORIGINAL_DISTANCE + 1;
 123
 124   return qs;
 125 }
 126
 127 Real
 128 Beam_scoring_problem::y_at (Real x, Beam_configuration const *p) const
 129 {
 130   return p->y[LEFT] + (x - x_span[LEFT]) * p->y.delta () / x_span.delta ();
 131 }
 132
 133 /****************************************************************/
 134
 135 /*
 136   TODO:
 137
 138   - Make all demerits customisable
 139
 140   - Add demerits for quants per se, as to forbid a specific quant
 141   entirely
 142 */
 143
 144 // This is a temporary hack to see how much we can gain by using a
 145 // priority queue on the beams to score.
 146 static int score_count = 0;
 147 LY_DEFINE (ly_beam_score_count, "ly:beam-score-count", 0, 0, 0,
 148            (),
 149            "count number of beam scores.")
 150 {
 151   return scm_from_int (score_count);
 152 }
 153
 154 void Beam_scoring_problem::add_collision (Real x, Interval y,
 155                                           Real score_factor)
 156 {
 157   if (edge_dirs[LEFT] == edge_dirs[RIGHT])
 158     {
 159       Direction d = edge_dirs[LEFT];
 160
 161       Real quant_range_y = quant_range[LEFT][-d]
 162                            + (x - x_span[LEFT]) * (quant_range[RIGHT][-d] - quant_range[LEFT][-d]) / x_span.delta ();
 163
 164       if (d * (quant_range_y - minmax (d, y[UP], y[DOWN])) > 0)
 165         {
 166           return;
 167         }
 168     }
 169
 170   Beam_collision c;
 171   c.beam_y_.set_empty ();
 172
 173   for (vsize j = 0; j < segments_.size (); j++)
 174     {
 175       if (segments_[j].horizontal_.contains (x))
 176         c.beam_y_.add_point (segments_[j].vertical_count_ * beam_translation);
 177       if (segments_[j].horizontal_[LEFT] > x)
 178         break;
 179     }
 180   c.beam_y_.widen (0.5 * beam_thickness);
 181
 182   c.x_ = x;
 183
 184   y *= 1 / staff_space;
 185   c.y_ = y;
 186   c.base_penalty_ = score_factor;
 187   collisions_.push_back (c);
 188 }
 189
 190 void Beam_scoring_problem::init_collisions (vector<Grob *> grobs)
 191 {
 192   Grob *common_x = NULL;
 193   segments_ = Beam::get_beam_segments (beam, &common_x);
 194   vector_sort (segments_, beam_segment_less);
 195   if (common[X_AXIS] != common_x)
 196     {
 197       programming_error ("Disagree on common x. Skipping collisions in beam scoring.");
 198       return;
 199     }
 200
 201   set<Grob *> stems;
 202   for (vsize i = 0; i < grobs.size (); i++)
 203     {
 204       Box b;
 205       for (Axis a = X_AXIS; a < NO_AXES; incr (a))
 206         b[a] = grobs[i]->extent (common[a], a);
 207
 208       Real width = b[X_AXIS].length ();
 209       Real width_factor = sqrt (width / staff_space);
 210
 211       Direction d = LEFT;
 212       do
 213         add_collision (b[X_AXIS][d], b[Y_AXIS], width_factor);
 214       while (flip (&d) != LEFT);
 215
 216       Grob *stem = unsmob_grob (grobs[i]->get_object ("stem"));
 217       if (stem && Stem::has_interface (stem) && Stem::is_normal_stem (stem))
 218         {
 219           stems.insert (stem);
 220         }
 221     }
 222
 223   for (set<Grob *>::const_iterator it (stems.begin ()); it != stems.end (); it++)
 224     {
 225       Grob *s = *it;
 226       Real x = s->extent (common[X_AXIS], X_AXIS).center ();
 227
 228       Direction stem_dir = get_grob_direction (*it);
 229       Interval y;
 230       y.set_full ();
 231       y[-stem_dir] = Stem::chord_start_y (*it) + (*it)->relative_coordinate (common[Y_AXIS], Y_AXIS)
 232                      - beam->relative_coordinate (common[Y_AXIS], Y_AXIS);
 233
 234       Real factor = parameters.STEM_COLLISION_FACTOR;
 235       if (!unsmob_grob (s->get_object ("beam"))
 236           && !Stem::flag (s).is_empty ())
 237         factor = 1.0;
 238       add_collision (x, y, factor);
 239     }
 240 }
 241
 242 void Beam_scoring_problem::init_stems ()
 243 {
 244   extract_grob_set (beam, "covered-grobs", collisions);
 245   extract_grob_set (beam, "stems", stems);
 246   for (int a = 2; a--;)
 247     {
 248       common[a] = common_refpoint_of_array (stems, beam, Axis (a));
 249       common[a] = common_refpoint_of_array (collisions, common[a], Axis (a));
 250     }
 251
 252   Drul_array<Grob *> edge_stems (Beam::first_normal_stem (beam),
 253                                  Beam::last_normal_stem (beam));
 254   Direction d = LEFT;
 255   do
 256     x_span[d] = edge_stems[d] ? edge_stems[d]->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 257   while (flip (&d) != LEFT);
 258
 259   Drul_array<bool> dirs_found (0, 0);
 260   for (vsize i = 0; i < stems.size (); i++)
 261     {
 262       Grob *s = stems[i];
 263       if (!Stem::is_normal_stem (s))
 264         continue;
 265
 266       Stem_info si (Stem::get_stem_info (s));
 267       si.scale (1 / staff_space);
 268       stem_infos.push_back (si);
 269       dirs_found[si.dir_] = true;
 270
 271       bool f = to_boolean (s->get_property ("french-beaming"))
 272                && s != edge_stems[LEFT] && s != edge_stems[RIGHT];
 273
 274       Real y = Beam::calc_stem_y (beam, s, common, x_span[LEFT], x_span[RIGHT], CENTER,
 275                                   Interval (0, 0), f);
 276       base_lengths.push_back (y / staff_space);
 277       stem_xpositions.push_back (s->relative_coordinate (common[X_AXIS], X_AXIS));
 278     }
 279
 280   edge_dirs = Drul_array<Direction> (CENTER, CENTER);
 281   if (stem_infos.size ())
 282     {
 283       edge_dirs = Drul_array<Direction> (stem_infos[0].dir_,
 284                                          stem_infos.back ().dir_);
 285     }
 286
 287   is_xstaff = Align_interface::has_interface (common[Y_AXIS]);
 288   is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
 289
 290   staff_radius = Staff_symbol_referencer::staff_radius (beam);
 291   edge_beam_counts = Drul_array<int>
 292                      (Stem::beam_multiplicity (stems[0]).length () + 1,
 293                       Stem::beam_multiplicity (stems.back ()).length () + 1);
 294
 295   // TODO - why are we dividing by staff_space?
 296   beam_translation = Beam::get_beam_translation (beam) / staff_space;
 297
 298   d = LEFT;
 299   do
 300     {
 301       quant_range[d].set_full ();
 302       if (!edge_stems[d])
 303         continue;
 304
 305       Real stem_offset = edge_stems[d]->relative_coordinate (common[Y_AXIS], Y_AXIS)
 306                          - beam->relative_coordinate (common[Y_AXIS], Y_AXIS);
 307       Interval heads = Stem::head_positions (edge_stems[d]) * 0.5 * staff_space;
 308
 309       Direction ed = edge_dirs[d];
 310       heads.widen (0.5 * staff_space
 311                    + (edge_beam_counts[d] - 1) * beam_translation + beam_thickness * .5);
 312       quant_range[d][-ed] = heads[ed] + stem_offset;
 313     }
 314   while (flip (&d) != LEFT);
 315
 316   init_collisions (collisions);
 317 }
 318
 319 Beam_scoring_problem::Beam_scoring_problem (Grob *me, Drul_array<Real> ys)
 320 {
 321   beam = me;
 322   unquanted_y = ys;
 323
 324   /*
 325     Calculations are relative to a unit-scaled staff, i.e. the quants are
 326     divided by the current staff_space.
 327   */
 328   staff_space = Staff_symbol_referencer::staff_space (me);
 329   beam_thickness = Beam::get_beam_thickness (me) / staff_space;
 330   line_thickness = Staff_symbol_referencer::line_thickness (me) / staff_space;
 331
 332   // This is the least-squares DY, corrected for concave beams.
 333   musical_dy = robust_scm2double (me->get_property ("least-squares-dy"), 0);
 334
 335   parameters.fill (me);
 336   init_stems ();
 337 }
 338
 339 void
 340 Beam_scoring_problem::generate_quants (vector<Beam_configuration *> *scores) const
 341 {
 342   int region_size = (int) parameters.REGION_SIZE;
 343
 344   // Knees and collisions are harder, lets try some more possibilities
 345   if (is_knee)
 346     region_size += 2;
 347   if (collisions_.size ())
 348     region_size += 2;
 349
 350   Real straddle = 0.0;
 351   Real sit = (beam_thickness - line_thickness) / 2;
 352   Real inter = 0.5;
 353   Real hang = 1.0 - (beam_thickness - line_thickness) / 2;
 354   Real base_quants [] = {straddle, sit, inter, hang};
 355   int num_base_quants = int (sizeof (base_quants) / sizeof (Real));
 356
 357   /*
 358     Asymetry ? should run to <= region_size ?
 359   */
 360   vector<Real> unshifted_quants;
 361   for (int i = -region_size; i < region_size; i++)
 362     for (int j = 0; j < num_base_quants; j++)
 363       {
 364         unshifted_quants.push_back (i + base_quants[j]);
 365       }
 366
 367   for (vsize i = 0; i < unshifted_quants.size (); i++)
 368     for (vsize j = 0; j < unshifted_quants.size (); j++)
 369       {
 370         Beam_configuration *c
 371           = Beam_configuration::new_config (unquanted_y,
 372                                             Interval (unshifted_quants[i],
 373                                                       unshifted_quants[j]));
 374
 375         Direction d = LEFT;
 376         do
 377           {
 378             if (!quant_range[d].contains (c->y[d]))
 379               {
 380                 delete c;
 381                 c = NULL;
 382                 break;
 383               }
 384           }
 385         while (flip (&d) != LEFT);
 386         if (c)
 387           scores->push_back (c);
 388       }
 389
 390 }
 391
 392 void Beam_scoring_problem::one_scorer (Beam_configuration *config) const
 393 {
 394   score_count++;
 395   switch (config->next_scorer_todo)
 396     {
 397     case SLOPE_IDEAL:
 398       score_slope_ideal (config);
 399       break;
 400     case SLOPE_DIRECTION:
 401       score_slope_direction (config);
 402       break;
 403     case SLOPE_MUSICAL:
 404       score_slope_musical (config);
 405       break;
 406     case FORBIDDEN:
 407       score_forbidden_quants (config);
 408       break;
 409     case STEM_LENGTHS:
 410       score_stem_lengths (config);
 411       break;
 412     case COLLISIONS:
 413       score_collisions (config);
 414       break;
 415     case HORIZONTAL_INTER:
 416       score_horizontal_inter_quants (config);
 417       break;
 418
 419     case NUM_SCORERS:
 420     case ORIGINAL_DISTANCE:
 421     default:
 422       assert (false);
 423     }
 424   config->next_scorer_todo++;
 425 }
 426
 427 Beam_configuration *
 428 Beam_scoring_problem::force_score (SCM inspect_quants, const vector<Beam_configuration *> &configs) const
 429 {
 430   Drul_array<Real> ins = ly_scm2interval (inspect_quants);
 431   Real mindist = 1e6;
 432   Beam_configuration *best = NULL;
 433   for (vsize i = 0; i < configs.size (); i++)
 434     {
 435       Real d = fabs (configs[i]->y[LEFT] - ins[LEFT]) + fabs (configs[i]->y[RIGHT] - ins[RIGHT]);
 436       if (d < mindist)
 437         {
 438           best = configs[i];
 439           mindist = d;
 440         }
 441     }
 442   if (mindist > 1e5)
 443     programming_error ("cannot find quant");
 444
 445   while (!best->done ())
 446     one_scorer (best);
 447
 448   return best;
 449 }
 450
 451 Drul_array<Real>
 452 Beam_scoring_problem::solve () const
 453 {
 454   vector<Beam_configuration *> configs;
 455   generate_quants (&configs);
 456
 457   if (configs.empty ())
 458     {
 459       programming_error ("No viable beam quanting found.  Using unquanted y value.");
 460       return unquanted_y;
 461     }
 462
 463   Beam_configuration *best = NULL;
 464
 465   bool debug
 466     = to_boolean (beam->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring")));
 467   SCM inspect_quants = beam->get_property ("inspect-quants");
 468   if (scm_is_pair (inspect_quants))
 469     {
 470       debug = true;
 471       best = force_score (inspect_quants, configs);
 472     }
 473   else
 474     {
 475       std::priority_queue < Beam_configuration *, std::vector<Beam_configuration *>,
 476           Beam_configuration_less > queue;
 477       for (vsize i = 0; i < configs.size (); i++)
 478         queue.push (configs[i]);
 479
 480       /*
 481         TODO
 482
 483         It would be neat if we generated new configurations on the
 484         fly, depending on the best complete score so far, eg.
 485
 486         if (best->done()) {
 487           if (best->demerits < sqrt(queue.size())
 488             break;
 489           while (best->demerits > sqrt(queue.size()) {
 490             generate and insert new configuration
 491           }
 492         }
 493
 494         that would allow us to do away with region_size altogether.
 495       */
 496       while (true)
 497         {
 498           best = queue.top ();
 499           if (best->done ())
 500             break;
 501
 502           queue.pop ();
 503           one_scorer (best);
 504           queue.push (best);
 505         }
 506     }
 507
 508   Interval final_positions = best->y;
 509
 510 #if DEBUG_BEAM_SCORING
 511   if (debug)
 512     {
 513       // debug quanting
 514       int completed = 0;
 515       for (vsize i = 0; i < configs.size (); i++)
 516         {
 517           if (configs[i]->done ())
 518             completed++;
 519         }
 520
 521       string card = best->score_card_ + to_string (" c%d/%d", completed, configs.size ());
 522       beam->set_property ("annotation", ly_string2scm (card));
 523     }
 524 #endif
 525
 526   junk_pointers (configs);
 527   return final_positions;
 528 }
 529
 530 void
 531 Beam_scoring_problem::score_stem_lengths (Beam_configuration *config) const
 532 {
 533   Real limit_penalty = parameters.STEM_LENGTH_LIMIT_PENALTY;
 534   Drul_array<Real> score (0, 0);
 535   Drul_array<int> count (0, 0);
 536
 537   for (vsize i = 0; i < stem_xpositions.size (); i++)
 538     {
 539       Real x = stem_xpositions[i];
 540       Real dx = x_span.delta ();
 541       Real beam_y = dx
 542                     ? config->y[RIGHT] * (x - x_span[LEFT]) / dx + config->y[LEFT] * (x_span[RIGHT] - x) / dx
 543                     : (config->y[RIGHT] + config->y[LEFT]) / 2;
 544       Real current_y = beam_y + base_lengths[i];
 545       Real length_pen = parameters.STEM_LENGTH_DEMERIT_FACTOR;
 546
 547       Stem_info info = stem_infos[i];
 548       Direction d = info.dir_;
 549
 550       score[d] += limit_penalty * max (0.0, (d * (info.shortest_y_ - current_y)));
 551
 552       Real ideal_diff = d * (current_y - info.ideal_y_);
 553       Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
 554
 555       /* We introduce a power, to make the scoring strictly
 556          convex. Otherwise a symmetric knee beam (up/down/up/down)
 557          does not have an optimum in the middle. */
 558       if (is_knee)
 559         ideal_score = pow (ideal_score, 1.1);
 560
 561       score[d] += length_pen * ideal_score;
 562       count[d]++;
 563     }
 564
 565   /* Divide by number of stems, to make the measure scale-free. */
 566   Direction d = DOWN;
 567   do
 568     score[d] /= max (count[d], 1);
 569   while (flip (&d) != DOWN);
 570
 571   config->add (score[LEFT] + score[RIGHT], "L");
 572 }
 573
 574 void
 575 Beam_scoring_problem::score_slope_direction (Beam_configuration *config) const
 576 {
 577   Real dy = config->y.delta ();
 578   Real damped_dy = unquanted_y.delta ();
 579   Real dem = 0.0;
 580   /*
 581     DAMPING_DIRECTION_PENALTY is a very harsh measure, while for
 582     complex beaming patterns, horizontal is often a good choice.
 583
 584     TODO: find a way to incorporate the complexity of the beam in this
 585     penalty.
 586   */
 587   if (sign (damped_dy) != sign (dy))
 588     {
 589       if (!dy)
 590         {
 591           if (fabs (damped_dy / x_span.delta ()) > parameters.ROUND_TO_ZERO_SLOPE)
 592             dem += parameters.DAMPING_DIRECTION_PENALTY;
 593           else
 594             dem += parameters.HINT_DIRECTION_PENALTY;
 595         }
 596       else
 597         dem += parameters.DAMPING_DIRECTION_PENALTY;
 598     }
 599
 600   config->add (dem, "Sd");
 601 }
 602
 603 // Score for going against the direction of the musical pattern
 604 void
 605 Beam_scoring_problem::score_slope_musical (Beam_configuration *config) const
 606 {
 607   Real dy = config->y.delta ();
 608   Real dem = parameters.MUSICAL_DIRECTION_FACTOR
 609              * max (0.0, (fabs (dy) - fabs (musical_dy)));
 610   config->add (dem, "Sm");
 611 }
 612
 613 // Score deviation from calculated ideal slope.
 614 void
 615 Beam_scoring_problem::score_slope_ideal (Beam_configuration *config) const
 616 {
 617   Real dy = config->y.delta ();
 618   Real damped_dy = unquanted_y.delta ();
 619   Real dem = 0.0;
 620
 621   Real slope_penalty = parameters.IDEAL_SLOPE_FACTOR;
 622
 623   /* Xstaff beams tend to use extreme slopes to get short stems. We
 624      put in a penalty here. */
 625   if (is_xstaff)
 626     slope_penalty *= 10;
 627
 628   /* Huh, why would a too steep beam be better than a too flat one ? */
 629   dem += shrink_extra_weight (fabs (damped_dy) - fabs (dy), 1.5)
 630          * slope_penalty;
 631
 632   config->add (dem, "Si");
 633 }
 634
 635 static Real
 636 my_modf (Real x)
 637 {
 638   return x - floor (x);
 639 }
 640
 641 // TODO - there is some overlap with forbidden quants, but for
 642 // horizontal beams, it is much more serious to have stafflines
 643 // appearing in the wrong place, so we have a separate scorer.
 644 void
 645 Beam_scoring_problem::score_horizontal_inter_quants (Beam_configuration *config) const
 646 {
 647   if (config->y.delta () == 0.0
 648       && abs (config->y[LEFT]) < staff_radius * staff_space)
 649     {
 650       Real yshift = config->y[LEFT] - 0.5 * staff_space;
 651       if (fabs (my_round (yshift) - yshift) < 0.01 * staff_space)
 652         config->add (parameters.HORIZONTAL_INTER_QUANT_PENALTY, "H");
 653     }
 654 }
 655
 656 /*
 657   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
 658   because for 32nd and 64th beams the forbidden quants are relatively
 659   more important than stem lengths.
 660 */
 661 void
 662 Beam_scoring_problem::score_forbidden_quants (Beam_configuration *config) const
 663 {
 664   Real dy = config->y.delta ();
 665
 666   Real extra_demerit = parameters.SECONDARY_BEAM_DEMERIT
 667                        / max (edge_beam_counts[LEFT], edge_beam_counts[RIGHT]);
 668
 669   Direction d = LEFT;
 670   Real dem = 0.0;
 671   Real eps = parameters.BEAM_EPS;
 672
 673   do
 674     {
 675       for (int j = 1; j <= edge_beam_counts[d]; j++)
 676         {
 677           Direction stem_dir = edge_dirs[d];
 678
 679           /*
 680             The 2.2 factor is to provide a little leniency for
 681             borderline cases. If we do 2.0, then the upper outer line
 682             will be in the gap of the (2, sit) quant, leading to a
 683             false demerit.
 684           */
 685           Real gap1 = config->y[d] - stem_dir * ((j - 1) * beam_translation + beam_thickness / 2 - line_thickness / 2.2);
 686           Real gap2 = config->y[d] - stem_dir * (j * beam_translation - beam_thickness / 2 + line_thickness / 2.2);
 687
 688           Interval gap;
 689           gap.add_point (gap1);
 690           gap.add_point (gap2);
 691
 692           for (Real k = -staff_radius;
 693                k <= staff_radius + eps; k += 1.0)
 694             if (gap.contains (k))
 695               {
 696                 Real dist = min (fabs (gap[UP] - k), fabs (gap[DOWN] - k));
 697
 698                 /*
 699                   this parameter is tuned to grace-stem-length.ly
 700                 */
 701                 Real fixed_demerit = 0.4;
 702
 703                 dem += extra_demerit
 704                        * (fixed_demerit
 705                           + (1 - fixed_demerit) * (dist / gap.length ()) * 2);
 706               }
 707         }
 708     }
 709   while ((flip (&d)) != LEFT);
 710
 711   if (max (edge_beam_counts[LEFT], edge_beam_counts[RIGHT]) >= 2)
 712     {
 713       Real straddle = 0.0;
 714       Real sit = (beam_thickness - line_thickness) / 2;
 715       Real inter = 0.5;
 716       Real hang = 1.0 - (beam_thickness - line_thickness) / 2;
 717
 718       Direction d = LEFT;
 719       do
 720         {
 721           if (edge_beam_counts[d] >= 2
 722               && fabs (config->y[d] - edge_dirs[d] * beam_translation) < staff_radius + inter)
 723             {
 724               // TODO up/down symmetry.
 725               if (edge_dirs[d] == UP && dy <= eps
 726                   && fabs (my_modf (config->y[d]) - sit) < eps)
 727                 dem += extra_demerit;
 728
 729               if (edge_dirs[d] == DOWN && dy >= eps
 730                   && fabs (my_modf (config->y[d]) - hang) < eps)
 731                 dem += extra_demerit;
 732             }
 733
 734           if (edge_beam_counts[d] >= 3
 735               && fabs (config->y[d] - 2 * edge_dirs[d] * beam_translation) < staff_radius + inter)
 736             {
 737               // TODO up/down symmetry.
 738               if (edge_dirs[d] == UP && dy <= eps
 739                   && fabs (my_modf (config->y[d]) - straddle) < eps)
 740                 dem += extra_demerit;
 741
 742               if (edge_dirs[d] == DOWN && dy >= eps
 743                   && fabs (my_modf (config->y[d]) - straddle) < eps)
 744                 dem += extra_demerit;
 745             }
 746         }
 747       while (flip (&d) != LEFT);
 748     }
 749
 750   config->add (dem, "F");
 751 }
 752
 753 void
 754 Beam_scoring_problem::score_collisions (Beam_configuration *config) const
 755 {
 756   Real demerits = 0.0;
 757   for (vsize i = 0; i < collisions_.size (); i++)
 758     {
 759       Interval collision_y = collisions_[i].y_;
 760       Real x = collisions_[i].x_;
 761
 762       Real center_beam_y = y_at (x, config);
 763       Interval beam_y = center_beam_y + collisions_[i].beam_y_;
 764
 765       Real dist = infinity_f;
 766       if (!intersection (beam_y, collision_y).is_empty ())
 767         dist = 0.0;
 768       else
 769         dist = min (beam_y.distance (collision_y[DOWN]),
 770                     beam_y.distance (collision_y[UP]));
 771
 772       Real scale_free
 773         = max (parameters.COLLISION_PADDING - dist, 0.0) /
 774           parameters.COLLISION_PADDING;
 775       demerits
 776       += collisions_[i].base_penalty_ *
 777          pow (scale_free, 3) * parameters.COLLISION_PENALTY;
 778     }
 779
 780   config->add (demerits, "C");
 781 }