From 05f902d270633d1342bc0083b80812e87d594d94 Mon Sep 17 00:00:00 2001 From: fred Date: Sun, 24 Mar 2002 19:49:52 +0000 Subject: [PATCH] lilypond-0.0.78 --- flower/diagonal-storage.cc | 236 ++++++++++++++++++++++++++++++++++++ flower/include/matrix.hh | 24 +++- flower/matrix-debug.cc | 15 ++- flower/test/mat-test.cc | 60 +++++++++ lily/ineq-constrained-qp.cc | 198 ++++++++++++++++++++++++++++++ 5 files changed, 528 insertions(+), 5 deletions(-) create mode 100644 flower/diagonal-storage.cc create mode 100644 flower/test/mat-test.cc create mode 100644 lily/ineq-constrained-qp.cc diff --git a/flower/diagonal-storage.cc b/flower/diagonal-storage.cc new file mode 100644 index 0000000000..3cf2f93d20 --- /dev/null +++ b/flower/diagonal-storage.cc @@ -0,0 +1,236 @@ +/* + diagonal-storage.cc -- implement Diagonal_storage + + source file of the Flower Library + + (c) 1997 Han-Wen Nienhuys +*/ + + +#include "diagonal-storage.hh" + +int +Diagonal_storage::dim()const +{ + return band_.rows(); +} + +Diagonal_storage::Diagonal_storage() +{ +} + +int +Diagonal_storage::rows() const +{ + return band_.rows(); +} + +int +Diagonal_storage::cols() const +{ + return band_.rows(); +} + +int +Diagonal_storage::band_size_i()const +{ + return (band_.cols()-1)/2; +} + +void +Diagonal_storage::set_band_size(int s) +{ + Full_storage f(dim(), 2*s+1); + for (int i=0; i < dim(); i++) { + int k=-s; + for ( ; k < -band_size_i(); k++) + f.elem(i,k + s) = 0.0; + for ( ; k <= band_size_i()&& k<=s ; k++) + f.elem(i, k + s) = band_.elem(i,k+ band_size_i()); + for( ; k <= s; k++) + f.elem(i, k + s) =0.0; + } + + band_ = f; +} + + + +/* + any takers? + */ +void +Diagonal_storage::insert_row(int ) +{ + assert(false); +} + +void +Diagonal_storage::delete_row(int ) +{ + assert(false); +} + +void +Diagonal_storage::resize(int,int) +{ +} + +void +Diagonal_storage::resize(int) +{ +} + +void +Diagonal_storage::delete_column(int ) +{ + assert(false); +} + +Diagonal_storage::~Diagonal_storage() +{ +} + + +bool +Diagonal_storage::band_elt_b(int i,int j)const +{ + return abs(i-j) <= band_size_i(); +} + +void +Diagonal_storage::assert_valid(int i,int j )const +{ + assert( band_elt_b(i,j) ); + assert( i >=0 && j >=0 && i < dim() && j < dim()); +} + + +void +Diagonal_storage::resize_dim(int d) +{ + Full_storage f(d, 2*band_size_i()+1); + for (int i=0; i < d&& i < dim(); i++) { + for ( int k=0; k < 2*band_size_i(); k++) + f.elem(i,k) = elem(i,k); + } + + band_ = f; +} + + + +bool +Diagonal_storage::mult_ok(int i,int )const +{ + return i < dim(); +} + +void +Diagonal_storage::mult_next(int &i, int &j)const +{ + j++; + if ( j < i - band_size_i() ) + j = i- band_size_i(); + if ( j > i + band_size_i() || j >= dim() ) { + i++; + j = i - band_size_i(); + if (j < 0) + j=0; + } +} + +bool +Diagonal_storage::trans_ok(int ,int j)const +{ + return j < dim(); +} + +void +Diagonal_storage::trans_next(int &i, int& j)const +{ + i++; + if ( i < j - band_size_i()) + i = j-band_size_i(); + + if ( i >= dim() || i > j + band_size_i() ) { + j++; + i = j - band_size_i(); + if (i < 0) + i=0; + } +} + +static Real nul_entry=0.0; + +Real +Diagonal_storage::elem(int i, int j)const +{ + if (abs ( i-j ) > band_size_i()) + return 0; + else + return band_.elem(i, j - i +band_size_i()); +} + +Real & +Diagonal_storage::elem(int i, int j) +{ + /* + if this fails, the previous call fucked up + */ + assert(nul_entry); + if (abs ( i-j ) > band_size_i()) + return nul_entry; + else + return band_.elem(i, j - i + band_size_i()); +} + +/* + Hairy routine to try to save some fp-ops + */ + +bool +Diagonal_storage::try_right_multiply(Matrix_storage*dest, + const Matrix_storage*right)const +{ + if ( right->name() != Diagonal_storage::static_name() ) + return false; + + const Diagonal_storage* diag = (Diagonal_storage const*)right; + int band2 = diag->band_size_i(); + int n = dim(); + /* + should check if dest is a Diagonal_storage of sufficient size too. + */ + for (int i=0; i < n; i++) { + for (int j = 0; j < n; j++) { + int startk = i - band_size_i() >? 0 >? j - band2; + int stopk = i + band_size_i() elem(relk, j); + dest->elem(i, j) = sum; + + } + } + return true; +} + +IMPLEMENT_IS_TYPE_B1(Diagonal_storage, Matrix_storage); + + +Diagonal_storage::Diagonal_storage(Matrix_storage*stor_l, int band_i) +{ + set_band_size(band_i); + resize_dim(stor_l->dim()); + + for ( int i=0,j=0; mult_ok(i,j); mult_next(i,j)) + band_.elem(i, j + band_i -i ) = stor_l->elem(i,j); +} + +void +Diagonal_storage::OK() const +{ + band_.OK(); +} diff --git a/flower/include/matrix.hh b/flower/include/matrix.hh index df013e0849..b6e084efda 100644 --- a/flower/include/matrix.hh +++ b/flower/include/matrix.hh @@ -24,8 +24,12 @@ class Matrix { + friend Matrix operator *(Matrix const &m1, Matrix const &m2); + +protected: Matrix_storage *dat; - + void set(Matrix_storage*); + Matrix(Matrix_storage*); public: void OK() const { dat->OK(); } int cols() const { return dat->cols(); } @@ -37,7 +41,16 @@ public: */ int dim() const; - // Matrix() { dat = 0; } + /** + the band size of the matrix. + @ret + + 0 <= band_i() <= dim + */ + int band_i() const; + bool band_b()const; + void set_full()const; + void try_set_band()const; ~Matrix() { delete dat; } /// set entries to r @@ -76,7 +89,7 @@ public: square n matrix, initialised to null */ Matrix(int n); - + /** n x m matrix, init to 0 */ @@ -140,4 +153,9 @@ inline Matrix operator -(Matrix m1,const Matrix m2) m1 -= m2; return m1; } +inline Matrix operator +(Matrix m1,const Matrix m2) +{ + m1 += m2; + return m1; +} #endif diff --git a/flower/matrix-debug.cc b/flower/matrix-debug.cc index 47a66896a5..099a05beb9 100644 --- a/flower/matrix-debug.cc +++ b/flower/matrix-debug.cc @@ -1,15 +1,26 @@ +/* + matrix-debug.cc -- implement Matrix print routines + + source file of the Flower Library + + (c) 1997 Han-Wen Nienhuys +*/ + + #include "flower-debug.hh" #include "matrix.hh" +#include "matrix-storage.hh" Matrix::operator String() const { String s; #ifndef NPRINT - s="matrix {\n"; + Matrix_storage const * stor_c_l = dat; + s=String("matrix { (") + dat->name() + ")\n"; for (int i=0; i< rows(); i++){ for (int j = 0; j < cols(); j++) { - s+= String(dat->elem(i,j), "%6f "); + s+= String(stor_c_l->elem(i,j), "%6f "); } s+="\n"; } diff --git a/flower/test/mat-test.cc b/flower/test/mat-test.cc new file mode 100644 index 0000000000..f7ce01bfc5 --- /dev/null +++ b/flower/test/mat-test.cc @@ -0,0 +1,60 @@ +/* + mat-test.cc -- test Matrix + + source file of the Flower Library + + (c) 1997 Han-Wen Nienhuys +*/ + +#include +#include "matrix.hh" +#include "string.hh" +#include "flower-test.hh" +#include "choleski.hh" + +void +matrix() +{ + int N=10; + Matrix m(N,N), q(N,N); + Vector v(N); + + for (int i=0; i < N; i++) { + v(i) =i; + for (int j=0; j < N; j++) { + m(i,j) = i+j; + q(i,j) = (abs(i-j) > 3) ?0 :i-j; + } + } + + cout << "v: " << String(v); + cout << "m: " << String(m ); + cout << "q: " << String(q); + cout << "m*q; " << String(m*q); + cout << "m*m: " << String(m*m); + m.OK(); + cout << "m: " << String(m); + cout << "q.band " << q.band_i() << endl; + q.try_set_band(); + cout << "q(B): " << q; + q.OK(); + Matrix sum(q+q); + cout << "q + q " << sum; + q.OK(); + cout << "q*q: " << q*q; + q.OK(); + + Matrix hilbert(N,N), h2(hilbert); + for (int i=0; i < N; i++) { + for (int j=0; j < N; j++) { + hilbert(i,j) = 1/(i+j+1); + h2 (i,j) = (abs(i-j) > 3) ?0 : hilbert(i,j); + } + } + h2.try_set_band(); + Choleski_decomposition ch(h2); + cout << "red Hilbert " << h2; + cout << "choleski " << ch.L; +} + +ADD_TEST(matrix); diff --git a/lily/ineq-constrained-qp.cc b/lily/ineq-constrained-qp.cc new file mode 100644 index 0000000000..8fed514623 --- /dev/null +++ b/lily/ineq-constrained-qp.cc @@ -0,0 +1,198 @@ +/* + ineq-constrained-qp.cc -- implement Ineq_constrained_qp + + source file of the GNU LilyPond music typesetter + + (c) 1996,1997 Han-Wen Nienhuys +*/ +#include "ineq-constrained-qp.hh" +#include "qlpsolve.hh" +#include "const.hh" +#include "debug.hh" +/* + MAy be this also should go into a library + */ + +const int MAXITER=100; // qlpsolve.hh + +/* + assume x(idx) == value, and adjust constraints, lin and quad accordingly + + TODO: add const_term + */ +void +Ineq_constrained_qp::eliminate_var(int idx, Real value) +{ + Vector row(quad.row(idx)); + row*= value; + + quad.delete_row(idx); + + quad.delete_column(idx); + + lin.del(idx); + row.del(idx); + lin +=row ; + + for (int i=0; i < cons.size(); i++) { + consrhs[i] -= cons[i](idx) *value; + cons[i].del(idx); + } +} + +void +Ineq_constrained_qp::add_inequality_cons(Vector c, double r) +{ + cons.push(c); + consrhs.push(r); +} + +Ineq_constrained_qp::Ineq_constrained_qp(int novars): + quad(novars), + lin(novars), + const_term (0.0) +{ +} + +void +Ineq_constrained_qp::OK() const +{ +#if !defined(NDEBUG) && defined(PARANOID) + assert(cons.size() == consrhs.size()); + Matrix Qdif= quad - quad.transposed(); + assert(Qdif.norm()/quad.norm() < EPS); +#endif +} + + +Real +Ineq_constrained_qp::eval (Vector v) +{ + return v * quad * v + lin * v + const_term; +} + + +int +min_elt_index(Vector v) +{ + Real m=INFTY; int idx=-1; + for (int i = 0; i < v.dim(); i++){ + if (v(i) < m) { + idx = i; + m = v(i); + } + assert(v(i) <= INFTY); + } + return idx; +} + + +/**the numerical solving. Mordecai Avriel, Nonlinear Programming: analysis and methods (1976) + Prentice Hall. + + Section 13.3 + + This is a "projected gradient" algorithm. Starting from a point x + the next point is found in a direction determined by projecting + the gradient onto the active constraints. (well, not really the + gradient. The optimal solution obeying the active constraints is + tried. This is why H = Q^-1 in initialisation) ) + + + */ +Vector +Ineq_constrained_qp::solve(Vector start) const +{ + if (!dim()) + return Vector(0); + // experimental +// quad.try_set_band(); + + Active_constraints act(this); + + + act.OK(); + + + Vector x(start); + Vector gradient=quad*x+lin; +// Real fvalue = x*quad*x/2 + lin*x + const_term; +// it's no use. + + Vector last_gradient(gradient); + int iterations=0; + + while (iterations++ < MAXITER) { + Vector direction= - act.find_active_optimum(gradient); + + mtor << "gradient "<< gradient<< "\ndirection " << direction<<"\n"; + + if (direction.norm() > EPS) { + mtor << act.status() << '\n'; + + Real minalf = INFTY; + + Inactive_iter minidx(act); + + + /* + we know the optimum on this "hyperplane". Check if we + bump into the edges of the simplex + */ + + for (Inactive_iter ia(act); ia.ok(); ia++) { + + if (ia.vec() * direction >= 0) + continue; + Real alfa= - (ia.vec()*x - ia.rhs())/ + (ia.vec()*direction); + + if (minalf > alfa) { + minidx = ia; + minalf = alfa; + } + } + Real unbounded_alfa = 1.0; + Real optimal_step = min(minalf, unbounded_alfa); + + Vector deltax=direction * optimal_step; + x += deltax; + gradient += optimal_step * (quad * deltax); + + mtor << "step = " << optimal_step<< " (|dx| = " << + deltax.norm() << ")\n"; + + if (minalf < unbounded_alfa) { + /* bumped into an edge. try again, in smaller space. */ + act.add(minidx.idx()); + mtor << "adding cons "<< minidx.idx()<<'\n'; + continue; + } + /*ASSERT: we are at optimal solution for this "plane"*/ + + + } + + Vector lagrange_mult=act.get_lagrange(gradient); + int m= min_elt_index(lagrange_mult); + + if (m>=0 && lagrange_mult(m) > 0) { + break; // optimal sol. + } else if (m<0) { + assert(gradient.norm() < EPS) ; + + break; + } + + mtor << "dropping cons " << m<<'\n'; + act.drop(m); + } + if (iterations >= MAXITER) + WARN<<"didn't converge!\n"; + + mtor << ": found " << x<<" in " << iterations <<" iterations\n"; + assert_solution(x); + return x; +} + + -- 2.39.5