lily/qlp.cc

   1 /*
   2   qlp.cc -- implement Mixed_qp
   3
   4   source file of the GNU LilyPond music typesetter
   5
   6   (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
   7 */
   8
   9 #include "debug.hh"
  10 #include "const.hh"
  11 #include "qlp.hh"
  12 #include "choleski.hh"
  13
  14
  15 void
  16 Mixed_qp::add_equality_cons(Vector , double )
  17 {
  18     assert(false);
  19 }
  20
  21 void
  22 Mixed_qp::add_fixed_var(int i, Real r)
  23 {
  24     eq_cons.push(i);
  25     eq_consrhs.push(r);
  26 }
  27
  28
  29 /**
  30     eliminate appropriate variables, until we have a Ineq_constrained_qp
  31     then solve that.
  32
  33     PRE
  34     cons should be ascending
  35     */
  36 Vector
  37 Mixed_qp::solve(Vector start) const
  38 {
  39     if (!dim())
  40         return Vector(0);
  41
  42     print();
  43     Ineq_constrained_qp pure(*this);
  44
  45     for  (int i= eq_cons.size()-1; i>=0; i--) {
  46         pure.eliminate_var(eq_cons[i], eq_consrhs[i]);
  47         start.del(eq_cons[i]);
  48     }
  49     Vector sol = pure.solve(start);
  50     for (int i= 0; i < eq_cons.size(); i++) {
  51         sol.insert( eq_consrhs[i],eq_cons[i]);
  52     }
  53     return sol;
  54 }
  55
  56
  57 void
  58 Ineq_constrained_qp::assert_solution(Vector sol) const
  59 {
  60     Array<int> binding;
  61     for (int i=0; i < cons.size(); i++) {
  62         Real R=cons[i] * sol- consrhs[i];
  63         assert(R> -EPS);
  64         if (R < EPS)
  65             binding.push(i);
  66     }
  67     // KKT check...
  68     // todo
  69 }
  70
  71 void
  72 Ineq_constrained_qp::print() const
  73 {
  74 #ifndef NPRINT
  75     mtor << "Quad " << quad;
  76     mtor << "lin " << lin <<"\n"
  77         << "const " << const_term<<"\n";
  78     for (int i=0; i < cons.size(); i++) {
  79         mtor << "constraint["<<i<<"]: " << cons[i] << " >= " << consrhs[i];
  80         mtor << "\n";
  81     }
  82 #endif
  83 }
  84
  85 Mixed_qp::Mixed_qp(int n)
  86     : Ineq_constrained_qp(n)
  87 {
  88 }
  89
  90 void
  91 Mixed_qp::OK() const
  92 {
  93 #ifndef NDEBUG
  94     Ineq_constrained_qp::OK();
  95     assert(eq_consrhs.size() == eq_cons.size());
  96 #endif
  97 }
  98
  99 void
 100 Mixed_qp::print() const
 101 {
 102 #ifndef NPRINT
 103     Ineq_constrained_qp::print();
 104     for (int i=0; i < eq_cons.size(); i++) {
 105         mtor << "eq cons "<<i<<": x["<<eq_cons[i]<<"] == " << eq_consrhs[i]<<"\n";
 106     }
 107 #endif
 108 }
 109