src/linespace.cc

   1 #include <math.h>
   2 #include "linespace.hh"
   3 #include "pcol.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(const PCol *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(const PCol *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(const PCol *col, bool fixed, Real fixpos)
 198 {
 199     Colinfo c(col,(fixed)? &fixpos :  0);
 200     cols.push(c);
 201 }
 202
 203 void
 204 Spacing_problem::add_ideal(const Idealspacing *i)
 205 {
 206     const PCol *l =i->left;
 207     const PCol *r= i->right;
 208
 209     if (!contains(l) || !contains(r)) {
 210         return;
 211     }
 212     ideals.push(i);
 213 }
 214
 215 void
 216 Spacing_problem::print_ideal(const Idealspacing*id)const
 217 {
 218 #ifndef NPRINT
 219     int l = col_id(id->left);
 220     int r = col_id(id->right);
 221
 222     mtor << "between " << l <<","<<r<<":" ;
 223     id->print();
 224 #endif
 225 }
 226
 227 void
 228 Spacing_problem::print() const
 229 {
 230 #ifndef NPRINT
 231     for (int i=0; i < cols.size(); i++) {
 232         mtor << "col " << i<<' ';
 233         cols[i].print();
 234     }
 235     for (int i=0; i < ideals.size(); i++) {
 236         print_ideal(ideals[i]);
 237     }
 238 #endif
 239
 240 }
 241
 242 /* **************** */
 243
 244 void
 245 Colinfo::print() const
 246 {
 247 #ifndef NPRINT
 248     mtor << "column { ";
 249     if (fixed())
 250         mtor << "fixed at " << fixed_position()<<", ";
 251     assert(pcol_);
 252     mtor << "[" << minleft() << ", " << minright() << "]";
 253     mtor <<"}\n";
 254 #endif
 255 }
 256
 257 Colinfo::Colinfo(Colinfo const&c)
 258 {
 259     fixpos = (c.fixpos)?new Real(*c.fixpos):0;
 260     pcol_ = c.pcol_;
 261     width = c.width;
 262 }
 263
 264 Colinfo::Colinfo(const PCol*col_p, const Real*fixed_r_p )
 265 {
 266     fixpos = (fixed_r_p)? new Real(*fixed_r_p) : 0;
 267     pcol_ = col_p;
 268     width = pcol_->width();
 269 }
 270
 271 Colinfo::~Colinfo()
 272 {
 273     delete fixpos;
 274 }
 275
 276 Colinfo::Colinfo()
 277 {
 278     pcol_=0;
 279     fixpos = 0;
 280 }
 281 void
 282 Colinfo::operator=(Colinfo const&c )
 283 {
 284     delete fixpos;
 285     fixpos = (c.fixpos)?new Real(*c.fixpos):0;
 286     pcol_ = c.pcol_;
 287     width = c.width;
 288 }