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