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
   9 const Real COLFUDGE=1e-3;
  10
  11
  12 bool
  13 Spacing_problem::contains(PCol const *w)
  14 {
  15     for (int i=0; i< cols.size(); i++)
  16         if (cols[i].pcol_ == w)
  17             return true;
  18     return false;
  19 }
  20
  21 int
  22 Spacing_problem::col_id(PCol const *w)const
  23 {
  24     for (int i=0; i< cols.size(); i++)
  25         if (cols[i].pcol_ == w)
  26             return i;
  27     assert(false);
  28     return -1;
  29 }
  30
  31 void
  32 Spacing_problem::OK() const
  33 {
  34 #ifndef NDEBUG
  35     Union_find connected(cols.size());
  36     Array<int> fixed;
  37     for (int i=0; i < ideals.size(); i++) {
  38         assert(ideals[i]->hooke > 0);
  39         int l = col_id(ideals[i]->left);
  40         int r = col_id(ideals[i]->right);
  41         connected.connect(l,r);
  42     }
  43     for (int i = 0; i < cols.size(); i++)
  44         if (cols[i].fixed())
  45             fixed.push(i);
  46     for (int i = 0; i < cols.size(); i++) {
  47         bool c=false;
  48         for (int j =0; j<fixed.size(); j++)
  49             c |=  connected.equiv(fixed[j],i);
  50         if (!c)
  51             WARN << "You have unconnected columns. \n"
  52                 "Check if bars and music fit each other\n"
  53                 "(crashing :-)\n";
  54         assert(c);
  55     }
  56 #endif
  57 }
  58
  59 bool
  60 Spacing_problem::check_constraints(Vector v) const
  61 {
  62     int dim=v.dim();
  63     for (int i=0; i < dim; i++) {
  64
  65         if (cols[i].fixed()&&
  66             abs(cols[i].fixed_position() - v(i)) > COLFUDGE)
  67             return false;
  68
  69         if (!i)
  70             continue;
  71
  72         Real mindist=cols[i-1].minright()
  73             +cols[i].minleft();
  74
  75         // ugh... compares
  76         Real dif =v(i) - v(i-1)- mindist;
  77         bool b = (dif > - COLFUDGE);
  78
  79
  80         if (!b)
  81             return false;
  82
  83     }
  84     return true;
  85 }
  86
  87 bool
  88 Spacing_problem::check_feasible() const
  89 {
  90     Vector sol(try_initial_solution());
  91     return check_constraints(sol);
  92 }
  93
  94 // generate a solution which obeys the min distances and fixed positions
  95 Vector
  96 Spacing_problem::try_initial_solution() const
  97 {
  98     int dim=cols.size();
  99     Vector initsol(dim);
 100     for (int i=0; i < dim; i++) {
 101         if (cols[i].fixed()) {
 102             initsol(i)=cols[i].fixed_position();
 103         } else {
 104             Real mindist=cols[i-1].minright()
 105                 +cols[i].minleft();
 106             assert(mindist >= 0.0);
 107             initsol(i)=initsol(i-1)+mindist;
 108
 109             //nog niet
 110             //if (i>0)
 111             //  assert(initsol(i) > initsol(i-1));
 112         }
 113     }
 114
 115     return initsol;
 116 }
 117 Vector
 118 Spacing_problem::find_initial_solution() const
 119 {
 120     Vector v(try_initial_solution());
 121     assert(check_constraints(v));
 122     return v;
 123 }
 124
 125 // generate the matrices
 126 void
 127 Spacing_problem::make_matrices(Matrix &quad, Vector &lin, Real &c) const
 128 {
 129     quad.fill(0);
 130     lin.fill(0);
 131     c = 0;
 132     for (int j=0; j < ideals.size(); j++){
 133         Idealspacing const*i=ideals[j];
 134         int l = col_id(i->left);
 135         int r = col_id(i->right);
 136
 137         quad(r,r) += i->hooke;
 138         quad(r,l) -= i->hooke;
 139         quad(l,r) -= i->hooke;
 140         quad(l,l) += i->hooke;
 141
 142         lin(r) -= i->space*i->hooke;
 143         lin(l) += i->space*i->hooke;
 144
 145         c += sqr(i->space);
 146     }
 147 }
 148
 149 // put the constraints into the LP problem
 150 void
 151 Spacing_problem::make_constraints(Mixed_qp& lp) const
 152 {
 153     int dim=cols.size();
 154     for (int j=0; j < dim; j++) {
 155         Colinfo *c=&(cols[j]);
 156         if (c->fixed()) {
 157             lp.add_fixed_var(j,c->fixed_position());
 158         }
 159         if (j > 0){
 160             Vector c1(dim);
 161
 162             c1(j)=1.0 ;
 163             c1(j-1)=-1.0 ;
 164             lp.add_inequality_cons(c1, cols[j-1].minright() +
 165                                    cols[j].minleft());
 166         }
 167     }
 168 }
 169
 170 Array<Real>
 171 Spacing_problem::solve() const
 172 {
 173     print();
 174     OK();
 175     assert(check_feasible());
 176
 177     /* optimalisatiefunctie */
 178     Mixed_qp lp(cols.size());
 179     make_matrices(lp.quad,lp.lin, lp.const_term);
 180     make_constraints(lp);
 181     Vector start=find_initial_solution();
 182     Vector sol(lp.solve(start));
 183     if (!check_constraints(sol)) {
 184         WARN << "solution doesn't satisfy constraints.\n" ;
 185     }
 186
 187
 188     Array<Real> posns(sol);
 189     posns.push(lp.eval(sol));
 190     return posns;
 191 }
 192
 193 /**
 194     add one column to the problem.
 195 */
 196 void
 197 Spacing_problem::add_column(PCol const *col, bool fixed, Real fixpos)
 198 {
 199     if (!col->used_b() )
 200         return;
 201     Colinfo c(col,(fixed)? &fixpos :  0);
 202     cols.push(c);
 203 }
 204
 205 void
 206 Spacing_problem::add_ideal(Idealspacing const *i)
 207 {
 208     PCol const *l =i->left;
 209     PCol const *r= i->right;
 210
 211     if (!contains(l) || !contains(r)) {
 212         return;
 213     }
 214     ideals.push(i);
 215 }
 216
 217 void
 218 Spacing_problem::print_ideal(Idealspacing const *id)const
 219 {
 220 #ifndef NPRINT
 221     int l = col_id(id->left);
 222     int r = col_id(id->right);
 223
 224     mtor << "between " << l <<","<<r<<":" ;
 225     id->print();
 226 #endif
 227 }
 228
 229 void
 230 Spacing_problem::print() const
 231 {
 232 #ifndef NPRINT
 233     for (int i=0; i < cols.size(); i++) {
 234         mtor << "col " << i<<' ';
 235         cols[i].print();
 236     }
 237     for (int i=0; i < ideals.size(); i++) {
 238         print_ideal(ideals[i]);
 239     }
 240 #endif
 241
 242 }
 243
 244 /* **************** */
 245
 246 void
 247 Colinfo::print() const
 248 {
 249 #ifndef NPRINT
 250     mtor << "column { ";
 251     if (fixed())
 252         mtor << "fixed at " << fixed_position()<<", ";
 253     assert(pcol_);
 254     mtor << "[" << minleft() << ", " << minright() << "]";
 255     mtor <<"}\n";
 256 #endif
 257 }
 258
 259 Colinfo::Colinfo(Colinfo const&c)
 260 {
 261     fixpos = (c.fixpos)?new Real(*c.fixpos):0;
 262     pcol_ = c.pcol_;
 263     width = c.width;
 264 }
 265
 266 Colinfo::Colinfo(PCol const *col_p, Real const *fixed_r_p )
 267 {
 268     fixpos = (fixed_r_p)? new Real(*fixed_r_p) : 0;
 269     pcol_ = col_p;
 270     width = pcol_->width();
 271 }
 272
 273 Colinfo::~Colinfo()
 274 {
 275     delete fixpos;
 276 }
 277
 278 Colinfo::Colinfo()
 279 {
 280     pcol_=0;
 281     fixpos = 0;
 282 }
 283 void
 284 Colinfo::operator=(Colinfo const&c )
 285 {
 286     delete fixpos;
 287     fixpos = (c.fixpos)?new Real(*c.fixpos):0;
 288     pcol_ = c.pcol_;
 289     width = c.width;
 290 }