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, Real &c) 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         c += sqr(i->space);
 139     }
 140 }
 141
 142 // put the constraints into the LP problem
 143 void
 144 Spacing_problem::make_constraints(Optimisation_problem& lp) const
 145 {
 146     int dim=cols.sz();
 147     for (int j=0; j < dim; j++) {
 148         Colinfo *c=&(cols[j]);
 149         if (c->fixed) {
 150             lp.add_fixed_var(j,c->fixpos);
 151         }
 152         if (j > 0){
 153             Vector c1(dim);
 154
 155
 156             c1(j)=1.0 ;
 157             c1(j-1)=-1.0 ;
 158             lp.add_inequality_cons(c1, cols[j-1].minright() +
 159                                    cols[j].minleft());
 160         }
 161     }
 162 }
 163
 164 svec<Real>
 165 Spacing_problem::solve() const
 166 {
 167     OK();
 168     assert(check_feasible());
 169     print();
 170
 171     /* optimalisatiefunctie */
 172     Optimisation_problem lp(cols.sz());
 173     make_matrices(lp.quad,lp.lin, lp.const_term);
 174     make_constraints(lp);
 175     Vector start=find_initial_solution();
 176     Vector sol(lp.solve(start));
 177     if (!check_constraints(sol)) {
 178         WARN << "solution doesn't satisfy constraints.\n" ;
 179     }
 180
 181
 182     svec<Real> posns(sol);
 183     posns.add(lp.eval(sol));
 184     return posns;
 185 }
 186
 187 /*
 188     add one column to the problem.
 189 */
 190 void
 191 Spacing_problem::add_column(const PCol *col, bool fixed, Real fixpos)
 192 {
 193     Colinfo c;
 194     c.fixed=fixed;
 195     c.fixpos=fixpos;
 196     c.col=col;
 197     cols.add(c);
 198 }
 199
 200 void
 201 Spacing_problem::add_ideal(const Idealspacing *i)
 202 {
 203     const PCol *l =i->left;
 204     const PCol *r= i->right;
 205
 206     if (!contains(l) || !contains(r)) {
 207         return;
 208     }
 209     ideals.add(i);
 210 }
 211
 212 void
 213 Spacing_problem::print_ideal(const Idealspacing*id)const
 214 {
 215     int l = col_id(id->left);
 216     int r = col_id(id->right);
 217
 218     mtor << "idealspacing { between " << l <<","<<r<<'\n';
 219     mtor << "distance "<<id->space<< " strength " << id->hooke << "}\n";
 220 }
 221
 222 void
 223 Spacing_problem::print() const
 224 {
 225     for (int i=0; i < cols.sz(); i++) {
 226         mtor << "col " << i<<' ';
 227         cols[i].print();
 228     }
 229     for (int i=0; i < ideals.sz(); i++) {
 230         print_ideal(ideals[i]);
 231     }
 232 }
 233
 234 void
 235 Colinfo::print() const
 236 {
 237     mtor << "column { ";
 238     if (fixed)
 239         mtor << "fixed at " << fixpos<<", ";
 240     mtor << "[" << minleft() << ", " << minright() << "]";
 241     mtor <<"}\n";
 242 }