lily/spring-spacer.cc

   1 /*
   2   spring-spacer.cc -- implement Spring_spacer
   3
   4   source file of the GNU LilyPond music typesetter
   5
   6   (c) 1996--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
   7 */
   8
   9
  10 #include <math.h>
  11 #include <limits.h>
  12 #include "killing-cons.tcc"
  13 #include "spring-spacer.hh"
  14 #include "p-col.hh"
  15 #include "debug.hh"
  16 #include "dimensions.hh"
  17 #include "qlp.hh"
  18 #include "unionfind.hh"
  19 #include "idealspacing.hh"
  20 #include "pointer.tcc"
  21 #include "score-column.hh"
  22 #include "paper-def.hh"
  23 #include "colhpos.hh"
  24 #include "main.hh"
  25
  26 Vector
  27 Spring_spacer::default_solution() const
  28 {
  29   return try_initial_solution();
  30 }
  31
  32 Score_column*
  33 Spring_spacer::scol_l (int i)
  34 {
  35   return dynamic_cast<Score_column*>(cols_[i].pcol_l_);
  36 }
  37
  38 const Real COLFUDGE=1e-3;
  39 template class P<Real>;         // ugh.
  40
  41 bool
  42 Spring_spacer::contains_b (Paper_column const *w)
  43 {
  44   for (int i=0; i< cols_.size(); i++)
  45     if (cols_[i].pcol_l_ == w)
  46       return true;
  47   return false;
  48 }
  49
  50
  51 void
  52 Spring_spacer::OK() const
  53 {
  54 #ifndef NDEBUG
  55   for (int i = 1; i < cols_.size(); i++)
  56     assert (cols_[i].rank_i_ > cols_[i-1].rank_i_);
  57   for (int i = 1; i < loose_col_arr_.size(); i++)
  58     assert (loose_col_arr_[i].rank_i_ > loose_col_arr_[i-1].rank_i_);
  59 #endif
  60 }
  61
  62 /**
  63   Make sure no unconnected columns happen.
  64  */
  65 void
  66 Spring_spacer::handle_loose_cols()
  67 {
  68   Union_find connected (cols_.size());
  69   Array<int> fixed;
  70
  71   for (Cons<Idealspacing> *i  = ideal_p_list_; i; i = i->next_)
  72     {
  73       connected.connect (i->car_->cols_drul_[LEFT],i->car_->cols_drul_[RIGHT]);
  74     }
  75   for (int i = 0; i < cols_.size(); i++)
  76     if (cols_[i].fixed_b())
  77       fixed.push (i);
  78   for (int i=1; i < fixed.size(); i++)
  79     connected.connect (fixed[i-1], fixed[i]);
  80
  81   for (int i = cols_.size(); i--;)
  82     {
  83       if (! connected.equiv (fixed[0], i))
  84         {
  85           warning (_f ("unconnected column: %d", i));
  86           loosen_column (i);
  87         }
  88     }
  89   OK();
  90 }
  91
  92
  93 /**
  94   Guess a stupid position for loose columns.  Put loose columns at
  95   regular distances from enclosing calced columns
  96   */
  97 void
  98 Spring_spacer::position_loose_cols (Vector &sol_vec) const
  99 {
 100   if (!loose_col_arr_.size())
 101     return ;
 102   assert (sol_vec.dim());
 103   Array<bool> fix_b_arr;
 104   fix_b_arr.set_size (cols_.size() + loose_col_arr_.size ());
 105   Real utter_right_f=-infinity_f;
 106   Real utter_left_f =infinity_f;
 107   for (int i=0; i < loose_col_arr_.size(); i++)
 108     {
 109       fix_b_arr[loose_col_arr_[i].rank_i_] = false;
 110     }
 111   for (int i=0; i < cols_.size(); i++)
 112     {
 113       int r= cols_[i].rank_i_;
 114       fix_b_arr[r] = true;
 115       utter_right_f = utter_right_f >? sol_vec (i);
 116       utter_left_f = utter_left_f <? sol_vec (i);
 117     }
 118   Vector v (fix_b_arr.size());
 119   int j =0;
 120   int k =0;
 121   for (int i=0; i < v.dim(); i++)
 122     {
 123       if (fix_b_arr[i])
 124         {
 125           assert (cols_[j].rank_i_ == i);
 126           v (i) = sol_vec (j++);
 127         }
 128       else
 129         {
 130           Real left_pos_f =
 131             (j>0) ?sol_vec (j-1) : utter_left_f;
 132           Real right_pos_f =
 133             (j < sol_vec.dim()) ? sol_vec (j) : utter_right_f;
 134           int left_rank = (j>0) ? cols_[j-1].rank_i_ : 0;
 135           int right_rank = (j<sol_vec.dim()) ? cols_[j].rank_i_ : sol_vec.dim ();
 136
 137           int d_r = right_rank - left_rank;
 138           Column_info loose=loose_col_arr_[k++];
 139           int r = loose.rank_i_ ;
 140           assert (r > left_rank && r < right_rank);
 141
 142           v (i) =  (r - left_rank)*left_pos_f/ d_r +
 143             (right_rank - r) *right_pos_f /d_r;
 144         }
 145     }
 146   sol_vec = v;
 147 }
 148
 149 bool
 150 Spring_spacer::check_constraints (Vector v) const
 151 {
 152   int dim=v.dim();
 153   assert (dim == cols_.size());
 154   DOUT << "checking " << v;
 155   for (int i=0; i < dim; i++)
 156     {
 157       if (cols_[i].fixed_b() &&
 158           abs (cols_[i].fixed_position() - v (i)) > COLFUDGE)
 159         {
 160           DOUT << "Fixpos broken\n";
 161           return false;
 162         }
 163       Array<Spacer_rod> const &rods (cols_[i].rods_[RIGHT]);
 164       for (int j =0; j < rods.size (); j++)
 165         {
 166           int other =rods[j].other_idx_;
 167           Real diff =v (other) - v (i) ;
 168           if (COLFUDGE +diff <  rods[j].distance_f_)
 169             {
 170               DOUT << "i, other_i: " << i << "  " << other << '\n';
 171               DOUT << "dist, minimal = " << diff << " "
 172                    << rods[j].distance_f_ << '\n';
 173               return false;
 174             }
 175         }
 176
 177     }
 178   return true;
 179 }
 180
 181 /** try to generate a solution which obeys the min
 182     distances and fixed positions
 183  */
 184 Vector
 185 Spring_spacer::try_initial_solution() const
 186 {
 187   Vector v;
 188   if (!try_initial_solution_and_tell (v))
 189     {
 190       warning (_ ("I'm too fat; call Oprah"));
 191     }
 192   return v;
 193
 194 }
 195
 196 bool
 197 Spring_spacer::try_initial_solution_and_tell (Vector &v) const
 198 {
 199   int dim=cols_.size();
 200   bool succeeded = true;
 201   Vector initsol (dim);
 202
 203   assert (cols_[0].fixed_b ());
 204   DOUT << "fixpos 0 " << cols_[0].fixed_position ();
 205   for (int i=0; i < dim; i++)
 206     {
 207       Real min_x = i ?  initsol (i-1) : cols_[0].fixed_position ();
 208       Array<Spacer_rod> const &sr_arr(cols_[i].rods_[LEFT]);
 209       for (int j=0; j < sr_arr.size (); j++)
 210         {
 211           min_x = min_x >? (initsol (sr_arr[j].other_idx_) + sr_arr[j].distance_f_);
 212         }
 213       initsol (i) = min_x;
 214
 215       if (cols_[i].fixed_b())
 216         {
 217           initsol (i)=cols_[i].fixed_position();
 218           if (initsol (i) < min_x )
 219             {
 220               DOUT << "failing: init, min : " << initsol (i) << " " << min_x << '\n';
 221               initsol (i) = min_x;
 222               succeeded = false;
 223             }
 224         }
 225     }
 226   v = initsol;
 227
 228   DOUT << "tried and told solution: " << v;
 229   if (!succeeded)
 230     DOUT << "(failed)\n";
 231   return succeeded;
 232 }
 233
 234
 235
 236 // generate the matrices
 237 void
 238 Spring_spacer::make_matrices (Matrix &quad, Vector &lin, Real &c) const
 239 {
 240   quad.fill (0);
 241   lin.fill (0);
 242   c = 0;
 243
 244   for (Cons<Idealspacing> *p =ideal_p_list_; p; p = p->next_)
 245     {
 246       Idealspacing *i = p->car_;
 247       int l = i->cols_drul_[LEFT];
 248       int r = i->cols_drul_[RIGHT];
 249
 250       quad (r,r) += i->hooke_f_;
 251       quad (r,l) -= i->hooke_f_;
 252       quad (l,r) -= i->hooke_f_;
 253       quad (l,l) += i->hooke_f_;
 254
 255       lin (r) -= i->space_f_*i->hooke_f_;
 256       lin (l) += i->space_f_*i->hooke_f_;
 257
 258       c += sqr (i->space_f_);
 259     }
 260
 261   if (quad.dim() > 10)
 262     quad.set_band();
 263
 264
 265 }
 266
 267 void
 268 Spring_spacer::set_fixed_cols (Mixed_qp &qp) const
 269 {
 270   for (int j=0; j < cols_.size(); j++)
 271     if (cols_[j].fixed_b())
 272       qp.add_fixed_var (j,cols_[j].fixed_position());
 273 }
 274
 275 // put the constraints into the LP problem
 276 void
 277 Spring_spacer::make_constraints (Mixed_qp& lp) const
 278 {
 279   int dim=cols_.size();
 280
 281   for (int j=0; j < dim -1; j++)
 282     {
 283       Array<Spacer_rod> const&rod_arr (cols_[j].rods_[RIGHT]);
 284       for (int i = 0; i < rod_arr.size (); i++)
 285         {
 286           Vector c1(dim);
 287           c1(rod_arr[i].other_idx_)=1.0 ;
 288           c1(j)=-1.0 ;
 289
 290           lp.add_inequality_cons (c1, rod_arr[i].distance_f_);
 291         }
 292     }
 293 }
 294
 295
 296 Real
 297 Spring_spacer::calculate_energy_f (Vector solution) const
 298 {
 299   Real e = 0.0;
 300   for (Cons<Idealspacing>*p =ideal_p_list_; p; p = p->next_)
 301     {
 302       Idealspacing * i = p->car_;
 303       e += i->energy_f(solution(i->cols_drul_[RIGHT]) - solution(i->cols_drul_[LEFT]));
 304     }
 305
 306   return e;
 307 }
 308 void
 309 Spring_spacer::lower_bound_solution (Column_x_positions*positions) const
 310 {
 311   Mixed_qp lp (cols_.size());
 312   make_matrices (lp.quad_,lp.lin_, lp.const_term_);
 313   set_fixed_cols (lp);
 314
 315   Vector start (cols_.size());
 316   start.fill (0.0);
 317   Vector solution_vec (lp.solve (start));
 318
 319   DOUT << "Lower bound sol: " << solution_vec;
 320   positions->energy_f_ = calculate_energy_f (solution_vec);
 321   positions->config = solution_vec;
 322   positions->satisfies_constraints_b_ = check_constraints (solution_vec);
 323 }
 324
 325 Spring_spacer::Spring_spacer ()
 326 {
 327   ideal_p_list_ =0;
 328   energy_normalisation_f_ = 1.0;
 329 }
 330
 331 void
 332 Spring_spacer::solve (Column_x_positions*positions) const
 333 {
 334   DOUT << "Spring_spacer::solve ()...";
 335
 336   Vector solution_try;
 337
 338   bool constraint_satisfaction = try_initial_solution_and_tell (solution_try);
 339   if  (constraint_satisfaction)
 340     {
 341       Mixed_qp lp (cols_.size());
 342       make_matrices (lp.quad_,lp.lin_, lp.const_term_);
 343       make_constraints (lp);
 344       set_fixed_cols (lp);
 345
 346       Vector solution_vec (lp.solve (solution_try));
 347
 348       positions->satisfies_constraints_b_ = check_constraints (solution_vec);
 349       if (!positions->satisfies_constraints_b_)
 350         {
 351           WARN << _ ("solution doesn't satisfy constraints") << '\n' ;
 352         }
 353       position_loose_cols (solution_vec);
 354       positions->energy_f_ = calculate_energy_f (solution_vec);
 355       positions->config = solution_vec;
 356       positions->error_col_l_arr_ = error_pcol_l_arr();
 357     }
 358   else
 359     {
 360       positions->set_stupid_solution (solution_try);
 361     }
 362
 363   DOUT << "Finished Spring_spacer::solve ()...";
 364 }
 365
 366 /**
 367   add one column to the problem.
 368 */
 369 void
 370 Spring_spacer::add_column (Paper_column  *col, bool fixed, Real fixpos)
 371 {
 372   Column_info c (col,(fixed)? &fixpos :  0);
 373   int this_rank =  cols_.size();
 374   c.rank_i_ = this_rank;
 375
 376   for (int i=0; i < col->minimal_dists_arr_drul_[LEFT].size (); i++)
 377     {
 378       Column_rod &cr = col->minimal_dists_arr_drul_[LEFT][i];
 379       int left_idx = cr.other_l_->rank_i () - cols_[0].pcol_l_->rank_i ();
 380       if (left_idx < 0)
 381         continue;
 382
 383       if (cols_[left_idx].pcol_l_ != cr.other_l_)
 384         continue;
 385
 386       Spacer_rod l_rod;
 387       l_rod.distance_f_ = cr.distance_f_;
 388       l_rod.other_idx_ = left_idx;
 389       c.rods_[LEFT].push (l_rod);
 390
 391       Spacer_rod r_rod;
 392       r_rod.distance_f_ = cr.distance_f_;
 393       r_rod.other_idx_ = this_rank;
 394       cols_[left_idx].rods_[RIGHT].push (r_rod);
 395     }
 396
 397   for (int i=0; i < col->spring_arr_drul_[LEFT].size (); i++)
 398     {
 399       Column_spring &cr = col->spring_arr_drul_[LEFT][i];
 400       int idx = cr.other_l_->rank_i () - cols_[0].pcol_l_->rank_i ();
 401       if (idx < 0)
 402         continue;
 403
 404       if (cols_[idx].pcol_l_ != cr.other_l_)
 405             continue;
 406
 407
 408       connect (idx, this_rank, cr.distance_f_,
 409                cr.strength_f_ / cr.distance_f_);
 410     }
 411
 412   cols_.push (c);
 413 }
 414
 415 Line_of_cols
 416 Spring_spacer::error_pcol_l_arr() const
 417 {
 418   Link_array<Paper_column> retval;
 419   for (int i=0; i< cols_.size(); i++)
 420     if (cols_[i].ugh_b_)
 421       retval.push (cols_[i].pcol_l_);
 422   for (int i=0;  i < loose_col_arr_.size(); i++)
 423     {
 424       retval.push (loose_col_arr_[i].pcol_l_);
 425     }
 426   return retval;
 427 }
 428
 429 /*
 430   Ugh. Should junk this.
 431  */
 432 void
 433 Spring_spacer::loosen_column (int idx)
 434 {
 435   Column_info c=cols_.get (idx);
 436
 437   Cons<Idealspacing> **pp = &ideal_p_list_;
 438
 439   while (*pp)
 440     {
 441       Idealspacing *j = (*pp)->car_;
 442       if (j->cols_drul_[LEFT] == idx|| j->cols_drul_[RIGHT] == idx)
 443         {
 444           delete remove_cons (pp);
 445         }
 446       else
 447         {
 448           pp = &(*pp)->next_;
 449         }
 450     }
 451   c.ugh_b_ = true;
 452
 453   int j=0;
 454   for (; j < loose_col_arr_.size(); j++)
 455     {
 456       if (loose_col_arr_[j].rank_i_ > c.rank_i_)
 457         break;
 458     }
 459   loose_col_arr_.insert (c,j);
 460 }
 461
 462
 463 void
 464 Spring_spacer::print() const
 465 {
 466 #ifndef NPRINT
 467   for (int i=0; i < cols_.size(); i++)
 468     {
 469       DOUT << "col " << i << " ";
 470       cols_[i].print();
 471     }
 472
 473   for (Cons<Idealspacing> *p =ideal_p_list_; p; p = p->next_)
 474     {
 475       p->car_->print();
 476     }
 477 #endif
 478 }
 479
 480
 481 void
 482 Spring_spacer::connect (int i, int j, Real d, Real h)
 483 {
 484   assert(d >= 0 && d <= 100 CM);
 485   assert(h >=0);
 486
 487   Idealspacing * s = new Idealspacing;
 488
 489   s->cols_drul_[LEFT] = i ;
 490   s->cols_drul_[RIGHT] = j;
 491   s->space_f_ = d;
 492   s->hooke_f_ = h;
 493
 494   ideal_p_list_ = new Killing_cons<Idealspacing> (s, ideal_p_list_);
 495 }
 496
 497 void
 498 Spring_spacer::prepare()
 499 {
 500   handle_loose_cols();
 501   print();
 502 }
 503
 504 Line_spacer*
 505 Spring_spacer::constructor()
 506 {
 507   return new Spring_spacer;
 508 }
 509
 510
 511