src/linespace.cc

   1 #include <math.h>
   2 #include "linespace.hh"
   3 #include "debug.hh"
   4 #include "qlp.hh"
   5 #include "unionfind.hh"
   6
   7 const Real COLFUDGE=1e-3;
   8 //#define COLFUDGE 1e-3
   9 bool
  10 Spacing_problem::contains(const PCol *w)
  11 {
  12     for (int i=0; i< cols.sz(); i++)
  13         if (cols[i].pcol_ == w)
  14             return true;
  15     return false;
  16 }
  17
  18 int
  19 Spacing_problem::col_id(const PCol *w)const
  20 {
  21     for (int i=0; i< cols.sz(); i++)
  22         if (cols[i].pcol_ == w)
  23             return i;
  24     assert(false);
  25 }
  26
  27 void
  28 Spacing_problem::OK() const
  29 {
  30 #ifndef NDEBUG
  31     Union_find connected(cols.sz());
  32     svec<int> fixed;
  33     for (int i=0; i < ideals.sz(); i++) {
  34         assert(ideals[i]->hooke > 0);
  35         int l = col_id(ideals[i]->left);
  36         int r = col_id(ideals[i]->right);
  37         connected.connect(l,r);
  38     }
  39     for (int i = 0; i < cols.sz(); i++)
  40         if (cols[i].fixed)
  41             fixed.add(i);
  42     for (int i = 0; i < cols.sz(); i++) {
  43         bool c=false;
  44         for (int j =0; j<fixed.sz(); j++)
  45             c |=  connected.equiv(j,i);
  46         assert(c);
  47     }
  48 #endif
  49 }
  50
  51 bool
  52 Spacing_problem::check_constraints(Vector v) const
  53 {
  54     int dim=v.dim();
  55     // mtor << "checking solution " << v << '\n';
  56     for (int i=0; i < dim; i++) {
  57
  58         if (cols[i].fixed&& ABS(cols[i].fixpos - v(i)) > COLFUDGE) {
  59             return false;
  60         }
  61         if (!i)
  62             continue;
  63
  64         Real mindist=cols[i-1].minright()
  65             +cols[i].minleft();
  66
  67         // ugh... compares
  68         Real dif =v(i) - v(i-1)- mindist;
  69         bool b = (dif > - COLFUDGE);
  70
  71
  72 #if 1
  73         if (!b)
  74             return false;
  75
  76 #else
  77         mtor << "dif= "<<dif<<" fudge= " << COLFUDGE<< " dif >fudge= "<<
  78             b << "\n";
  79
  80         /* fucks up for unknown reasons */
  81         if (dif  < -COLFUDGE)
  82             return false;
  83 #endif
  84
  85     }
  86     return true;
  87 }
  88
  89 bool
  90 Spacing_problem::check_feasible() const
  91 {
  92     Vector sol(try_initial_solution());
  93     return check_constraints(sol);
  94 }
  95
  96 // generate a solution which obeys the min distances and fixed positions
  97 Vector
  98 Spacing_problem::try_initial_solution() const
  99 {
 100     int dim=cols.sz();
 101     Vector initsol(dim);
 102     for (int i=0; i < dim; i++) {
 103         if (cols[i].fixed) {
 104             initsol(i)=cols[i].fixpos;
 105         } else {
 106             Real mindist=cols[i-1].minright()
 107                 +cols[i].minleft();
 108             assert(mindist >= 0.0);
 109             initsol(i)=initsol(i-1)+mindist;
 110
 111             //nog niet
 112             //if (i>0)
 113             //  assert(initsol(i) > initsol(i-1));
 114         }
 115     }
 116
 117     return initsol;
 118 }
 119 Vector
 120 Spacing_problem::find_initial_solution() const
 121 {
 122     Vector v(try_initial_solution());
 123     assert(check_constraints(v));
 124     return v;
 125 }
 126 // generate the matrices
 127 void
 128 Spacing_problem::make_matrices(Matrix &quad, Vector &lin, Real &c) const
 129 {
 130     quad.fill(0);
 131     lin.fill(0);
 132     for (int j=0; j < ideals.sz(); 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.sz();
 154     for (int j=0; j < dim; j++) {
 155         Colinfo *c=&(cols[j]);
 156         if (c->fixed) {
 157             lp.add_fixed_var(j,c->fixpos);
 158         }
 159         if (j > 0){
 160             Vector c1(dim);
 161
 162
 163             c1(j)=1.0 ;
 164             c1(j-1)=-1.0 ;
 165             lp.add_inequality_cons(c1, cols[j-1].minright() +
 166                                    cols[j].minleft());
 167         }
 168     }
 169 }
 170
 171 svec<Real>
 172 Spacing_problem::solve() const
 173 {
 174     print();
 175     OK();
 176     assert(check_feasible());
 177
 178
 179     /* optimalisatiefunctie */
 180     Mixed_qp lp(cols.sz());
 181     make_matrices(lp.quad,lp.lin, lp.const_term);
 182     make_constraints(lp);
 183     Vector start=find_initial_solution();
 184     Vector sol(lp.solve(start));
 185     if (!check_constraints(sol)) {
 186         WARN << "solution doesn't satisfy constraints.\n" ;
 187     }
 188
 189
 190     svec<Real> posns(sol);
 191     posns.add(lp.eval(sol));
 192     return posns;
 193 }
 194
 195 /*
 196     add one column to the problem.
 197 */
 198 void
 199 Spacing_problem::add_column(const PCol *col, bool fixed, Real fixpos)
 200 {
 201     Colinfo c;
 202     c.fixed=fixed;
 203     c.fixpos=fixpos;
 204     assert(col);
 205     c.pcol_=col;
 206     cols.add(c);
 207 }
 208
 209 void
 210 Spacing_problem::add_ideal(const Idealspacing *i)
 211 {
 212     const PCol *l =i->left;
 213     const PCol *r= i->right;
 214
 215     if (!contains(l) || !contains(r)) {
 216         return;
 217     }
 218     ideals.add(i);
 219 }
 220
 221 void
 222 Spacing_problem::print_ideal(const Idealspacing*id)const
 223 {
 224 #ifndef NPRINT
 225     int l = col_id(id->left);
 226     int r = col_id(id->right);
 227
 228     mtor << "between " << l <<","<<r<<":" ;
 229     id->print();
 230 #endif
 231 }
 232
 233 void
 234 Spacing_problem::print() const
 235 {
 236 #ifndef NPRINT
 237     for (int i=0; i < cols.sz(); i++) {
 238         mtor << "col " << i<<' ';
 239         cols[i].print();
 240     }
 241     for (int i=0; i < ideals.sz(); i++) {
 242         print_ideal(ideals[i]);
 243     }
 244 #endif
 245
 246 }
 247
 248 void
 249 Colinfo::print() const
 250 {
 251 #ifndef NPRINT
 252     mtor << "column { ";
 253     if (fixed)
 254         mtor << "fixed at " << fixpos<<", ";
 255     assert(pcol_);
 256     mtor << "[" << minleft() << ", " << minright() << "]";
 257     mtor <<"}\n";
 258 #endif
 259 }
 260
 261 Colinfo::Colinfo()
 262 {
 263     fixed=false;
 264     pcol_=0;
 265 }