lily/linespace.cc

   1 #include <math.h>
   2 #include "linespace.hh"
   3 #include "p-col.hh"
   4 #include "debug.hh"
   5 #include "qlp.hh"
   6 #include "unionfind.hh"
   7 #include "idealspacing.hh"
   8 #include "pointer.tcc"
   9
  10 const Real COLFUDGE=1e-3;
  11 template class P<Real>;         // ugh.
  12
  13 bool
  14 Spacing_problem::contains(PCol const *w)
  15 {
  16     for (int i=0; i< cols.size(); i++)
  17         if (cols[i].pcol_l_ == w)
  18             return true;
  19     return false;
  20 }
  21
  22 int
  23 Spacing_problem::col_id(PCol const *w)const
  24 {
  25     for (int i=0; i< cols.size(); i++)
  26         if (cols[i].pcol_l_ == w)
  27             return i;
  28     assert(false);
  29     return -1;
  30 }
  31
  32 void
  33 Spacing_problem::OK() const
  34 {
  35 #ifndef NDEBUG
  36     for (int i = 1; i < cols.size(); i++)
  37         assert(cols[i].rank_i_ > cols[i-1].rank_i_);
  38     for (int i = 1; i < loose_col_arr_.size(); i++)
  39         assert(loose_col_arr_[i].rank_i_ > loose_col_arr_[i-1].rank_i_);
  40 #endif
  41 }
  42
  43 /**
  44   Make sure no unconnected columns happen.
  45  */
  46 void
  47 Spacing_problem::handle_loose_cols()
  48 {
  49     Union_find connected(cols.size());
  50     Array<int> fixed;
  51     for (int i=0; i < ideals.size(); i++) {
  52         assert(ideals[i]->hooke > 0);
  53         int l = col_id(ideals[i]->left);
  54         int r = col_id(ideals[i]->right);
  55         connected.connect(l,r);
  56     }
  57     for (int i = 0; i < cols.size(); i++)
  58         if (cols[i].fixed())
  59             fixed.push(i);
  60     for (int i=1; i < fixed.size(); i++)
  61         connected.connect(fixed[i-1], fixed[i]);
  62
  63     for (int i = cols.size(); i--; ) {
  64         if (! connected.equiv(fixed[0], i)) {
  65             warning("unconnected column: " + String(i));
  66             loosen_column(i);
  67         }
  68     }
  69     OK();
  70 }
  71
  72
  73 /** Guess a stupid position for loose columns.  Put loose columns at
  74   regular distances from enclosing calced columns */
  75 void
  76 Spacing_problem::position_loose_cols(Vector &sol_vec)const
  77 {
  78     if (!loose_col_arr_.size())
  79         return ;
  80     assert(sol_vec.dim());
  81     Array<bool> fix_b_arr;
  82     fix_b_arr.set_size(cols.size() + loose_col_arr_.size());
  83     Real utter_right_f=-INFTY;
  84     Real utter_left_f =INFTY;
  85     for (int i=0; i < loose_col_arr_.size(); i++) {
  86         fix_b_arr[loose_col_arr_[i].rank_i_] = false;
  87     }
  88     for (int i=0; i < cols.size(); i++) {
  89         int r= cols[i].rank_i_;
  90         fix_b_arr[r] = true;
  91         utter_right_f = utter_right_f >? sol_vec(i);
  92         utter_left_f = utter_left_f <? sol_vec(i);
  93     }
  94     Vector v(fix_b_arr.size());
  95     int j =0;
  96     int k =0;
  97     for (int i=0; i < v.dim(); i++) {
  98         if (fix_b_arr[i]) {
  99             assert(cols[j].rank_i_ == i);
 100             v(i) = sol_vec(j++);
 101         } else {
 102             Real left_pos_f =
 103                 (j>0) ?sol_vec(j-1) : utter_left_f;
 104             Real right_pos_f =
 105                 (j < sol_vec.dim()) ? sol_vec(j) : utter_right_f;
 106             int left_rank = (j>0) ? cols[j-1].rank_i_ : 0;
 107             int right_rank = (j<sol_vec.dim()) ? cols[j].rank_i_ : sol_vec.dim();
 108
 109             int d_r = right_rank - left_rank;
 110             Colinfo loose=loose_col_arr_[k++];
 111             int r = loose.rank_i_ ;
 112             assert(r > left_rank && r < right_rank);
 113
 114             v(i) =  (r - left_rank)*left_pos_f/ d_r +
 115                 (right_rank - r) *right_pos_f /d_r;
 116         }
 117     }
 118     sol_vec = v;
 119 }
 120
 121 bool
 122 Spacing_problem::check_constraints(Vector v) const
 123 {
 124     int dim=v.dim();
 125     assert(dim == cols.size());
 126
 127     for (int i=0; i < dim; i++) {
 128
 129         if (cols[i].fixed()&&
 130             abs(cols[i].fixed_position() - v(i)) > COLFUDGE)
 131             return false;
 132
 133         if (!i)
 134             continue;
 135
 136         Real mindist=cols[i-1].minright()
 137             +cols[i].minleft();
 138
 139         // ugh... compares
 140         Real dif =v(i) - v(i-1)- mindist;
 141         bool b = (dif > - COLFUDGE);
 142
 143
 144         if (!b)
 145             return false;
 146
 147     }
 148     return true;
 149 }
 150
 151 void
 152 Spacing_problem::prepare()
 153 {
 154     handle_loose_cols();
 155 }
 156
 157 bool
 158 Spacing_problem::check_feasible() const
 159 {
 160     Vector sol(try_initial_solution());
 161     return check_constraints(sol);
 162 }
 163
 164 /// generate a solution which obeys the min distances and fixed positions
 165 Vector
 166 Spacing_problem::try_initial_solution() const
 167 {
 168     int dim=cols.size();
 169     Vector initsol(dim);
 170     for (int i=0; i < dim; i++) {
 171         if (cols[i].fixed()) {
 172             initsol(i)=cols[i].fixed_position();
 173
 174             if (i > 0) {
 175                 Real r =initsol(i-1)  + cols[i-1].minright();
 176                 if (initsol(i) < r ) {
 177                     warning("overriding fixed position");
 178                     initsol(i) =r;
 179                 }
 180             }
 181
 182         } else {
 183             Real mindist=cols[i-1].minright()
 184                 +cols[i].minleft();
 185             if (mindist < 0.0)
 186                 warning("Excentric column");
 187             initsol(i)=initsol(i-1)+mindist;
 188         }
 189     }
 190
 191     return initsol;
 192 }
 193
 194
 195
 196 Vector
 197 Spacing_problem::find_initial_solution() const
 198 {
 199     Vector v(try_initial_solution());
 200     assert(check_constraints(v));
 201     return v;
 202 }
 203
 204 // generate the matrices
 205 void
 206 Spacing_problem::make_matrices(Matrix &quad, Vector &lin, Real &c) const
 207 {
 208     quad.fill(0);
 209     lin.fill(0);
 210     c = 0;
 211     for (int j=0; j < ideals.size(); j++){
 212         Idealspacing const*i=ideals[j];
 213         int l = col_id(i->left);
 214         int r = col_id(i->right);
 215
 216         quad(r,r) += i->hooke;
 217         quad(r,l) -= i->hooke;
 218         quad(l,r) -= i->hooke;
 219         quad(l,l) += i->hooke;
 220
 221         lin(r) -= i->space*i->hooke;
 222         lin(l) += i->space*i->hooke;
 223
 224         c += sqr(i->space);
 225     }
 226 }
 227
 228 // put the constraints into the LP problem
 229 void
 230 Spacing_problem::make_constraints(Mixed_qp& lp) const
 231 {
 232     int dim=cols.size();
 233     for (int j=0; j < dim; j++) {
 234         Colinfo c=cols[j];
 235         if (c.fixed()) {
 236             lp.add_fixed_var(j,c.fixed_position());
 237         }
 238         if (j > 0){
 239             Vector c1(dim);
 240
 241             c1(j)=1.0 ;
 242             c1(j-1)=-1.0 ;
 243             lp.add_inequality_cons(c1, cols[j-1].minright() +
 244                                    cols[j].minleft());
 245         }
 246     }
 247 }
 248
 249 Array<Real>
 250 Spacing_problem::solve() const
 251 {
 252     assert(check_feasible());
 253     print();
 254
 255     Mixed_qp lp(cols.size());
 256     make_matrices(lp.quad,lp.lin, lp.const_term);
 257     make_constraints(lp);
 258     Vector start=find_initial_solution();
 259     Vector sol(lp.solve(start));
 260     if (!check_constraints(sol)) {
 261         WARN << "solution doesn't satisfy constraints.\n" ;
 262     }
 263     Real energy_f =lp.eval(sol);
 264     position_loose_cols(sol);
 265
 266     Array<Real> posns(sol);
 267
 268     posns.push(energy_f);
 269     return posns;
 270 }
 271
 272 /**
 273     add one column to the problem.
 274 */
 275 void
 276 Spacing_problem::add_column(PCol  *col, bool fixed, Real fixpos)
 277 {
 278     Colinfo c(col,(fixed)? &fixpos :  0);
 279     if (cols.size())
 280         c.rank_i_ = cols.top().rank_i_+1;
 281     else
 282         c.rank_i_ = 0;
 283     cols.push(c);
 284 }
 285
 286 Array<PCol*>
 287 Spacing_problem::error_pcol_l_arr()const
 288 {
 289     Array<PCol*> retval;
 290     for (int i=0; i< cols.size(); i++)
 291         if (cols[i].ugh_b_)
 292             retval.push(cols[i].pcol_l_);
 293     for (int i=0;  i < loose_col_arr_.size(); i++) {
 294         retval.push(loose_col_arr_[i].pcol_l_);
 295     }
 296     return retval;
 297 }
 298
 299 void
 300 Spacing_problem::loosen_column(int i)
 301 {
 302     Colinfo c=cols.get(i);
 303     for (int i=0; i < ideals.size(); ) {
 304         Idealspacing const *i_l =ideals[i];
 305         if (i_l->left == c.pcol_l_ || i_l->right == c.pcol_l_)
 306             ideals.del(i);
 307         else
 308             i++;
 309     }
 310     c.ugh_b_ = true;
 311
 312     int i=0;
 313     for (; i < loose_col_arr_.size(); i++) {
 314         if (loose_col_arr_[i].rank_i_ > c.rank_i_)
 315             break;
 316     }
 317     loose_col_arr_.insert(c,i);
 318 }
 319
 320 void
 321 Spacing_problem::add_ideal(Idealspacing const *i)
 322 {
 323     PCol const *l =i->left;
 324     PCol const *r= i->right;
 325
 326     if (!contains(l) || !contains(r)) {
 327         return;
 328     }
 329     ideals.push(i);
 330 }
 331
 332 void
 333 Spacing_problem::print_ideal(Idealspacing const *id)const
 334 {
 335 #ifndef NPRINT
 336     int l = col_id(id->left);
 337     int r = col_id(id->right);
 338
 339     mtor << "between " << l <<","<<r<<":" ;
 340     id->print();
 341 #endif
 342 }
 343
 344 void
 345 Spacing_problem::print() const
 346 {
 347 #ifndef NPRINT
 348     for (int i=0; i < cols.size(); i++) {
 349         mtor << "col " << i<<' ';
 350         cols[i].print();
 351     }
 352     for (int i=0; i < ideals.size(); i++) {
 353         print_ideal(ideals[i]);
 354     }
 355 #endif
 356
 357 }
 358
 359 /* **************** */
 360
 361 void
 362 Colinfo::print() const
 363 {
 364 #ifndef NPRINT
 365     mtor << "column { ";
 366     if (fixed())
 367         mtor << "fixed at " << fixed_position()<<", ";
 368     assert(pcol_l_);
 369     mtor << "[" << minleft() << ", " << minright() << "]";
 370     mtor <<"}\n";
 371 #endif
 372 }
 373
 374 Colinfo::Colinfo(PCol *col_l, Real const *fixed_C)
 375 {
 376     if (fixed_C)
 377         fixpos_p_.set_l(fixed_C);
 378     ugh_b_ = false;
 379     pcol_l_ = col_l;
 380     width = pcol_l_->width();
 381 }
 382
 383
 384 Colinfo::Colinfo()
 385 {
 386     ugh_b_ = false;
 387     pcol_l_ =0;
 388 }