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