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