--- /dev/null
+/*
+ diagonal-storage.cc -- implement Diagonal_storage
+
+ source file of the Flower Library
+
+ (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+
+#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() <? n-1 <? j + band2;
+ int relk = startk + band_size_i() -i;
+ Real sum =0.0;
+ for ( int k = startk; k <= stopk; k++)
+ sum += band_.elem(i, relk) * diag->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();
+}
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(); }
*/
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
square n matrix, initialised to null
*/
Matrix(int n);
-
+
/**
n x m matrix, init to 0
*/
m1 -= m2;
return m1;
}
+inline Matrix operator +(Matrix m1,const Matrix m2)
+{
+ m1 += m2;
+ return m1;
+}
#endif
+/*
+ matrix-debug.cc -- implement Matrix print routines
+
+ source file of the Flower Library
+
+ (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+
#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";
}
--- /dev/null
+/*
+ mat-test.cc -- test Matrix
+
+ source file of the Flower Library
+
+ (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+#include <iostream.h>
+#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);
--- /dev/null
+/*
+ ineq-constrained-qp.cc -- implement Ineq_constrained_qp
+
+ source file of the GNU LilyPond music typesetter
+
+ (c) 1996,1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+#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;
+}
+
+
projects / lilypond.git / commitdiff