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_cons_p_)
  72     {
  73       connected.connect (i->car_p_->cols_drul_[LEFT],i->car_p_->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 distances and fixed positions
 182  */
 183 Vector
 184 Spring_spacer::try_initial_solution() const
 185 {
 186   Vector v;
 187   if (!try_initial_solution_and_tell (v))
 188     {
 189       warning (_ ("I'm too fat; call Oprah"));
 190     }
 191   return v;
 192
 193 }
 194
 195 bool
 196 Spring_spacer::try_initial_solution_and_tell (Vector &v) const
 197 {
 198   int dim=cols_.size();
 199   bool succeeded = true;
 200   Vector initsol (dim);
 201
 202   assert (cols_[0].fixed_b ());
 203   DOUT << "fixpos 0 " << cols_[0].fixed_position ();
 204   for (int i=0; i < dim; i++)
 205     {
 206       Real min_x = i ?  initsol (i-1) : cols_[0].fixed_position ();
 207       Array<Spacer_rod> const &sr_arr(cols_[i].rods_[LEFT]);
 208       for (int j=0; j < sr_arr.size (); j++)
 209         {
 210           min_x = min_x >? (initsol (sr_arr[j].other_idx_) + sr_arr[j].distance_f_);
 211         }
 212       initsol (i) = min_x;
 213
 214       if (cols_[i].fixed_b())
 215         {
 216           initsol (i)=cols_[i].fixed_position();
 217           if (initsol (i) < min_x )
 218             {
 219               DOUT << "failing: init, min : " << initsol (i) << " " << min_x << '\n';
 220               initsol (i) = min_x;
 221               succeeded = false;
 222             }
 223         }
 224     }
 225   v = initsol;
 226
 227   DOUT << "tried and told solution: " << v;
 228   if (!succeeded)
 229     DOUT << "(failed)\n";
 230   return succeeded;
 231 }
 232
 233
 234
 235 // generate the matrices
 236 void
 237 Spring_spacer::make_matrices (Matrix &quad, Vector &lin, Real &c) const
 238 {
 239   quad.fill (0);
 240   lin.fill (0);
 241   c = 0;
 242
 243   for (Cons<Idealspacing> *p =ideal_p_list_; p; p = p->next_cons_p_)
 244     {
 245       Idealspacing *i = p->car_p_;
 246       int l = i->cols_drul_[LEFT];
 247       int r = i->cols_drul_[RIGHT];
 248
 249       quad (r,r) += i->hooke_f_;
 250       quad (r,l) -= i->hooke_f_;
 251       quad (l,r) -= i->hooke_f_;
 252       quad (l,l) += i->hooke_f_;
 253
 254       lin (r) -= i->space_f_*i->hooke_f_;
 255       lin (l) += i->space_f_*i->hooke_f_;
 256
 257       c += sqr (i->space_f_);
 258     }
 259
 260   if (quad.dim() > 10)
 261     quad.set_band();
 262
 263
 264 }
 265
 266 void
 267 Spring_spacer::set_fixed_cols (Mixed_qp &qp) const
 268 {
 269   for (int j=0; j < cols_.size(); j++)
 270     if (cols_[j].fixed_b())
 271       qp.add_fixed_var (j,cols_[j].fixed_position());
 272 }
 273
 274 // put the constraints into the LP problem
 275 void
 276 Spring_spacer::make_constraints (Mixed_qp& lp) const
 277 {
 278   int dim=cols_.size();
 279
 280   for (int j=0; j < dim -1; j++)
 281     {
 282       Array<Spacer_rod> const&rod_arr (cols_[j].rods_[RIGHT]);
 283       for (int i = 0; i < rod_arr.size (); i++)
 284         {
 285           Vector c1(dim);
 286           c1(rod_arr[i].other_idx_)=1.0 ;
 287           c1(j)=-1.0 ;
 288
 289           lp.add_inequality_cons (c1, rod_arr[i].distance_f_);
 290         }
 291     }
 292 }
 293
 294
 295 Real
 296 Spring_spacer::calculate_energy_f (Vector solution) const
 297 {
 298   Real e = 0.0;
 299   for (Cons<Idealspacing>*p =ideal_p_list_; p; p = p->next_cons_p_)
 300     {
 301       Idealspacing * i = p->car_p_;
 302       e += i->energy_f(solution(i->cols_drul_[RIGHT]) - solution(i->cols_drul_[LEFT]));
 303     }
 304
 305   return e;
 306 }
 307 void
 308 Spring_spacer::lower_bound_solution (Column_x_positions*positions) const
 309 {
 310   Mixed_qp lp (cols_.size());
 311   make_matrices (lp.quad_,lp.lin_, lp.const_term_);
 312   set_fixed_cols (lp);
 313
 314   Vector start (cols_.size());
 315   start.fill (0.0);
 316   Vector solution_vec (lp.solve (start));
 317
 318   DOUT << "Lower bound sol: " << solution_vec;
 319   positions->energy_f_ = calculate_energy_f (solution_vec);
 320   positions->config = solution_vec;
 321   positions->satisfies_constraints_b_ = check_constraints (solution_vec);
 322 }
 323
 324 Spring_spacer::Spring_spacer ()
 325 {
 326   ideal_p_list_ =0;
 327   energy_normalisation_f_ = 1.0;
 328 }
 329
 330 void
 331 Spring_spacer::solve (Column_x_positions*positions) const
 332 {
 333   DOUT << "Spring_spacer::solve ()...";
 334
 335   Vector solution_try;
 336
 337   bool constraint_satisfaction = try_initial_solution_and_tell (solution_try);
 338   if  (constraint_satisfaction)
 339     {
 340       Mixed_qp lp (cols_.size());
 341       make_matrices (lp.quad_,lp.lin_, lp.const_term_);
 342       make_constraints (lp);
 343       set_fixed_cols (lp);
 344
 345       Vector solution_vec (lp.solve (solution_try));
 346
 347       positions->satisfies_constraints_b_ = check_constraints (solution_vec);
 348       if (!positions->satisfies_constraints_b_)
 349         {
 350           WARN << _ ("solution doesn't satisfy constraints") << '\n' ;
 351         }
 352       position_loose_cols (solution_vec);
 353       positions->energy_f_ = calculate_energy_f (solution_vec);
 354       positions->config = solution_vec;
 355       positions->error_col_l_arr_ = error_pcol_l_arr();
 356     }
 357   else
 358     {
 359       positions->set_stupid_solution (solution_try);
 360     }
 361
 362   DOUT << "Finished Spring_spacer::solve ()...";
 363 }
 364
 365 /**
 366   add one column to the problem.
 367 */
 368 void
 369 Spring_spacer::add_column (Paper_column  *col, bool fixed, Real fixpos)
 370 {
 371   Column_info c (col,(fixed)? &fixpos :  0);
 372   int this_rank =  cols_.size();
 373   c.rank_i_ = this_rank;
 374
 375   for (int i=0; i < col->minimal_dists_arr_drul_[LEFT].size (); i++)
 376     {
 377       Column_rod &cr = col->minimal_dists_arr_drul_[LEFT][i];
 378       int left_idx = cr.other_l_->rank_i () - cols_[0].pcol_l_->rank_i ();
 379       if (left_idx < 0)
 380         continue;
 381
 382       if (cols_[left_idx].pcol_l_ != cr.other_l_)
 383         continue;
 384
 385       Spacer_rod l_rod;
 386       l_rod.distance_f_ = cr.distance_f_;
 387       l_rod.other_idx_ = left_idx;
 388       c.rods_[LEFT].push (l_rod);
 389
 390       Spacer_rod r_rod;
 391       r_rod.distance_f_ = cr.distance_f_;
 392       r_rod.other_idx_ = this_rank;
 393       cols_[left_idx].rods_[RIGHT].push (r_rod);
 394     }
 395
 396   cols_.push (c);
 397 }
 398
 399 Line_of_cols
 400 Spring_spacer::error_pcol_l_arr() const
 401 {
 402   Link_array<Paper_column> retval;
 403   for (int i=0; i< cols_.size(); i++)
 404     if (cols_[i].ugh_b_)
 405       retval.push (cols_[i].pcol_l_);
 406   for (int i=0;  i < loose_col_arr_.size(); i++)
 407     {
 408       retval.push (loose_col_arr_[i].pcol_l_);
 409     }
 410   return retval;
 411 }
 412 /*
 413   Ugh. Should junk this.
 414  */
 415 void
 416 Spring_spacer::loosen_column (int idx)
 417 {
 418   Column_info c=cols_.get (idx);
 419
 420   Cons<Idealspacing> **pp = &ideal_p_list_;
 421
 422   while (*pp)
 423     {
 424       Idealspacing *j = (*pp)->car_p_;
 425       if (j->cols_drul_[LEFT] == idx|| j->cols_drul_[RIGHT] == idx)
 426         {
 427           delete remove_cons_p (pp);
 428         }
 429       else
 430         {
 431           pp = &(*pp)->next_cons_p_;
 432         }
 433     }
 434   c.ugh_b_ = true;
 435
 436   int j=0;
 437   for (; j < loose_col_arr_.size(); j++)
 438     {
 439       if (loose_col_arr_[j].rank_i_ > c.rank_i_)
 440         break;
 441     }
 442   loose_col_arr_.insert (c,j);
 443 }
 444
 445
 446 void
 447 Spring_spacer::print() const
 448 {
 449 #ifndef NPRINT
 450   for (int i=0; i < cols_.size(); i++)
 451     {
 452       DOUT << "col " << i << " ";
 453       cols_[i].print();
 454     }
 455
 456   for (Cons<Idealspacing> *p =ideal_p_list_; p; p = p->next_cons_p_)
 457     {
 458       p->car_p_->print();
 459     }
 460 #endif
 461 }
 462
 463
 464 void
 465 Spring_spacer::connect (int i, int j, Real d, Real h)
 466 {
 467   assert(d >= 0 && d <= 100 CM);
 468   assert(h >=0);
 469
 470   Idealspacing * s = new Idealspacing;
 471
 472   s->cols_drul_[LEFT] = i ;
 473   s->cols_drul_[RIGHT] = j;
 474   s->space_f_ = d;
 475   s->hooke_f_ = h;
 476
 477   ideal_p_list_ = new Killing_cons<Idealspacing> (s, ideal_p_list_);
 478 }
 479
 480
 481
 482
 483 void
 484 Spring_spacer::prepare()
 485 {
 486   DOUT << "Preparing..";
 487   calc_idealspacing();
 488   handle_loose_cols();
 489   print();
 490   DOUT << "finished preparing.\n";
 491 }
 492
 493 Line_spacer*
 494 Spring_spacer::constructor()
 495 {
 496   return new Spring_spacer;
 497 }
 498
 499
 500
 501 /**
 502   get the shortest_playing running note at a time. */
 503 void
 504 Spring_spacer::get_ruling_durations(Array<Moment> &shortest_playing_arr,
 505                                     Array<Moment> &context_shortest_arr)
 506 {
 507   for (int i=0; i < cols_.size(); i++)
 508     {
 509       scol_l (i)->preprocess();
 510       scol_l (i)->print ();
 511     }
 512   int start_context_i=0;
 513   Moment context_shortest;
 514   context_shortest.set_infinite (1);
 515   context_shortest_arr.set_size(cols_.size());
 516
 517   for (int i=0; i < cols_.size(); i++)
 518     {
 519       Score_column * sc = scol_l(i);
 520       Moment now = scol_l (i)->when();
 521       Moment shortest_playing;
 522       shortest_playing.set_infinite (1);
 523
 524       if (!sc->musical_b ())
 525         {
 526           for (int ji=i; ji >= start_context_i; ji--)
 527             context_shortest_arr[ji] = context_shortest;
 528           start_context_i = i;
 529           context_shortest.set_infinite (1);
 530         }
 531       if (sc->durations.size())
 532         {
 533           context_shortest = context_shortest <? sc->durations[0];
 534         }
 535
 536       // ji was j, but triggered ICE
 537       for (int ji=i+1; ji --;)
 538         {
 539           if (scol_l(ji)->durations.size() &&
 540               now - scol_l(ji)->when() >= shortest_playing)
 541             break;
 542
 543           for (int k =  scol_l (ji)->durations.size();
 544                k-- && scol_l(ji)->durations[k] + scol_l(ji)->when() > now;
 545                )
 546             {
 547               shortest_playing = shortest_playing <? scol_l(ji)->durations[k];
 548             }
 549         }
 550       shortest_playing_arr.push(shortest_playing);
 551     }
 552
 553 #ifndef NPRINT
 554   DOUT << "shortest_playing/:[ ";
 555   for (int i=0; i < shortest_playing_arr.size(); i++)
 556     {
 557       DOUT << shortest_playing_arr[i] << " ";
 558       DOUT << context_shortest_arr[i] << ", ";
 559     }
 560   DOUT << "]\n";
 561 #endif
 562 }
 563
 564 /*
 565   TODO: take out the refs to width
 566
 567  */
 568 /**
 569   generate springs between columns.
 570
 571   TODO: This needs rethinking.......
 572
 573   *  Spacing should take optical
 574   effects into account
 575
 576   *  Should be decentralised
 577
 578   The algorithm is taken from :
 579
 580   John S. Gourlay. ``Spacing a Line of Music,'' Technical Report
 581   OSU-CISRC-10/87-TR35, Department of Computer and Information
 582   Science, The Ohio State University, 1987.
 583
 584   */
 585 void
 586 Spring_spacer::calc_idealspacing()
 587 {
 588   Array<Moment> shortest_playing_arr;
 589   Array<Moment> context_shortest_arr;
 590   get_ruling_durations(shortest_playing_arr, context_shortest_arr);
 591
 592   Real interline_f = paper_l ()->get_realvar (interline_scm_sym);
 593
 594   Array<Real> ideal_arr;
 595   Array<Real> hooke_arr;
 596   for (int i=0; i < cols_.size() - 1; i++){
 597     ideal_arr.push (-1.0);
 598     hooke_arr.push (1.0);
 599   }
 600
 601   /*
 602      First do all non-musical columns
 603   */
 604   for (int i=0; i < cols_.size(); i++)
 605     {
 606       if (!scol_l (i)->musical_b() && i+1 < cols_.size())
 607         {
 608           Real symbol_distance =cols_[i].width_[RIGHT] + 2 PT;
 609           Real durational_distance = 0;
 610           Moment delta_t =  scol_l (i+1)->when() - scol_l (i)->when () ;
 611
 612           /*
 613             ugh should use shortest_playing distance
 614           */
 615           if (delta_t)
 616             {
 617               Real k=  paper_l()->arithmetic_constant (context_shortest_arr[i]);
 618               durational_distance =  paper_l()->length_mom_to_dist (delta_t,k);
 619             }
 620           symbol_distance += -cols_[i+1].width_[LEFT];
 621
 622
 623           ideal_arr[i] = symbol_distance >? durational_distance;
 624           hooke_arr[i] = 1; //2.0;
 625         }
 626     }
 627
 628   /*
 629      Then musicals
 630   */
 631   for (int i=0; i < cols_.size(); i++)
 632     {
 633       if (scol_l (i)->musical_b())
 634         {
 635           Moment shortest_playing_len = shortest_playing_arr[i];
 636           Moment context_shortest = context_shortest_arr[i];
 637           if (! shortest_playing_len)
 638             {
 639               warning (_f ("can't find a ruling note at %s",
 640                 scol_l (i)->when().str ()));
 641               shortest_playing_len = 1;
 642             }
 643           if (! context_shortest)
 644             {
 645               warning (_f ("no minimum in measure at %s",
 646                       scol_l (i)->when().str ()));
 647               context_shortest = 1;
 648             }
 649           Moment delta_t = scol_l (i+1)->when() - scol_l (i)->when ();
 650           Real k=  paper_l()->arithmetic_constant(context_shortest);
 651           Real dist = paper_l()->length_mom_to_dist (shortest_playing_len, k);
 652           dist *= (double)(delta_t / shortest_playing_len);
 653
 654           /*
 655             According to [Ross] and [Wanske], and from what i've seen:
 656
 657             * whitespace at the begin of the bar should be fixed at
 658             (about) one interline.
 659
 660             [Ross]:
 661             when spacing gets real tight, a smaller fixed value may be
 662             used, so that there are two discrete amounts of whitespace
 663             possible at the begin of a bar; but this is not implemented
 664             right now.
 665
 666             * whitespace at the end of the bar is the normal amount of
 667             "hinterfleish" that would have been used, had there been
 668             yet another note in the bar.
 669
 670             [Ross]:
 671             some editors argue that the bar line should not take any
 672             space, not to hinder the flow of music spaced around a bar
 673             line.
 674
 675             [Ross] and [Wanske] do not suggest this, however.  Further,
 676             it introduces some spacing problems and I think that it is ugly
 677             too.
 678
 679             [jcn]
 680           */
 681
 682           /*
 683              first musical column of bar
 684           */
 685           if (i && !scol_l (i - 1)->musical_b ())
 686             {
 687               // one interline minimum at start of bar
 688
 689               // cols_[i].width_[RIGHT] += interline_f;
 690               cols_[i].width_[RIGHT] = cols_[i].width_[RIGHT] >? interline_f;
 691
 692               // should adjust dist too?
 693               ideal_arr[i-1] = ideal_arr[i-1] >? (2 * interline_f);
 694             }
 695
 696           /*
 697              last musical column of bar
 698           */
 699           if (i + 1 < cols_.size () && !scol_l(i+1)->musical_b ())
 700             {
 701               // two interline minimum ok for last column?
 702               dist = dist >? 2 * interline_f;
 703
 704               // set minimum rod
 705               /*
 706                 urg: simply *adding* an interline leaves big gaps at
 707                 end of measure in star-spangled-banner (after lyrics
 708                 at eom).
 709
 710                    cols_[i].width_[RIGHT] += interline_f; // before
 711
 712                 having a minimum of one interline solves this problem
 713                 in most (but not all??) cases.
 714
 715                 for music without lyrics (esp. when set very tightly),
 716                 adding an interline looks good: probably because this
 717                 hides a bug that allows the last note's "hinterfleish"
 718                 to be removed (e.g., see wtk1-fugue2: that's ugly now).
 719                 -- jcn
 720                */
 721
 722               cols_[i].width_[RIGHT] = cols_[i].width_[RIGHT] >? 2 * interline_f;
 723             }
 724
 725           // ugh, do we need this?
 726           if (i < cols_.size () - 1 && !scol_l (i + 1)->musical_b ())
 727             {
 728               Real minimum = -cols_[i + 1].width_[LEFT] + cols_[i].width_[RIGHT]
 729                 + interline_f / 2;
 730               dist = dist >? minimum;
 731             }
 732           ideal_arr[i] = dist;
 733         }
 734     }
 735
 736   /*
 737     shorter distances should stretch less.
 738
 739     (and how bout
 740
 741       hooke[i] = 2 * max_ideal_space - ideal[i]
 742
 743     ?)
 744   */
 745   for (int i=0; i < ideal_arr.size(); i++)
 746     hooke_arr[i] = 1/ideal_arr[i];
 747
 748   for (int i=0; i < ideal_arr.size(); i++)
 749     {
 750       assert (ideal_arr[i] >=0 && hooke_arr[i] >=0);
 751       connect (i, i+1, ideal_arr[i], hooke_arr[i]);
 752     }
 753 }