src/linespace.cc

   1 #include <math.h>
   2 #include "linespace.hh"
   3 #include "pcol.hh"
   4 #include "debug.hh"
   5 #include "qlp.hh"
   6 #include "unionfind.hh"
   7 #include "idealspacing.hh"
   8
   9 const Real COLFUDGE=1e-3;
  10
  11 //#define COLFUDGE 1e-3
  12 bool
  13 Spacing_problem::contains(const PCol *w)
  14 {
  15     for (int i=0; i< cols.sz(); i++)
  16         if (cols[i].pcol_ == w)
  17             return true;
  18     return false;
  19 }
  20
  21 int
  22 Spacing_problem::col_id(const PCol *w)const
  23 {
  24     for (int i=0; i< cols.sz(); i++)
  25         if (cols[i].pcol_ == w)
  26             return i;
  27     assert(false);
  28     return -1;
  29 }
  30
  31 void
  32 Spacing_problem::OK() const
  33 {
  34 #ifndef NDEBUG
  35     Union_find connected(cols.sz());
  36     svec<int> fixed;
  37     for (int i=0; i < ideals.sz(); i++) {
  38         assert(ideals[i]->hooke > 0);
  39         int l = col_id(ideals[i]->left);
  40         int r = col_id(ideals[i]->right);
  41         connected.connect(l,r);
  42     }
  43     for (int i = 0; i < cols.sz(); i++)
  44         if (cols[i].fixed())
  45             fixed.add(i);
  46     for (int i = 0; i < cols.sz(); i++) {
  47         bool c=false;
  48         for (int j =0; j<fixed.sz(); j++)
  49             c |=  connected.equiv(j,i);
  50         assert(c);
  51     }
  52 #endif
  53 }
  54
  55 bool
  56 Spacing_problem::check_constraints(Vector v) const
  57 {
  58     int dim=v.dim();
  59     // mtor << "checking solution " << v << '\n';
  60     for (int i=0; i < dim; i++) {
  61
  62         if (cols[i].fixed()&& abs(cols[i].fixed_position() - v(i)) > COLFUDGE) {
  63             return false;
  64         }
  65         if (!i)
  66             continue;
  67
  68         Real mindist=cols[i-1].minright()
  69             +cols[i].minleft();
  70
  71         // ugh... compares
  72         Real dif =v(i) - v(i-1)- mindist;
  73         bool b = (dif > - COLFUDGE);
  74
  75
  76 #if 1
  77         if (!b)
  78             return false;
  79
  80 #else
  81         mtor << "dif= "<<dif<<" fudge= " << COLFUDGE<< " dif >fudge= "<<
  82             b << "\n";
  83
  84         /* fucks up for unknown reasons */
  85         if (dif  < -COLFUDGE)
  86             return false;
  87 #endif
  88
  89     }
  90     return true;
  91 }
  92
  93 bool
  94 Spacing_problem::check_feasible() const
  95 {
  96     Vector sol(try_initial_solution());
  97     return check_constraints(sol);
  98 }
  99
 100 // generate a solution which obeys the min distances and fixed positions
 101 Vector
 102 Spacing_problem::try_initial_solution() const
 103 {
 104     int dim=cols.sz();
 105     Vector initsol(dim);
 106     for (int i=0; i < dim; i++) {
 107         if (cols[i].fixed()) {
 108             initsol(i)=cols[i].fixed_position();
 109         } else {
 110             Real mindist=cols[i-1].minright()
 111                 +cols[i].minleft();
 112             assert(mindist >= 0.0);
 113             initsol(i)=initsol(i-1)+mindist;
 114
 115             //nog niet
 116             //if (i>0)
 117             //  assert(initsol(i) > initsol(i-1));
 118         }
 119     }
 120
 121     return initsol;
 122 }
 123 Vector
 124 Spacing_problem::find_initial_solution() const
 125 {
 126     Vector v(try_initial_solution());
 127     assert(check_constraints(v));
 128     return v;
 129 }
 130
 131 // generate the matrices
 132 void
 133 Spacing_problem::make_matrices(Matrix &quad, Vector &lin, Real &c) const
 134 {
 135     quad.fill(0);
 136     lin.fill(0);
 137     c = 0;
 138     for (int j=0; j < ideals.sz(); j++){
 139         Idealspacing const*i=ideals[j];
 140         int l = col_id(i->left);
 141         int r = col_id(i->right);
 142
 143         quad(r,r) += i->hooke;
 144         quad(r,l) -= i->hooke;
 145         quad(l,r) -= i->hooke;
 146         quad(l,l) += i->hooke;
 147
 148         lin(r) -= i->space*i->hooke;
 149         lin(l) += i->space*i->hooke;
 150
 151         c += sqr(i->space);
 152     }
 153 }
 154
 155 // put the constraints into the LP problem
 156 void
 157 Spacing_problem::make_constraints(Mixed_qp& lp) const
 158 {
 159     int dim=cols.sz();
 160     for (int j=0; j < dim; j++) {
 161         Colinfo *c=&(cols[j]);
 162         if (c->fixed()) {
 163             lp.add_fixed_var(j,c->fixed_position());
 164         }
 165         if (j > 0){
 166             Vector c1(dim);
 167
 168             c1(j)=1.0 ;
 169             c1(j-1)=-1.0 ;
 170             lp.add_inequality_cons(c1, cols[j-1].minright() +
 171                                    cols[j].minleft());
 172         }
 173     }
 174 }
 175
 176 svec<Real>
 177 Spacing_problem::solve() const
 178 {
 179     print();
 180     OK();
 181     assert(check_feasible());
 182
 183     /* optimalisatiefunctie */
 184     Mixed_qp lp(cols.sz());
 185     make_matrices(lp.quad,lp.lin, lp.const_term);
 186     make_constraints(lp);
 187     Vector start=find_initial_solution();
 188     Vector sol(lp.solve(start));
 189     if (!check_constraints(sol)) {
 190         WARN << "solution doesn't satisfy constraints.\n" ;
 191     }
 192
 193
 194     svec<Real> posns(sol);
 195     posns.add(lp.eval(sol));
 196     return posns;
 197 }
 198
 199 /*
 200     add one column to the problem.
 201 */
 202 void
 203 Spacing_problem::add_column(const PCol *col, bool fixed, Real fixpos)
 204 {
 205     Colinfo c(col,(fixed)? &fixpos :  0);
 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 /*****************/
 249
 250 void
 251 Colinfo::print() const
 252 {
 253 #ifndef NPRINT
 254     mtor << "column { ";
 255     if (fixed())
 256         mtor << "fixed at " << fixed_position()<<", ";
 257     assert(pcol_);
 258     mtor << "[" << minleft() << ", " << minright() << "]";
 259     mtor <<"}\n";
 260 #endif
 261 }
 262
 263 Colinfo::Colinfo(Colinfo const&c)
 264 {
 265     fixpos = (c.fixpos)?new Real(*c.fixpos):0;
 266     pcol_ = c.pcol_;
 267     width = c.width;
 268 }
 269
 270 Colinfo::Colinfo(const PCol*col_p, const Real*fixed_r_p )
 271 {
 272     fixpos = (fixed_r_p)? new Real(*fixed_r_p) : 0;
 273     pcol_ = col_p;
 274     width = pcol_->width();
 275 }
 276
 277 Colinfo::~Colinfo()
 278 {
 279     delete fixpos;
 280 }
 281
 282 Colinfo::Colinfo()
 283 {
 284     pcol_=0;
 285     fixpos = 0;
 286 }
 287 void
 288 Colinfo::operator=(Colinfo const&c )
 289 {
 290     delete fixpos;
 291     fixpos = (c.fixpos)?new Real(*c.fixpos):0;
 292     pcol_ = c.pcol_;
 293     width = c.width;
 294 }