1 /*
2 This file is part of LilyPond, the GNU music typesetter.
3
4 Copyright (C) 1997--2021 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 "config.hh"
22
23 #include "beam-scoring-problem.hh"
24
25 #include "align-interface.hh"
26 #include "beam.hh"
27 #include "direction.hh"
28 #include "directional-element-interface.hh"
29 #include "grob.hh"
30 #include "grob-array.hh"
31 #include "item.hh"
32 #include "international.hh"
33 #include "interval-minefield.hh"
34 #include "least-squares.hh"
35 #include "libc-extension.hh"
36 #include "note-head.hh"
37 #include "output-def.hh"
38 #include "pointer-group-interface.hh"
39 #include "spanner.hh"
40 #include "staff-symbol-referencer.hh"
41 #include "stencil.hh"
42 #include "stem.hh"
43 #include "warn.hh"
44 #include "string-convert.hh"
45
46 #include <algorithm>
47 #include <cmath>
48 #include <memory>
49 #include <queue>
50 #include <set>
51 #include <vector>
52
53 using std::set;
54 using std::string;
55 using std::unique_ptr;
56 using std::vector;
57
58 // Compute the increase from dr.front () to dr.back ().
59 static constexpr Real
delta(const Drul_array<Real> & dr)60 delta (const Drul_array<Real> &dr)
61 {
62 return dr.back () - dr.front ();
63 }
64
65 Real
get_detail(SCM alist,SCM sym,Real def)66 get_detail (SCM alist, SCM sym, Real def)
67 {
68 SCM entry = scm_assq (sym, alist);
69
70 if (scm_is_pair (entry))
71 return from_scm<double> (scm_cdr (entry), def);
72 return def;
73 }
74
75 void
fill(Grob * him)76 Beam_quant_parameters::fill (Grob *him)
77 {
78 SCM details = get_property (him, "details");
79
80 // General
81 BEAM_EPS = get_detail (details, ly_symbol2scm ("beam-eps"), 1e-3);
82 REGION_SIZE = get_detail (details, ly_symbol2scm ("region-size"), 2);
83
84 // forbidden quants
85 SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0)
86 // For stems that are non-standard, the forbidden beam quanting
87 // doesn't really work, so decrease their importance.
88 * exp (- 8 * fabs (1.0 - from_scm<double> (get_property (him, "length-fraction"), 1.0)));
89 STEM_LENGTH_DEMERIT_FACTOR = get_detail (details, ly_symbol2scm ("stem-length-demerit-factor"), 5);
90 HORIZONTAL_INTER_QUANT_PENALTY = get_detail (details, ly_symbol2scm ("horizontal-inter-quant"), 500);
91
92 STEM_LENGTH_LIMIT_PENALTY = get_detail (details, ly_symbol2scm ("stem-length-limit-penalty"), 5000);
93 DAMPING_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("damping-direction-penalty"), 800);
94 HINT_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("hint-direction-penalty"), 20);
95 MUSICAL_DIRECTION_FACTOR = get_detail (details, ly_symbol2scm ("musical-direction-factor"), 400);
96 IDEAL_SLOPE_FACTOR = get_detail (details, ly_symbol2scm ("ideal-slope-factor"), 10);
97 ROUND_TO_ZERO_SLOPE = get_detail (details, ly_symbol2scm ("round-to-zero-slope"), 0.02);
98
99 // Collisions
100 COLLISION_PENALTY = get_detail (details, ly_symbol2scm ("collision-penalty"), 500);
101
102 /* For grace notes, beams get scaled down to 80%, but glyphs go down
103 to 63% (magstep -4 for accidentals). To make the padding
104 commensurate with glyph size for grace notes, we take the square
105 of the length fraction, yielding a 64% decrease.
106 */
107 COLLISION_PADDING = get_detail (details, ly_symbol2scm ("collision-padding"), 0.5)
108 * sqr (from_scm<double> (get_property (him, "length-fraction"), 1.0));
109 STEM_COLLISION_FACTOR = get_detail (details, ly_symbol2scm ("stem-collision-factor"), 0.1);
110 }
111
112 // Add x if x is positive, add |x|*fac if x is negative.
113 static Real
shrink_extra_weight(Real x,Real fac)114 shrink_extra_weight (Real x, Real fac)
115 {
116 return fabs (x) * ((x < 0) ? fac : 1.0);
117 }
118
119 /****************************************************************/
120
Beam_configuration()121 Beam_configuration::Beam_configuration ()
122 {
123 demerits = 0.0;
124 next_scorer_todo = ORIGINAL_DISTANCE;
125 }
126
done() const127 bool Beam_configuration::done () const
128 {
129 return next_scorer_todo >= NUM_SCORERS;
130 }
131
add(Real demerit,const string & reason)132 void Beam_configuration::add (Real demerit, const string &reason)
133 {
134 demerits += demerit;
135
136 if (demerit)
137 score_card_ += to_string (" %s %.2f", reason.c_str (), demerit);
138 }
139
140 unique_ptr<Beam_configuration>
new_config(Drul_array<Real> start,Drul_array<Real> offset)141 Beam_configuration::new_config (Drul_array<Real> start,
142 Drul_array<Real> offset)
143 {
144 unique_ptr<Beam_configuration> qs (new Beam_configuration);
145 qs->y = Drul_array<Real> (int (start[LEFT]) + offset[LEFT],
146 int (start[RIGHT]) + offset[RIGHT]);
147
148 // This orders the sequence so we try combinations closest to the
149 // the ideal offset first.
150 Real start_score = abs (offset[RIGHT]) + abs (offset[LEFT]);
151 qs->demerits = start_score / 1000.0;
152 qs->next_scorer_todo = ORIGINAL_DISTANCE + 1;
153
154 return qs;
155 }
156
157 Real
y_at(Real x,Beam_configuration const * p) const158 Beam_scoring_problem::y_at (Real x, Beam_configuration const *p) const
159 {
160 return p->y[LEFT] + x * delta (p->y) / x_span_;
161 }
162
163 /****************************************************************/
164
165 /*
166 TODO:
167
168 - Make all demerits customisable
169
170 - Add demerits for quants per se, as to forbid a specific quant
171 entirely
172 */
173
add_collision(Real x,Interval y,Real score_factor)174 void Beam_scoring_problem::add_collision (Real x, Interval y,
175 Real score_factor)
176 {
177 // We used to screen for quant range, but no more.
178
179 Beam_collision c;
180 c.beam_y_.set_empty ();
181
182 for (vsize j = 0; j < segments_.size (); j++)
183 {
184 if (segments_[j].horizontal_.contains (x))
185 c.beam_y_.add_point (segments_[j].vertical_count_ * beam_translation_);
186 if (segments_[j].horizontal_[LEFT] > x)
187 break;
188 }
189 c.beam_y_.widen (0.5 * beam_thickness_);
190
191 c.x_ = x;
192
193 y *= 1 / staff_space_;
194 c.y_ = y;
195 c.base_penalty_ = score_factor;
196 collisions_.push_back (c);
197 }
198
init_instance_variables(Grob * me,Drul_array<Real> ys,bool align_broken_intos)199 void Beam_scoring_problem::init_instance_variables (Grob *me, Drul_array<Real> ys, bool align_broken_intos)
200 {
201 beam_ = dynamic_cast<Spanner *> (me);
202 unquanted_y_ = ys;
203
204 /*
205 If 'ys' are finite, use them as starting points for y-positions of the
206 ends of the beam, instead of the best-fit through the natural ends of
207 the stems. Otherwise, we want to do initial slope calculations.
208 */
209 do_initial_slope_calculations_ = false;
210 for (const auto d : {LEFT, RIGHT})
211 do_initial_slope_calculations_ |= !std::isfinite (unquanted_y_[d]);
212
213 /*
214 Calculations are relative to a unit-scaled staff, i.e. the quants are
215 divided by the current staff_space_.
216 */
217 staff_space_ = Staff_symbol_referencer::staff_space (beam_);
218 beam_thickness_ = Beam::get_beam_thickness (beam_) / staff_space_;
219 line_thickness_ = Staff_symbol_referencer::line_thickness (beam_) / staff_space_;
220 max_beam_count_ = Beam::get_beam_count (beam_);
221 length_fraction_
222 = from_scm<double> (get_property (beam_, "length-fraction"), 1.0);
223 // This is the least-squares DY, corrected for concave beams.
224 musical_dy_ = from_scm<double> (get_property (beam_, "least-squares-dy"), 0);
225
226 vector<Spanner *> beams;
227 align_broken_intos_ = align_broken_intos;
228 if (align_broken_intos_)
229 {
230 Spanner *orig = beam_->original ();
231 if (!orig)
232 align_broken_intos_ = false;
233 else if (!orig->broken_intos_.size ())
234 align_broken_intos_ = false;
235 else
236 beams.insert (beams.end (), orig->broken_intos_.begin (), orig->broken_intos_.end ());
237 }
238 if (!align_broken_intos_)
239 beams.push_back (beam_);
240
241 /*
242 x_span_ is a single scalar, cumulatively summing the length of all the
243 segments the parent beam was broken-into.
244 */
245 x_span_ = 0.0;
246 is_knee_ = false;
247 normal_stem_count_ = 0;
248 for (vsize i = 0; i < beams.size (); i++)
249 {
250 extract_grob_set (beams[i], "stems", stems);
251 extract_grob_set (beams[i], "covered-grobs", fake_collisions);
252 vector<Grob *> collisions;
253
254 for (vsize j = 0; j < fake_collisions.size (); j++)
255 if (fake_collisions[j]->get_system () == beams[i]->get_system ())
256 collisions.push_back (fake_collisions[j]);
257
258 Grob *common[NO_AXES];
259 for (const auto a : {X_AXIS, Y_AXIS})
260 common[a] = common_refpoint_of_array (stems, beams[i], a);
261
262 for (const auto d : {LEFT, RIGHT})
263 common[X_AXIS] = beams[i]->get_bound (d)->common_refpoint (common[X_AXIS], X_AXIS);
264
265 // positions of the endpoints of this beam segment, including any overhangs
266 const Interval x_pos = from_scm (get_property (beams[i], "X-positions"),
267 Interval (0.0, 0.0));
268
269 Drul_array<Grob *> edge_stems (Beam::first_normal_stem (beams[i]),
270 Beam::last_normal_stem (beams[i]));
271
272 Drul_array<bool> dirs_found;
273
274 Real my_y = beams[i]->relative_coordinate (common[Y_AXIS], Y_AXIS);
275
276 Interval beam_width (-1.0, -1.0);
277 for (vsize j = 0; j < stems.size (); j++)
278 {
279 Grob *s = stems[j];
280 beam_multiplicity_.push_back (Stem::beam_multiplicity (stems[j]));
281 head_positions_.push_back (Stem::head_positions (stems[j]));
282 is_normal_.push_back (Stem::is_normal_stem (stems[j]));
283
284 Stem_info si (Stem::get_stem_info (s));
285 si.scale (1 / staff_space_);
286 stem_infos_.push_back (si);
287 chord_start_y_.push_back (Stem::chord_start_y (s));
288 dirs_found[si.dir_] = true;
289
290 Beam_stem_end stem_end
291 = Beam::calc_stem_y (beams[i], s, common,
292 x_pos[LEFT], x_pos[RIGHT], CENTER,
293 Interval (0), 0);
294 Real y = stem_end.stem_y_;
295 /* Remark: French Beaming is irrelevant for beam quanting */
296 base_lengths_.push_back (y / staff_space_);
297 stem_xpositions_.push_back (s->relative_coordinate (common[X_AXIS], X_AXIS) - x_pos[LEFT] + x_span_);
298 stem_ypositions_.push_back (s->relative_coordinate (common[Y_AXIS], Y_AXIS) - my_y);
299
300 if (is_normal_.back ())
301 {
302 if (beam_width[LEFT] == -1.0)
303 beam_width[LEFT] = stem_xpositions_.back ();
304 beam_width[RIGHT] = stem_xpositions_.back ();
305 }
306 }
307
308 edge_dirs_ = {};
309 normal_stem_count_ += Beam::normal_stem_count (beams[i]);
310 if (normal_stem_count_)
311 edge_dirs_ = Drul_array<Direction> (stem_infos_[0].dir_,
312 stem_infos_.back ().dir_);
313
314 is_xstaff_ = has_interface<Align_interface> (common[Y_AXIS]);
315 is_knee_ |= dirs_found[DOWN] && dirs_found[UP];
316
317 staff_radius_ = Staff_symbol_referencer::staff_radius (beams[i]);
318 edge_beam_counts_ = Drul_array<int>
319 (Stem::beam_multiplicity (stems[0]).length () + 1,
320 Stem::beam_multiplicity (stems.back ()).length () + 1);
321
322 // TODO - why are we dividing by staff_space_?
323 beam_translation_ = Beam::get_beam_translation (beams[i]) / staff_space_;
324
325 for (const auto d : {LEFT, RIGHT})
326 {
327 quant_range_[d].set_full ();
328 if (!edge_stems[d])
329 continue;
330
331 Real stem_offset = edge_stems[d]->relative_coordinate (common[Y_AXIS], Y_AXIS)
332 - beams[i]->relative_coordinate (common[Y_AXIS], Y_AXIS);
333 Interval heads = Stem::head_positions (edge_stems[d]) * 0.5 * staff_space_;
334
335 Direction ed = edge_dirs_[d];
336 heads.widen (0.5 * staff_space_
337 + (edge_beam_counts_[d] - 1) * beam_translation_ + beam_thickness_ * .5);
338 quant_range_[d][-ed] = heads[ed] + stem_offset;
339 }
340
341 segments_ = Beam::get_beam_segments (beams[i]);
342 std::sort (segments_.begin (), segments_.end (), beam_segment_less);
343 for (vsize j = 0; j < segments_.size (); j++)
344 segments_[j].horizontal_ += (x_span_ - x_pos[LEFT]);
345
346 set<Grob *> colliding_stems;
347 for (vsize j = 0; j < collisions.size (); j++)
348 {
349 if (!collisions[j]->is_live ())
350 continue;
351
352 if (has_interface<Beam> (collisions[j]) && Beam::is_cross_staff (collisions[j]))
353 continue;
354
355 Box b;
356 for (const auto a : {X_AXIS, Y_AXIS})
357 b[a] = collisions[j]->extent (common[a], a);
358
359 if (b[X_AXIS][RIGHT] < x_pos[LEFT] || b[X_AXIS][LEFT] > x_pos[RIGHT])
360 continue;
361 if (b[X_AXIS].is_empty () || b[Y_AXIS].is_empty ())
362 continue;
363
364 b[X_AXIS] += (x_span_ - x_pos[LEFT]);
365 b[Y_AXIS] -= my_y;
366 Real width = b[X_AXIS].length ();
367 Real width_factor = sqrt (width / staff_space_);
368
369 for (const auto d : {LEFT, RIGHT})
370 add_collision (b[X_AXIS][d], b[Y_AXIS], width_factor);
371
372 Grob *stem = unsmob<Grob> (get_object (collisions[j], "stem"));
373 if (has_interface<Stem> (stem) && Stem::is_normal_stem (stem))
374 {
375 colliding_stems.insert (stem);
376 }
377 }
378
379 for (set<Grob *>::const_iterator it (colliding_stems.begin ()); it != colliding_stems.end (); it++)
380 {
381 Grob *s = *it;
382 Real x = (robust_relative_extent (s, common[X_AXIS], X_AXIS) - x_pos[LEFT] + x_span_).center ();
383
384 Direction stem_dir = get_grob_direction (*it);
385 Interval y;
386 y.set_full ();
387 y[-stem_dir] = Stem::chord_start_y (*it) + (*it)->relative_coordinate (common[Y_AXIS], Y_AXIS)
388 - my_y;
389
390 Real factor = parameters_.STEM_COLLISION_FACTOR;
391 if (!unsmob<Grob> (get_object (s, "beam")))
392 factor = 1.0;
393 add_collision (x, y, factor);
394 }
395 x_span_ += beams[i]->spanner_length ();
396 }
397 }
398
Beam_scoring_problem(Grob * me,Drul_array<Real> ys,bool align_broken_intos)399 Beam_scoring_problem::Beam_scoring_problem (Grob *me, Drul_array<Real> ys, bool align_broken_intos)
400 {
401 beam_ = dynamic_cast<Spanner *> (me);
402 unquanted_y_ = ys;
403 align_broken_intos_ = align_broken_intos;
404
405 parameters_.fill (me);
406 init_instance_variables (me, ys, align_broken_intos);
407 if (do_initial_slope_calculations_)
408 {
409 least_squares_positions ();
410 slope_damping ();
411 shift_region_to_valid ();
412 }
413 }
414
415 // Assuming V is not empty, pick a 'reasonable' point inside V.
416 static Real
point_in_interval(Interval v,Real dist)417 point_in_interval (Interval v, Real dist)
418 {
419 if (std::isinf (v[DOWN]))
420 return v[UP] - dist;
421 else if (std::isinf (v[UP]))
422 return v[DOWN] + dist;
423 else
424 return v.center ();
425 }
426
427 /* Set stem's shorten property if unset.
428
429 TODO:
430 take some y-position (chord/beam/nearest?) into account
431 scmify forced-fraction
432
433 This is done in beam because the shorten has to be uniform over the
434 entire beam.
435 */
436
437 void
set_minimum_dy(Grob * me,Real * dy)438 set_minimum_dy (Grob *me, Real *dy)
439 {
440 if (*dy)
441 {
442 /*
443 If dy is smaller than the smallest quant, we
444 get absurd direction-sign penalties.
445 */
446
447 Real ss = Staff_symbol_referencer::staff_space (me);
448 Real beam_thickness = Beam::get_beam_thickness (me) / ss;
449 Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
450 Real sit = (beam_thickness - slt) / 2;
451 Real inter = 0.5;
452 Real hang = 1.0 - (beam_thickness - slt) / 2;
453
454 *dy = sign (*dy) * std::max (fabs (*dy),
455 std::min (std::min (sit, inter), hang));
456 }
457 }
458
459 void
no_visible_stem_positions()460 Beam_scoring_problem::no_visible_stem_positions ()
461 {
462 if (!head_positions_.size ())
463 {
464 unquanted_y_ = {};
465 return;
466 }
467
468 Interval head_positions;
469 Slice multiplicity;
470 for (vsize i = 0; i < head_positions_.size (); i++)
471 {
472 head_positions.unite (head_positions_[i]);
473 multiplicity.unite (beam_multiplicity_[i]);
474 }
475
476 Direction dir = get_grob_direction (beam_);
477
478 if (!dir)
479 programming_error ("The beam should have a direction by now.");
480
481 Real y = head_positions.linear_combination (dir)
482 * 0.5 * staff_space_
483 + dir * beam_translation_ * (multiplicity.length () + 1);
484
485 unquanted_y_ = Drul_array<Real> (y, y);
486 }
487
488 vsize
first_normal_index()489 Beam_scoring_problem::first_normal_index ()
490 {
491 for (vsize i = 0; i < is_normal_.size (); i++)
492 if (is_normal_[i])
493 return i;
494
495 beam_->programming_error ("No normal stems, but asking for first normal stem index.");
496 return 0;
497 }
498
499 vsize
last_normal_index()500 Beam_scoring_problem::last_normal_index ()
501 {
502 for (vsize i = is_normal_.size (); i--;)
503 if (is_normal_[i])
504 return i;
505
506 beam_->programming_error ("No normal stems, but asking for first normal stem index.");
507 return 0;
508 }
509
510 void
least_squares_positions()511 Beam_scoring_problem::least_squares_positions ()
512 {
513 if (!normal_stem_count_)
514 {
515 no_visible_stem_positions ();
516 return;
517 }
518
519 if (stem_infos_.size () < 1)
520 return;
521
522 vsize fnx = first_normal_index ();
523 vsize lnx = last_normal_index ();
524
525 Drul_array<Real> ideal (stem_infos_[fnx].ideal_y_ + stem_ypositions_[fnx],
526 stem_infos_[lnx].ideal_y_ + stem_ypositions_[lnx]);
527
528 Real y = 0;
529 Real slope = 0;
530 Real dy = 0;
531 Real ldy = 0.0;
532 if (!delta (ideal))
533 {
534 Drul_array<Real> chord (chord_start_y_[0],
535 chord_start_y_.back ());
536
537 /* Simple beams (2 stems) on middle line should be allowed to be
538 slightly sloped.
539
540 However, if both stems reach middle line,
541 ideal[LEFT] == ideal[RIGHT] and delta (ideal) == 0.
542
543 For that case, we apply artificial slope */
544 if (!ideal[LEFT] && delta (chord) && stem_infos_.size () == 2)
545 {
546 const Direction d (delta (chord));
547 unquanted_y_[d] = Beam::get_beam_thickness (beam_) / 2;
548 unquanted_y_[-d] = -unquanted_y_[d];
549 }
550 else
551 unquanted_y_ = ideal;
552
553 ldy = unquanted_y_[RIGHT] - unquanted_y_[LEFT];
554 }
555 else
556 {
557 vector<Offset> ideals;
558 for (vsize i = 0; i < stem_infos_.size (); i++)
559 if (is_normal_[i])
560 ideals.push_back (Offset (stem_xpositions_[i],
561 stem_infos_[i].ideal_y_
562 + stem_ypositions_[i]));
563
564 minimise_least_squares (&slope, &y, ideals);
565
566 dy = slope * x_span_;
567
568 set_minimum_dy (beam_, &dy);
569
570 ldy = dy;
571 unquanted_y_ = Drul_array<Real> (y, (y + dy));
572 }
573
574 musical_dy_ = ldy;
575 set_property (beam_, "least-squares-dy", to_scm (musical_dy_));
576 }
577
578 /*
579 Determine whether a beam is concave.
580
581 A beam is concave when the middle notes get closer to the
582 beam than the left and right edge notes.
583
584 This is determined in two ways: by looking at the positions of the
585 middle notes, or by looking at the deviation of the inside notes
586 compared to the line connecting first and last.
587
588 The tricky thing is what to do with beams with chords. There are no
589 real guidelines in this case.
590 */
591
592 bool
is_concave_single_notes(vector<int> const & positions,Direction beam_dir)593 is_concave_single_notes (vector<int> const &positions, Direction beam_dir)
594 {
595 Interval covering;
596 covering.add_point (positions[0]);
597 covering.add_point (positions.back ());
598
599 bool above = false;
600 bool below = false;
601 bool concave = false;
602
603 /*
604 notes above and below the interval covered by 1st and last note.
605 */
606 for (vsize i = 1; i + 1 < positions.size (); i++)
607 {
608 above = above || (positions[i] > covering[UP]);
609 below = below || (positions[i] < covering[DOWN]);
610 }
611
612 concave = concave || (above && below);
613 /*
614 A note as close or closer to the beam than begin and end, but the
615 note is reached in the opposite direction as the last-first dy
616 */
617 int dy = positions.back () - positions[0];
618 int closest = std::max (beam_dir * positions.back (), beam_dir * positions[0]);
619 for (vsize i = 2; !concave && i + 1 < positions.size (); i++)
620 {
621 int inner_dy = positions[i] - positions[i - 1];
622 if (sign (inner_dy) != sign (dy)
623 && (beam_dir * positions[i] >= closest
624 || beam_dir * positions[i - 1] >= closest))
625 concave = true;
626 }
627
628 bool all_closer = true;
629 for (vsize i = 1; all_closer && i + 1 < positions.size (); i++)
630 {
631 all_closer = all_closer
632 && (beam_dir * positions[i] > closest);
633 }
634
635 concave = concave || all_closer;
636 return concave;
637 }
638
639 Real
calc_positions_concaveness(vector<int> const & positions,Direction beam_dir)640 calc_positions_concaveness (vector<int> const &positions, Direction beam_dir)
641 {
642 Real dy = positions.back () - positions[0];
643 Real slope = dy / static_cast<Real> (positions.size () - 1);
644 Real concaveness = 0.0;
645 for (vsize i = 1; i + 1 < positions.size (); i++)
646 {
647 Real line_y = slope * static_cast<Real> (i) + positions[0];
648 concaveness += std::max (beam_dir * (positions[i] - line_y), 0.0);
649 }
650
651 concaveness /= static_cast<Real> (positions.size ());
652
653 /*
654 Normalize. For dy = 0, the slope ends up as 0 anyway, so the
655 scaling of concaveness doesn't matter much.
656 */
657 if (dy)
658 concaveness /= fabs (dy);
659 return concaveness;
660 }
661
662 Real
calc_concaveness()663 Beam_scoring_problem::calc_concaveness ()
664 {
665 SCM conc = get_property (beam_, "concaveness");
666 if (scm_is_number (conc))
667 return scm_to_double (conc);
668
669 if (is_knee_ || is_xstaff_)
670 return 0.0;
671
672 Direction beam_dir = CENTER;
673 for (vsize i = is_normal_.size (); i--;)
674 if (is_normal_[i] && stem_infos_[i].dir_)
675 beam_dir = stem_infos_[i].dir_;
676
677 if (normal_stem_count_ <= 2)
678 return 0.0;
679
680 vector<int> close_positions;
681 vector<int> far_positions;
682 for (vsize i = 0; i < is_normal_.size (); i++)
683 if (is_normal_[i])
684 {
685 /*
686 For chords, we take the note head that is closest to the beam.
687
688 Hmmm.. wait, for the beams in the last measure of morgenlied,
689 this doesn't look so good. Let's try the heads farthest from
690 the beam.
691 */
692
693 auto close_pos = static_cast<int> (rint (head_positions_[i][beam_dir]));
694 close_positions.push_back (close_pos);
695 auto far_pos = static_cast<int> (rint (head_positions_[i][-beam_dir]));
696 far_positions.push_back (far_pos);
697 }
698
699 Real concaveness = 0.0;
700
701 if (is_concave_single_notes (beam_dir == UP ? close_positions : far_positions, beam_dir))
702 {
703 concaveness = 10000;
704 }
705 else
706 {
707 concaveness = (calc_positions_concaveness (far_positions, beam_dir)
708 + calc_positions_concaveness (close_positions, beam_dir)) / 2;
709 }
710
711 return concaveness;
712 }
713
714 void
slope_damping()715 Beam_scoring_problem::slope_damping ()
716 {
717 if (normal_stem_count_ <= 1)
718 return;
719
720 SCM s = get_property (beam_, "damping");
721 Real damping = scm_to_double (s);
722 Real concaveness = calc_concaveness ();
723 if ((concaveness >= 10000) || (damping >= 10000))
724 {
725 unquanted_y_[LEFT] = unquanted_y_[RIGHT];
726 musical_dy_ = 0;
727 damping = 0;
728 }
729
730 if ((damping) && (damping + concaveness))
731 {
732 Real dy = unquanted_y_[RIGHT] - unquanted_y_[LEFT];
733
734 Real slope = dy && x_span_ ? dy / x_span_ : 0;
735
736 slope = 0.6 * tanh (slope) / (damping + concaveness);
737
738 Real damped_dy = slope * x_span_;
739
740 set_minimum_dy (beam_, &damped_dy);
741
742 unquanted_y_[LEFT] += (dy - damped_dy) / 2;
743 unquanted_y_[RIGHT] -= (dy - damped_dy) / 2;
744 }
745 }
746
747 void
shift_region_to_valid()748 Beam_scoring_problem::shift_region_to_valid ()
749 {
750 if (!normal_stem_count_)
751 return;
752
753 Real beam_dy = unquanted_y_[RIGHT] - unquanted_y_[LEFT];
754 Real slope = x_span_ ? beam_dy / x_span_ : 0.0;
755
756 /*
757 Shift the positions so that we have a chance of finding good
758 quants (i.e. no short stem failures.)
759 */
760 Interval feasible_left_point;
761 feasible_left_point.set_full ();
762
763 for (vsize i = 0; i < stem_infos_.size (); i++)
764 {
765 // TODO - check for invisible here...
766 Real left_y
767 = stem_infos_[i].shortest_y_
768 - slope * stem_xpositions_ [i];
769
770 /*
771 left_y is now relative to the stem S. We want relative to
772 ourselves, so translate:
773 */
774 left_y += stem_ypositions_[i];
775 Interval flp;
776 flp.set_full ();
777 flp[-stem_infos_[i].dir_] = left_y;
778
779 feasible_left_point.intersect (flp);
780 }
781
782 vector<Grob *> filtered;
783 /*
784 We only update these for objects that are too large for quanting
785 to find a workaround. Typically, these are notes with
786 stems, and timesig/keysig/clef, which take out the entire area
787 inside the staff as feasible.
788
789 The code below disregards the thickness and multiplicity of the
790 beam. This should not be a problem, as the beam quanting will
791 take care of computing the impact those exactly.
792 */
793 Real min_y_size = 2.0;
794
795 // A list of intervals into which beams may not fall
796 vector<Interval> forbidden_intervals;
797
798 for (vsize i = 0; i < collisions_.size (); i++)
799 {
800 if (collisions_[i].x_ < 0 || collisions_[i].x_ > x_span_)
801 continue;
802
803 if (collisions_[i].y_.length () < min_y_size)
804 continue;
805
806 Real dy = slope * collisions_[i].x_;
807
808 Interval disallowed;
809 for (const auto yd : {DOWN, UP})
810 {
811 Real left_y = collisions_[i].y_[yd] - dy;
812 disallowed[yd] = left_y;
813 }
814
815 forbidden_intervals.push_back (disallowed);
816 }
817
818 std::sort (forbidden_intervals.begin (), forbidden_intervals.end (),
819 Interval::left_less);
820 Real beam_left_y = unquanted_y_[LEFT];
821 Interval feasible_beam_placements (beam_left_y, beam_left_y);
822
823 Interval_minefield minefield (feasible_beam_placements, 0.0);
824 for (vsize i = 0; i < forbidden_intervals.size (); i++)
825 minefield.add_forbidden_interval (forbidden_intervals[i]);
826 minefield.solve ();
827 feasible_beam_placements = minefield.feasible_placements ();
828
829 // if the beam placement falls out of the feasible region, we push it
830 // to infinity so that it can never be a feasible candidate below
831 for (const auto d : {DOWN, UP})
832 {
833 if (!feasible_left_point.contains (feasible_beam_placements[d]))
834 feasible_beam_placements[d] = d * infinity_f;
835 }
836
837 if ((feasible_beam_placements[UP] == infinity_f && feasible_beam_placements[DOWN] == -infinity_f) && !feasible_left_point.is_empty ())
838 {
839 // We are somewhat screwed: we have a collision, but at least
840 // there is a way to satisfy stem length constraints.
841 beam_left_y = point_in_interval (feasible_left_point, 2.0);
842 }
843 else if (!feasible_left_point.is_empty ())
844 {
845 // Only one of them offers is feasible solution. Pick that one.
846 if (abs (beam_left_y - feasible_beam_placements[DOWN]) > abs (beam_left_y - feasible_beam_placements[UP]))
847 beam_left_y = feasible_beam_placements[UP];
848 else
849 beam_left_y = feasible_beam_placements[DOWN];
850 }
851 else
852 {
853 // We are completely screwed.
854 beam_->warning (_ ("no viable initial configuration found: may not find good beam slope"));
855 }
856
857 unquanted_y_ = Drul_array<Real> (beam_left_y, (beam_left_y + beam_dy));
858 }
859
860 void
generate_quants(vector<unique_ptr<Beam_configuration>> * scores) const861 Beam_scoring_problem::generate_quants (vector<unique_ptr<Beam_configuration>> *scores) const
862 {
863 auto region_size = static_cast<int> (parameters_.REGION_SIZE);
864
865 // Knees and collisions are harder, lets try some more possibilities
866 if (is_knee_)
867 region_size += 2;
868 if (collisions_.size ())
869 region_size += 2;
870
871 Real straddle = 0.0;
872 Real sit = (beam_thickness_ - line_thickness_) / 2;
873 Real inter = 0.5;
874 Real hang = 1.0 - (beam_thickness_ - line_thickness_) / 2;
875 Real base_quants [] = {straddle, sit, inter, hang};
876 int num_base_quants = int (sizeof (base_quants) / sizeof (Real));
877
878 /* for normal-sized beams, in case of more than 4 beams, the outer beam
879 used for generating quants will never interfere with staff lines, but
880 prevent the inside-staff beams from being neatly positioned.
881 A correctional grid_shift has to be applied to compensate. */
882 Real grid_shift = 0.0;
883 /* grid shift only makes sense for widened normal-sized beams: */
884 if (!is_knee_ && max_beam_count_ > 4 && length_fraction_ == 1.0)
885 grid_shift = (max_beam_count_ - 4) * (1.0 - beam_translation_);
886
887 /*
888 Asymetry ? should run to <= region_size ?
889 */
890 vector<Real> unshifted_quants;
891 for (int i = -region_size; i < region_size; i++)
892 for (int j = 0; j < num_base_quants; j++)
893 {
894 unshifted_quants.push_back (i + base_quants[j]);
895 }
896
897 for (vsize i = 0; i < unshifted_quants.size (); i++)
898 for (vsize j = 0; j < unshifted_quants.size (); j++)
899 {
900 Interval corr (0.0, 0.0);
901 if (grid_shift)
902 for (const auto d : {LEFT, RIGHT})
903 /* apply grid shift if quant outside 5-line staff: */
904 if ((unquanted_y_[d] + unshifted_quants[i]) * edge_dirs_[d] > 2.5)
905 corr[d] = grid_shift * edge_dirs_[d];
906 auto c = Beam_configuration::new_config
907 (unquanted_y_,
908 Drul_array<Real> (unshifted_quants[i] - corr[LEFT],
909 unshifted_quants[j] - corr[RIGHT]));
910
911 for (const auto d : {LEFT, RIGHT})
912 {
913 if (!quant_range_[d].contains (c->y[d]))
914 {
915 c.reset ();
916 break;
917 }
918 }
919 if (c)
920 scores->push_back (std::move (c));
921 }
922
923 }
924
one_scorer(Beam_configuration * config) const925 void Beam_scoring_problem::one_scorer (Beam_configuration *config) const
926 {
927 switch (config->next_scorer_todo)
928 {
929 case SLOPE_IDEAL:
930 score_slope_ideal (config);
931 break;
932 case SLOPE_DIRECTION:
933 score_slope_direction (config);
934 break;
935 case SLOPE_MUSICAL:
936 score_slope_musical (config);
937 break;
938 case FORBIDDEN:
939 score_forbidden_quants (config);
940 break;
941 case STEM_LENGTHS:
942 score_stem_lengths (config);
943 break;
944 case COLLISIONS:
945 score_collisions (config);
946 break;
947 case HORIZONTAL_INTER:
948 score_horizontal_inter_quants (config);
949 break;
950
951 case NUM_SCORERS:
952 case ORIGINAL_DISTANCE:
953 default:
954 assert (false);
955 }
956 config->next_scorer_todo++;
957 }
958
959 Beam_configuration *
force_score(SCM inspect_quants,const vector<unique_ptr<Beam_configuration>> & configs) const960 Beam_scoring_problem::force_score (SCM inspect_quants,
961 const vector<unique_ptr<Beam_configuration>> &configs) const
962 {
963 Drul_array<Real> ins = from_scm<Drul_array<Real>> (inspect_quants);
964 Real mindist = 1e6;
965 Beam_configuration *best = NULL;
966 for (vsize i = 0; i < configs.size (); i++)
967 {
968 Real d = fabs (configs[i]->y[LEFT] - ins[LEFT]) + fabs (configs[i]->y[RIGHT] - ins[RIGHT]);
969 if (d < mindist)
970 {
971 best = configs[i].get ();
972 mindist = d;
973 }
974 }
975 if (mindist > 1e5)
976 programming_error ("cannot find quant");
977
978 while (!best->done ())
979 one_scorer (best);
980
981 return best;
982 }
983
984 Drul_array<Real>
solve() const985 Beam_scoring_problem::solve () const
986 {
987 vector<unique_ptr<Beam_configuration>> configs;
988 generate_quants (&configs);
989
990 if (configs.empty ())
991 {
992 programming_error ("No viable beam quanting found. Using unquanted y value.");
993 return unquanted_y_;
994 }
995
996 if (from_scm<bool> (get_property (beam_, "skip-quanting")))
997 return unquanted_y_;
998
999 Beam_configuration *best = NULL;
1000
1001 bool debug
1002 = from_scm<bool> (beam_->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring")));
1003 SCM inspect_quants = get_property (beam_, "inspect-quants");
1004 if (scm_is_pair (inspect_quants))
1005 {
1006 debug = true;
1007 best = force_score (inspect_quants, configs);
1008 }
1009 else
1010 {
1011 std::priority_queue < Beam_configuration *, std::vector<Beam_configuration *>,
1012 Beam_configuration_less > queue;
1013 for (vsize i = 0; i < configs.size (); i++)
1014 queue.push (configs[i].get ());
1015
1016 /*
1017 TODO
1018
1019 It would be neat if we generated new configurations on the
1020 fly, depending on the best complete score so far, eg.
1021
1022 if (best->done()) {
1023 if (best->demerits < sqrt(queue.size())
1024 break;
1025 while (best->demerits > sqrt(queue.size()) {
1026 generate and insert new configuration
1027 }
1028 }
1029
1030 that would allow us to do away with region_size altogether.
1031 */
1032 while (true)
1033 {
1034 best = queue.top ();
1035 if (best->done ())
1036 break;
1037
1038 queue.pop ();
1039 one_scorer (best);
1040 queue.push (best);
1041 }
1042 }
1043
1044 Drul_array<Real> final_positions = best->y;
1045
1046 if (debug)
1047 {
1048 // debug quanting
1049 int completed = 0;
1050 for (vsize i = 0; i < configs.size (); i++)
1051 {
1052 if (configs[i]->done ())
1053 completed++;
1054 }
1055
1056 string card = best->score_card_ + to_string (" c%d/%zu", completed, configs.size ());
1057 set_property (beam_, "annotation", ly_string2scm (card));
1058 }
1059
1060 configs.clear ();
1061 if (align_broken_intos_)
1062 {
1063 Interval normalized_endpoints = from_scm (get_property (beam_, "normalized-endpoints"), Interval (0, 1));
1064 Real y_length = final_positions[RIGHT] - final_positions[LEFT];
1065
1066 final_positions[LEFT] += normalized_endpoints[LEFT] * y_length;
1067 final_positions[RIGHT] -= (1 - normalized_endpoints[RIGHT]) * y_length;
1068 }
1069
1070 return final_positions;
1071 }
1072
1073 void
score_stem_lengths(Beam_configuration * config) const1074 Beam_scoring_problem::score_stem_lengths (Beam_configuration *config) const
1075 {
1076 Real limit_penalty = parameters_.STEM_LENGTH_LIMIT_PENALTY;
1077 Drul_array<Real> score;
1078 Drul_array<int> count;
1079
1080 for (vsize i = 0; i < stem_xpositions_.size (); i++)
1081 {
1082 if (!is_normal_[i])
1083 continue;
1084
1085 Real x = stem_xpositions_[i];
1086 Real dx = x_span_;
1087 Real beam_y = dx
1088 ? config->y[RIGHT] * x / dx + config->y[LEFT] * (x_span_ - x) / dx
1089 : (config->y[RIGHT] + config->y[LEFT]) / 2;
1090 Real current_y = beam_y + base_lengths_[i];
1091 Real length_pen = parameters_.STEM_LENGTH_DEMERIT_FACTOR;
1092
1093 Stem_info info = stem_infos_[i];
1094 Direction d = info.dir_;
1095
1096 score[d] += limit_penalty * std::max (0.0, (d * (info.shortest_y_ - current_y)));
1097
1098 Real ideal_diff = d * (current_y - info.ideal_y_);
1099 Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
1100
1101 /* We introduce a power, to make the scoring strictly
1102 convex. Otherwise a symmetric knee beam (up/down/up/down)
1103 does not have an optimum in the middle. */
1104 if (is_knee_)
1105 ideal_score = pow (ideal_score, 1.1);
1106
1107 score[d] += length_pen * ideal_score;
1108 count[d]++;
1109 }
1110
1111 /* Divide by number of stems, to make the measure scale-free. */
1112 for (const auto d : {DOWN, UP})
1113 score[d] /= std::max (count[d], 1);
1114
1115 /*
1116 sometimes, two perfectly symmetric kneed beams will have the same score
1117 and can either be quanted up or down.
1118
1119 we choose the quanting in the direction of the slope so that the first stem
1120 always seems longer, reaching to the second, rather than squashed.
1121 */
1122 if (is_knee_ && (count[LEFT] == count[RIGHT]) && (count[LEFT] == 1))
1123 {
1124 const Direction d (delta (unquanted_y_));
1125 if (d)
1126 score[d] += (score[d] < 1.0) ? 0.01 : 0.0;
1127 }
1128
1129 config->add (score[LEFT] + score[RIGHT], "L");
1130 }
1131
1132 void
score_slope_direction(Beam_configuration * config) const1133 Beam_scoring_problem::score_slope_direction (Beam_configuration *config) const
1134 {
1135 Real dy = delta (config->y);
1136 Real damped_dy = delta (unquanted_y_);
1137 Real dem = 0.0;
1138 /*
1139 DAMPING_DIRECTION_PENALTY is a very harsh measure, while for
1140 complex beaming patterns, horizontal is often a good choice.
1141
1142 TODO: find a way to incorporate the complexity of the beam in this
1143 penalty.
1144 */
1145 if (sign (damped_dy) != sign (dy))
1146 {
1147 if (!dy)
1148 {
1149 if (fabs (damped_dy / x_span_) > parameters_.ROUND_TO_ZERO_SLOPE)
1150 dem += parameters_.DAMPING_DIRECTION_PENALTY;
1151 else
1152 dem += parameters_.HINT_DIRECTION_PENALTY;
1153 }
1154 else
1155 dem += parameters_.DAMPING_DIRECTION_PENALTY;
1156 }
1157
1158 config->add (dem, "Sd");
1159 }
1160
1161 // Score for going against the direction of the musical pattern
1162 void
score_slope_musical(Beam_configuration * config) const1163 Beam_scoring_problem::score_slope_musical (Beam_configuration *config) const
1164 {
1165 Real dy = delta (config->y);
1166 Real dem = parameters_.MUSICAL_DIRECTION_FACTOR
1167 * std::max (0.0, (fabs (dy) - fabs (musical_dy_)));
1168 config->add (dem, "Sm");
1169 }
1170
1171 // Score deviation from calculated ideal slope.
1172 void
score_slope_ideal(Beam_configuration * config) const1173 Beam_scoring_problem::score_slope_ideal (Beam_configuration *config) const
1174 {
1175 Real dy = delta (config->y);
1176 Real damped_dy = delta (unquanted_y_);
1177 Real dem = 0.0;
1178
1179 Real slope_penalty = parameters_.IDEAL_SLOPE_FACTOR;
1180
1181 /* Xstaff beams tend to use extreme slopes to get short stems. We
1182 put in a penalty here. */
1183 if (is_xstaff_)
1184 slope_penalty *= 10;
1185
1186 /* Huh, why would a too steep beam be better than a too flat one ? */
1187 dem += shrink_extra_weight (fabs (damped_dy) - fabs (dy), 1.5)
1188 * slope_penalty;
1189
1190 config->add (dem, "Si");
1191 }
1192
1193 static Real
my_modf(Real x)1194 my_modf (Real x)
1195 {
1196 return x - floor (x);
1197 }
1198
1199 // TODO - there is some overlap with forbidden quants, but for
1200 // horizontal beams, it is much more serious to have stafflines
1201 // appearing in the wrong place, so we have a separate scorer.
1202 void
score_horizontal_inter_quants(Beam_configuration * config) const1203 Beam_scoring_problem::score_horizontal_inter_quants (Beam_configuration *config) const
1204 {
1205 if (delta (config->y) == 0.0
1206 && abs (config->y[LEFT]) < staff_radius_ * staff_space_)
1207 {
1208 Real yshift = config->y[LEFT] - 0.5 * staff_space_;
1209 if (fabs (round_halfway_up (yshift) - yshift) < 0.01 * staff_space_)
1210 config->add (parameters_.HORIZONTAL_INTER_QUANT_PENALTY, "H");
1211 }
1212 }
1213
1214 /*
1215 TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
1216 because for 32nd and 64th beams the forbidden quants are relatively
1217 more important than stem lengths.
1218 */
1219 void
score_forbidden_quants(Beam_configuration * config) const1220 Beam_scoring_problem::score_forbidden_quants (Beam_configuration *config) const
1221 {
1222 Real dy = delta (config->y);
1223
1224 Real extra_demerit
1225 = parameters_.SECONDARY_BEAM_DEMERIT
1226 / std::max (edge_beam_counts_[LEFT], edge_beam_counts_[RIGHT]);
1227
1228 Real dem = 0.0;
1229 Real eps = parameters_.BEAM_EPS;
1230
1231 for (const auto d : {LEFT, RIGHT})
1232 {
1233 for (int j = 1; j <= edge_beam_counts_[d]; j++)
1234 {
1235 Direction stem_dir = edge_dirs_[d];
1236
1237 /*
1238 The fudge_factor is to provide a little leniency for
1239 borderline cases. If we do 2.0, then the upper outer line
1240 will be in the gap of the (2, sit) quant, leading to a
1241 false demerit. By increasing the fudge factor to 2.2, we
1242 fix this case.
1243 */
1244 Real fudge_factor = 2.2;
1245 Real gap1 = config->y[d] - stem_dir * ((j - 1) * beam_translation_ + beam_thickness_ / 2 - line_thickness_ / fudge_factor);
1246 Real gap2 = config->y[d] - stem_dir * (j * beam_translation_ - beam_thickness_ / 2 + line_thickness_ / fudge_factor);
1247
1248 Interval gap;
1249 gap.add_point (gap1);
1250 gap.add_point (gap2);
1251
1252 for (Real k = -staff_radius_;
1253 k <= staff_radius_ + eps; k += 1.0)
1254 if (gap.contains (k))
1255 {
1256 Real dist = std::min (fabs (gap[UP] - k), fabs (gap[DOWN] - k));
1257
1258 /*
1259 this parameter is tuned to grace-stem-length.ly
1260 retuned from 0.40 to 0.39 by MS because of slight increases
1261 in gap.length () resulting from measuring beams at real ends
1262 instead of from the middle of stems.
1263
1264 TODO:
1265 This function needs better comments so we know what is forbidden
1266 and why.
1267 */
1268 Real fixed_demerit = 0.39;
1269
1270 dem += extra_demerit
1271 * (fixed_demerit
1272 + (1 - fixed_demerit) * (dist / gap.length ()) * 2);
1273 }
1274 }
1275 }
1276
1277 config->add (dem, "Fl");
1278 dem = 0.0;
1279 if (std::max (edge_beam_counts_[LEFT], edge_beam_counts_[RIGHT]) >= 2)
1280 {
1281 Real straddle = 0.0;
1282 Real sit = (beam_thickness_ - line_thickness_) / 2;
1283 Real inter = 0.5;
1284 Real hang = 1.0 - (beam_thickness_ - line_thickness_) / 2;
1285
1286 for (const auto d : {LEFT, RIGHT})
1287 {
1288 if (edge_beam_counts_[d] >= 2
1289 && fabs (config->y[d] - edge_dirs_[d] * beam_translation_) < staff_radius_ + inter)
1290 {
1291 // TODO up/down symmetry.
1292 if (edge_dirs_[d] == UP && dy <= eps
1293 && fabs (my_modf (config->y[d]) - sit) < eps)
1294 dem += extra_demerit;
1295
1296 if (edge_dirs_[d] == DOWN && dy >= eps
1297 && fabs (my_modf (config->y[d]) - hang) < eps)
1298 dem += extra_demerit;
1299 }
1300
1301 if (edge_beam_counts_[d] >= 3
1302 && fabs (config->y[d] - 2 * edge_dirs_[d] * beam_translation_) < staff_radius_ + inter)
1303 {
1304 // TODO up/down symmetry.
1305 if (edge_dirs_[d] == UP && dy <= eps
1306 && fabs (my_modf (config->y[d]) - straddle) < eps)
1307 dem += extra_demerit;
1308
1309 if (edge_dirs_[d] == DOWN && dy >= eps
1310 && fabs (my_modf (config->y[d]) - straddle) < eps)
1311 dem += extra_demerit;
1312 }
1313 }
1314 }
1315
1316 config->add (dem, "Fs");
1317 }
1318
1319 void
score_collisions(Beam_configuration * config) const1320 Beam_scoring_problem::score_collisions (Beam_configuration *config) const
1321 {
1322 Real demerits = 0.0;
1323 for (vsize i = 0; i < collisions_.size (); i++)
1324 {
1325 Interval collision_y = collisions_[i].y_;
1326 Real x = collisions_[i].x_;
1327
1328 Real center_beam_y = y_at (x, config);
1329 Interval beam_y = center_beam_y + collisions_[i].beam_y_;
1330
1331 Real dist = infinity_f;
1332 if (!intersection (beam_y, collision_y).is_empty ())
1333 dist = 0.0;
1334 else
1335 dist = std::min (beam_y.distance (collision_y[DOWN]),
1336 beam_y.distance (collision_y[UP]));
1337
1338 Real scale_free
1339 = std::max (parameters_.COLLISION_PADDING - dist, 0.0)
1340 / parameters_.COLLISION_PADDING;
1341 Real collision_demerit = collisions_[i].base_penalty_ *
1342 pow (scale_free, 3) * parameters_.COLLISION_PENALTY;
1343
1344 if (collision_demerit > 0)
1345 {
1346 demerits += collision_demerit;
1347 }
1348 }
1349
1350 config->add (demerits, "C");
1351 }
1352