+++ /dev/null
-/*
- choleski.cc -- implement Choleski_decomposition
-
- source file of the Flower Library
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-#include "choleski.hh"
-const Real EPS = 1e-7; // so sue me. Hard coded
-
-// for testing new Matrix_storage.
-//#define PARANOID
-
-void
-Choleski_decomposition::full_matrix_solve (Vector &out, Vector const &rhs) const
-{
- int n= rhs.dim();
- assert (n == L.dim());
- Vector y;
- y.set_dim (n);
- out.set_dim (n);
-
- // forward substitution
- for (int i=0; i < n; i++)
- {
- Real sum (0.0);
- for (int j=0; j < i; j++)
- sum += y (j) * L(i,j);
- y (i) = (rhs (i) - sum)/L(i,i);
- }
-
- for (int i=0; i < n; i++)
- y (i) /= D(i);
-
- // backward subst
- Vector &x (out); // using input as return val.
- for (int i=n-1; i >= 0; i--)
- {
- Real sum (0.0);
- for (int j=i+1; j < n; j++)
- sum += L(j,i)*x (j);
- x (i) = (y (i) - sum)/L(i,i);
- }
-}
-
-void
-Choleski_decomposition::band_matrix_solve (Vector &out, Vector const &rhs) const
-{
- int n= rhs.dim();
- int b = L.band_i();
- assert (n == L.dim());
-
- out.set_dim (n);
-
- Vector y;
- y.set_dim (n);
-
- // forward substitution
- for (int i=0; i < n; i++)
- {
- Real sum (0.0);
- for (int j= 0 >? i - b; j < i; j++)
- sum += y (j) * L(i,j);
- y (i) = (rhs (i) - sum)/L(i,i);
- }
- for (int i=0; i < n; i++)
- y (i) /= D(i);
-
- // backward subst
- Vector &x (out); // using input as return val.
- for (int i=n-1; i >= 0; i--)
- {
- Real sum (0.0);
- for (int j=i+1; j <= i + b&&j < n ; j++)
- sum += L(j,i)*x (j);
- x (i) = (y (i) - sum)/L(i,i);
- }
-}
-
-void
-Choleski_decomposition::solve (Vector &x, Vector const &rhs) const
-{
- if (band_b_)
- {
- band_matrix_solve (x,rhs);
- }
- else
- full_matrix_solve (x,rhs);
-}
-
-Vector
-Choleski_decomposition::solve (Vector rhs) const
-{
- Vector r;
- solve (r, rhs);
- return r;
-}
-
-void
-Choleski_decomposition::full_matrix_decompose (Matrix const & P)
-{
-
- int n = P.dim();
- L.unit();
- for (int k= 0; k < n; k++)
- {
- for (int j = 0; j < k; j++)
- {
- Real sum (0.0);
- for (int l=0; l < j; l++)
- sum += L(k,l)*L(j,l)*D(l);
- L(k,j) = (P(k,j) - sum)/D(j);
- }
- Real sum=0.0;
-
- for (int l=0; l < k; l++)
- sum += sqr (L(k,l))*D(l);
- Real d = P(k,k) - sum;
- D(k) = d;
- }
-
-}
-
-void
-Choleski_decomposition::band_matrix_decompose (Matrix const &P)
-{
- int n = P.dim();
- int b = P.band_i();
- L.unit();
-
- for (int i= 0; i < n; i++)
- {
- for (int j = 0 >? i - b; j < i; j++)
- {
- Real sum (0.0);
- for (int l=0 >? i - b; l < j; l++)
- sum += L(i,l)*L(j,l)*D(l);
- L(i,j) = (P(i,j) - sum)/D(j);
- }
- Real sum=0.0;
-
- for (int l=0 >? i - b; l < i; l++)
- sum += sqr (L(i,l))*D(l);
- Real d = P(i,i) - sum;
- D(i) = d;
- }
- L.set_band();
- band_b_ = true;
-}
-
-
-
-
-/*
- Standard matrix algorithm.
- */
-
-Choleski_decomposition::Choleski_decomposition (Matrix const & P)
- : L(P.dim()), D(P.dim ())
-{
-#ifdef PARANOID
- assert ((P-P.transposed()).norm ()/P.norm () < EPS);
-#endif
- band_b_ = false;
-
- int b = P.calc_band_i ();
-
- if (b <= P.dim ()/2)
- band_matrix_decompose (P);
- else
- full_matrix_decompose (P);
-
-#ifdef PARANOID
- assert ((original()-P).norm () / P.norm () < EPS);
-#endif
-}
-
-Matrix
-Choleski_decomposition::original() const
-{
- Matrix T(L.dim());
- T.set_diag (D);
- return L*T*L.transposed();
-}
-
-Matrix
-Choleski_decomposition::inverse() const
-{
- int n=L.dim();
- Matrix invm (n);
- Vector e_i (n);
- Vector inv (n);
- for (int i = 0; i < n; i++)
- {
- e_i.set_unit (i);
- solve (inv, e_i);
- for (int j = 0 ; j<n; j++)
- invm (i,j) = inv (j);
- }
-
-#ifdef PARANOID
- Matrix I1(n), I2(original());
- I1.unit();
- assert ((I1-I2*invm).norm()/I2.norm () < EPS);
-#endif
-
- return invm;
-}
+++ /dev/null
-/*
- full-storage.cc -- implement Full_storage
-
- source file of the Flower Library
-
- (c) 1996, 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-#include "full-storage.hh"
-
-
-void
-Full_storage::operator=(Full_storage const &fs)
-{
- resize (fs.height_i_, fs.width_i_);
- OK();
- fs.OK();
- for (int i=0; i<height_i_; i++)
- for (int j=0; j<width_i_; j++)
- els_p_p_[i][j]= fs.els_p_p_[i][j];
-}
-
-
-void
-Full_storage::OK() const
-{
-#ifndef NDEBUG
- assert (max_height_i_ >= height_i_ && max_width_i_ >= width_i_);
- assert (height_i_ >= 0 && width_i_ >= 0);
- assert (els_p_p_||!max_height_i_);
-#endif
-}
-
-
-
-
-Full_storage::~Full_storage()
-{
- for (int i=0; i < max_height_i_; i++)
- delete [] els_p_p_[i];
- delete[] els_p_p_;
-}
-
-void
-Full_storage::resize (int rows, int cols)
-{
- OK();
- resize_cols (rows);
- resize_rows (cols);
- band_i_ = rows >? cols;
-}
-
-void
-Full_storage::delete_column (int k)
-{
- assert (0 <= k &&k<width_i_);
- for (int i=0; i< height_i_ ; i++)
- for (int j=k+1; j <width_i_; j++)
- els_p_p_[i][j-1]=els_p_p_[i][j];
- width_i_--;
-}
-
-
-void
-Full_storage::delete_row (int k)
-{
- assert (0 <= k &&k<height_i_);
- for (int i=k+1; i < height_i_ ; i++)
- for (int j=0; j < width_i_; j++)
- els_p_p_[i-1][j]=els_p_p_[i][j];
- height_i_--;
-}
-
-
-
-void
-Full_storage::insert_row (int k)
-{
- assert (0 <= k&& k <=height_i_);
- resize_cols (height_i_+1);
- for (int i=height_i_-1; i > k ; i--)
- for (int j=0; j <width_i_; j++)
- els_p_p_[i][j]=els_p_p_[i-1][j];
-
-}
-Array<Real>
-Full_storage::row (int n) const
-{
- Array<Real> r;
- for (int j = 0; j < cols(); j++)
- r.push (elem (n,j));
- return r;
-}
-
-Array<Real>
-Full_storage::column (int n) const
-{
- Array<Real> r;
- for (int i = 0; i < rows(); i++)
- r.push (elem (i,n));
- return r;
-}
-
-void
-Full_storage::set_size (int rows, int cols)
-{
- resize (rows,cols);
-}
-
-void
-Full_storage::set_size (int rows)
-{
-
- resize (rows);
-}
-
-
-
-void
-Full_storage::resize_cols (int newh)
-{
- if (newh <= max_height_i_)
- {
- height_i_=newh;
- return;
- }
-
- Real ** newa=new Real*[newh];
- int j=0;
- for (; j < height_i_; j++)
- newa[j] = els_p_p_[j];
- for (; j < newh; j++)
- newa[j] = new Real[max_width_i_];
- delete[] els_p_p_;
- els_p_p_=newa;
-
- height_i_ = max_height_i_ = newh;
-}
-
-
-
-void
-Full_storage::resize_rows (int neww)
-{
- if (neww <= max_width_i_)
- {
- width_i_=neww;
- return;
- }
- for (int i=0; i < max_height_i_ ; i++)
- {
- Real* newa = new Real[neww];
- for (int k=0; k < width_i_; k++)
- newa[k] = els_p_p_[i][k];
-
- delete[] els_p_p_[i];
- els_p_p_[i] = newa;
- }
- width_i_ = max_width_i_ = neww;
-}
-
-#ifdef INLINE
-#undef INLINE
-#endif
-#define INLINE
-
-
-#include "full-storage.icc"
+++ /dev/null
-#ifndef CHOLESKI_HH
-#define CHOLESKI_HH
-
-#include "matrix.hh"
-
-/**
- Choleski decomposition of a matrix
- structure for using the LU decomposition of a positive definite matrix.
-
- #P# is split into
-
- LD transpose (L)
- */
-struct Choleski_decomposition {
-
- /// lower triangle of Choleski decomposition
- Matrix L;
-
- bool band_b_;
-
- /// diagonal
- Vector D;
-
- /** Create decomposition of P.
- PRE
- P needs to be symmetric positive definite
- */
-
- Choleski_decomposition (Matrix const &P);
-
- /**
- solve Px = rhs
- */
- Vector solve (Vector rhs) const;
- void solve (Vector &dest, Vector const &rhs) const;
- Vector operator * (Vector rhs) const { return solve (rhs); }
- /**
- return the inverse of the matrix P.
- */
- Matrix inverse() const;
- /**
- return P, calc'ed from L and D
- */
- Matrix original() const;
-private:
- void full_matrix_solve (Vector &,Vector const&) const;
- void band_matrix_solve (Vector &, Vector const&) const;
- void full_matrix_decompose (Matrix const & P);
- void band_matrix_decompose (Matrix const &P);
-
-};
-#endif
+++ /dev/null
-/*
- full-storage.hh -- declare Full_storage
-
- source file of the Flower Library
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#ifndef FULL_STORAGE_HH
-#define FULL_STORAGE_HH
-
-#include "array.hh"
-#include "real.hh"
-
-#ifndef INLINE
-#define INLINE inline
-#endif
-
-/// simplest matrix storage. refer to its baseclass for the doco.
-class Full_storage
-{
- /// height, width
- int height_i_;
- int width_i_;
- /// maxima.
- int max_height_i_;
- int max_width_i_;
-
- /// the storage
- Real** els_p_p_;
-
- INLINE void init() ;
-
- INLINE bool valid (int i, int j) const ;
-
- void resize_rows (int);
- void resize_cols (int);
-
-public:
- int band_i_; // ugh
-
- /// check invariants
- void OK() const;
-
- /// height of matrix
- INLINE int rows() const;
-
- /// width of matrix
- INLINE int cols() const;
-
- /// size if square
- INLINE int dim() const;
-
- /** set the size. contents lost.
- PRE
- i >=0, j>=0
- */
- void set_size (int rows, int cols) ;
-
-
- /**set the size to square dimen. contents lost
- PRE
- i>=0
- */
- void set_size (int i) ;
- /**set the size to i.
-
- keep contents. If enlarged contents unspecified
-
- PRE
- i>=0, j>=0
-
- */
- void resize (int rows, int cols);
-
- /**
- set the size to square dimen. contents kept
- Keep contents. If enlarged contents are unspecified
-
- PRE
- i>=0
- */
- void resize (int i);
-
-
- /**
- access an element.
-
- Generate an errormessage, if this happens
- in the 0-part of a sparse matrix.
- */
-
- INLINE Real& elem (int i,int j);
-
- /// access a element, no modify
- INLINE Real elem (int i, int j) const;
-
- Array<Real> row (int i) const ;
- Array<Real> column (int j) const;
-
-
- /**
- add a row to the matrix before row k. Contents
- of added row are unspecified
-
- 0 <= k <= rows()
- */
- void insert_row (int k);
-
-
- /**
- delete a row from this matrix.
-
- PRE
- 0 <= k < rows();
- */
- void delete_row (int k);
- void delete_column (int k);
-
-
-
- /**
- at end of matrix?. when doing loop
-
- for (i=0; i<h; i++)
- for (j=0; j<w; j++)
- ..
-
- */
- INLINE bool mult_ok (int i, int j) const;
-
- /**
- walk through matrix (regular multiply).
- get next j for row i, or get next row i and reset j.
- this will make sparse matrix implementation easy.
-
- PRE
- mult_ok (i,j)
- */
- INLINE void mult_next (int &i, int &j) const;
-
- /**
- valid matrix entry. return false if at end of row
- */
- INLINE bool trans_ok (int i, int j) const;
-
- /**
- walk through matrix (transposed multiply).
- Get next i (for column j)
-
- PRE
- ver_ok (i,j)
- */
-
- INLINE void trans_next (int &i, int &j) const;
-
- INLINE Full_storage();
- INLINE Full_storage (int i, int j);
- INLINE Full_storage (Full_storage const&);
- INLINE Full_storage (int i);
- void operator=(Full_storage const &);
-
-
- ~Full_storage();
-};
-
-#include "full-storage.icc"
-
-
-#endif // FULL_STORAGE_HH
+++ /dev/null
-/*
- matrix-storage.hh -- declare Matrix_storage
-
- source file of the Flower Library
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#ifndef MATRIX_STORAGE_HH
-#define MATRIX_STORAGE_HH
-
-#include "array.hh"
-#include "real.hh"
-
-/**
-
- base class for interface with matrix storageclasses. There are no
- iterators for matrixclasses, since matrices are (like arrays)
- explicitly int-indexed.
-
- Iteration is provided by *_next, *_ok, which update and check both
- index variables simultaneously.
-
- TODO
- determine type of product matrix.
-
-*/
-class Matrix_storage {
-
-
-public:
- /// check invariants
- void OK() const;
-
- /// height of matrix
- int rows() const;
-
- /// width of matrix
- int cols() const;
-
- /// size if square
- int dim() const;
-
- /** set the size. contents lost.
- PRE
- i >=0, j>=0
- */
- void set_size (int rows, int cols) ;
-
-
- /**set the size to square dimen. contents lost
- PRE
- i>=0
- */
- void set_size (int i) ;
- /**set the size to i.
-
- keep contents. If enlarged contents unspecified
-
- PRE
- i>=0, j>=0
-
- */
- void resize (int rows, int cols);
-
- /**
- set the size to square dimen. contents kept
- Keep contents. If enlarged contents are unspecified
-
- PRE
- i>=0
- */
- void resize (int i);
-
-
- /**
- access an element.
-
- Generate an errormessage, if this happens
- in the 0-part of a sparse matrix.
- */
-
- Real& elem (int i,int j);
-
- /// access a element, no modify
- Real elem (int i, int j) const;
-
- Array<Real> row (int i) const ;
- Array<Real> column (int j) const;
-
-
- /**
- add a row to the matrix before row k. Contents
- of added row are unspecified
-
- 0 <= k <= rows()
- */
- void insert_row (int k);
-
-
- /**
- delete a row from this matrix.
-
- PRE
- 0 <= k < rows();
- */
- void delete_row (int k);
- void delete_column (int k);
- ~Matrix_storage() { }
-
-
- /**
- at end of matrix?. when doing loop
-
- for (i=0; i<h; i++)
- for (j=0; j<w; j++)
- ..
-
- */
- bool mult_ok (int i, int j) const;
-
- /**
- walk through matrix (regular multiply).
- get next j for row i, or get next row i and reset j.
- this will make sparse matrix implementation easy.
-
- PRE
- mult_ok (i,j)
- */
- void mult_next (int &i, int &j) const;
-
- /**
- valid matrix entry. return false if at end of row
- */
- bool trans_ok (int i, int j) const;
-
- /**
- walk through matrix (transposed multiply).
- Get next i (for column j)
-
- PRE
- ver_ok (i,j)
- */
-
- void trans_next (int &i, int &j) const;
-
- /// generate a "Full_storage" matrix
- static void set_band (Matrix_storage*&, int band);
-};
-
-#endif // MATRIX_STORAGE_HH
-
+++ /dev/null
-/*
- matrix.hh -- declare Matrix
-
- source file of the Flower Library
-
- (c) 1996, 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#ifndef MATRIX_HH
-#define MATRIX_HH
-
-#include "full-storage.hh"
-#include "vector.hh"
-#include "offset.hh"
-
-/** a Real matrix. This is a class for a nonsquare block of #Real#s. The
- implementation of sparse matrices is done in the appropriate #smat#
- class. Matrix only does the mathematical actions (adding,
- multiplying, etc.)
-
-
- TODO
- implement ref counting? */
-
-
-class Matrix {
- friend Matrix operator *(Matrix const &m1, Matrix const &m2);
-
-protected:
- Full_storage *dat_;
- void set (Matrix_storage*);
-public:
- void OK() const { dat_->OK(); }
- void set_band ();
- int calc_band_i () const;
- int cols() const { return dat_->cols (); }
- int rows() const { return dat_->rows (); }
-
- /** return the size of a matrix.
- PRE
- the matrix needs to be square.
- */
- int dim() const;
-
- /**
- the band size of the matrix.
- @ret
-
- 0 <= band_i() <= dim
- */
- int band_i() const;
-
- /// set entries to r
- void fill (Real r);
-
- /// set diagonal to d
- void set_diag (Real d);
-
- void set_diag (Vector d);
- /// set unit matrix
- void unit() { set_diag (1.0); }
-
- void operator+=(Matrix const &m);
- void operator-=(Matrix const &m);
- void operator*=(Real a);
- void operator/=(Real a) { (*this) *= 1/a; }
-
- /** add a row.
- add a row to the matrix before row k
-
- PRE
- v.dim() == cols ()
- 0 <= k <= rows()
- */
- void insert_row (Vector v,int k);
- /** .
- delete a row from this matrix.
-
- PRE
- 0 <= k < rows();
- */
- void delete_row (int k) { dat_->delete_row (k); }
- void delete_column (int k) { dat_->delete_column (k); }
-
- /**
- square n matrix, initialised to null
- */
- Matrix (int n);
-
- /**
- n x m matrix, init to 0
- */
- Matrix (int n, int m);
- Matrix (Matrix const &m);
-
- /// dyadic product: v * w.transpose
- Matrix (Vector v, Vector w);
- void operator=(Matrix const &m);
-
- /// access an element
- Real operator()(int i,int j) const { return dat_->elem (i,j); }
-
- /// access an element
- Real &operator()(int i, int j) { return dat_->elem (i,j); }
-
- /// Matrix multiply with vec (from right)
- Vector operator *(Vector const &v) const;
-
- /// set this to m1*m2.
- void set_product (Matrix const &m1, Matrix const &m2);
-
- Vector left_multiply (Vector const &) const;
-
- Matrix operator-() const;
-
- /// transpose this.
- void transpose();
-
- /// return a transposed copy.
- Matrix transposed() const ;
-
- Real norm() const;
- /** swap.
- PRE
- 0 <= c1,c2 < cols()
- */
- void swap_columns (int c1, int c2);
-
- /** swap.
- PRE
- 0 <= c1,c2 < rows()
- */
- void swap_rows (int c1, int c2);
-
-
- Vector row (int) const;
- Vector col (int) const;
-
- String str () const;
- void print() const;
- ~Matrix ();
-};
-
-inline Vector
-operator *(Vector &v, Matrix const & m) { return m.left_multiply (v); }
-Matrix operator *(Matrix const & m1,Matrix const &m2);
-Matrix operator /(Matrix const &m1,Real a);
-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
+++ /dev/null
-#ifndef VECTOR_HH
-#define VECTOR_HH
-
-#include <math.h>
-#include "real.hh"
-#include "array.hh"
-
-class Dstream;
-class String;
-
-/** a row of numbers.
- a vector. Storage is handled in Array, Vector only does the mathematics.
- */
-class Vector {
- Array<Real> dat;
-public:
- void OK() const { dat.OK();}
- int dim() const { return dat.size (); }
- void fill (Real r) {
- for (int i=0; i < dim(); i++)
- dat[i] =r;
- }
-
- Vector() { }
- Vector (Array<Real> d);
- Vector (Vector const &n);
-
- Real &operator()(int i) { return dat[i]; }
- Real operator()(int i) const { return dat[i]; }
- Real elem (int i) { return dat[i]; }
- Vector (int n) {
- dat.set_size (n);
- fill (0);
- }
- void set_dim (int i)
- {
- dat.set_size (i);
- }
-
- void insert (Real v, int i) {
- dat.insert (v,i);
- }
- void del (int i) { dat.del (i); }
-
- String str () const;
-
- void operator +=(Vector v) {
- assert (v.dim() == dim ());
- for (int i=0; i < dim(); i++)
- dat[i] += v.dat[i];
- }
-
- void operator -=(Vector v) {
- assert (v.dim() == dim ());
- for (int i=0; i < dim(); i++)
- dat[i] -= v (i);
- }
-
- Real operator *(Vector v) const {
- Real ip=0;
- assert (v.dim() == dim ());
- for (int i=0; i < dim(); i++)
- ip += dat[i] *v (i);
- return ip;
- }
- void operator *=(Real a) {
- for (int i=0; i < dim(); i++)
- dat[i] *= a;
- }
- void operator /=(Real a) {
- (*this) *= 1/a;
- }
- Vector operator-() const;
- Real norm_sq() {
- return ((*this) * (*this));
- }
- Real norm() {
- return sqrt (norm_sq());
- }
- operator Array<Real>() { return dat; }
- Array<Real> const &to_array()const { return dat; }
- void print() const;
- /// set to j-th element of unit-base
- void set_unit (int j) ;
-};
-
-inline Vector
-operator+(Vector a, Vector const &b) {
- a += b;
- return a;
-}
-
-inline Vector
-operator-(Vector a, Vector const &b) {
- a -= b;
- return a;
-}
-
-inline Vector
-operator*(Vector v, Real a) {
- v *= a;
- return v;
-}
-
-inline Vector
-operator*(Real a,Vector v) {
- v *= a;
- return v;
-}
-
-inline Vector
-operator/(Vector v,Real a) {
- v *= 1/a;
- return v;
-}
-
-#endif
+++ /dev/null
-/*
- matrix-debug.cc -- implement Matrix print routines
-
- source file of the Flower Library
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#include "flower-debug.hh"
-#include "matrix.hh"
-
-String
-Matrix::str () const
-{
- String s;
-#ifndef NPRINT
- Full_storage const * stor_c_l = dat_;
- s = String ("matrix {");
- for (int i=0; i< rows(); i++)
- {
- for (int j = 0; j < cols(); j++)
- {
- s+= to_str (stor_c_l->elem (i,j), "%6f ");
- }
- s+="\n";
- }
- s+="}\n";
-#endif
- return s;
-}
-
-
-void
-Matrix::print () const
-{
-#ifndef NPRINT
- fdebug << *this;
-#endif
-}
-
-String
-Vector::str () const
-{
- String s;
-#ifndef NPRINT
- s = String ("vector (") + to_str (dim ()) + ") [";
- for (int i=0; i < dim(); i++)
- {
- s += to_str (dat[i], "%6f") + to_str (' ');
- }
- s += "]\n";
-#endif
- return s;
-}
-
-
-void
-Vector::print() const
-{
-#ifndef NDEBUG
- fdebug << *this << '\n';
-#endif
-}
+++ /dev/null
-/*
- matrix.cc -- implement Matrix
-
- source file of the Flower Library
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-#include "matrix.hh"
-#include "full-storage.hh"
-
-Real
-Matrix::norm() const
-{
- Real r =0.0;
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- r += sqr (dat_->elem (i,j));
- return sqrt (r);
-}
-
-void
-Matrix::fill (Real r)
-{
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dat_->elem (i,j)=r;
-}
-
-void
-Matrix::set_diag (Real r)
-{
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dat_->elem (i,j)=(i==j) ? r: 0.0;
-}
-
-void
-Matrix::set_diag (Vector d)
-{
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dat_->elem (i,j)=(i==j) ? d (i): 0.0;
-}
-
-void
-Matrix::operator+=(Matrix const &m)
-{
- assert (m.cols() == cols ());
- assert (m.rows() == rows ());
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dat_->elem (i,j) += m (i,j);
-}
-
-void
-Matrix::operator-=(Matrix const &m)
-{
- assert (m.cols() == cols ());
- assert (m.rows() == rows ());
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dat_->elem (i,j) -= m (i,j);
-}
-
-
-void
-Matrix::operator*=(Real a)
-{
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dat_->elem (i,j) *= a;
-}
-
-void
-Matrix::operator=(Matrix const &m)
-{
- if (&m == this)
- return ;
- delete dat_;
- dat_ = new Full_storage (*m.dat_);
-}
-
-
-int
-Matrix::band_i () const
-{
- return dat_->band_i_;
-}
-
-void
-Matrix::set_band ()
-{
- dat_->band_i_ = calc_band_i ();
-}
-
-int
-Matrix::calc_band_i() const
-{
- int starty = dim();
- while (starty >= 0)
- {
- for (int i = starty, j = 0; i < dim(); i++, j++)
- if (dat_->elem (i,j))
- goto gotcha;
- for (int i=0, j = starty; j < dim(); i++,j++)
- if (dat_->elem (i,j))
- goto gotcha;
- starty --;
- }
- gotcha:
- return starty;
-}
-
-Matrix::Matrix (Matrix const &m)
-{
- dat_ = new Full_storage (*m.dat_);
-}
-
-
-Matrix::Matrix (int n, int m)
-{
- dat_ = new Full_storage(n,m);
- fill (0);
-}
-
-
-Matrix::Matrix (int n)
-{
- dat_= new Full_storage (n,n);
- fill (0);
-}
-
-Matrix::Matrix (Vector v, Vector w)
-{
- int n = v.dim();
- int m = w.dim ();
- dat_= new Full_storage (n,m);
- for (int i=0; i < n; i++)
- for (int j=0; j < m ; j++)
- dat_->elem (i,j)=v (i)*w (j);
-}
-
-
-Vector
-Matrix::row (int k) const
-{
- int n=cols();
-
-
- Vector v (n);
- for (int i=0; i < n; i++)
- v (i)=dat_->elem (k,i);
-
- return v;
-}
-
-Vector
-Matrix::col (int k) const
-{
- int n=rows();
- Vector v (n);
- for (int i=0; i < n; i++)
- v (i)=dat_->elem (i,k);
- return v;
-}
-
-Vector
-Matrix::left_multiply (Vector const & v) const
-{
- Vector dest (v.dim());
- assert (dat_->cols()==v.dim ());
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dest (i)+= dat_->elem (j,i)*v (j);
- return dest;
-}
-
-Vector
-Matrix::operator *(Vector const & v) const
-{
- Vector dest (rows());
- assert (dat_->cols()==v.dim ());
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- dest (i)+= dat_->elem (i,j)*v (j);
- return dest;
-}
-
-Matrix
-operator /(Matrix const& m1,Real a)
-{
- Matrix m (m1);
- m /= a;
- return m;
-}
-
-/*
- ugh. Only works for square matrices.
- */
-void
-Matrix::transpose()
-{
- for (int i=0, j=0; dat_->mult_ok (i,j); dat_->mult_next (i,j))
- {
- if (i >= j)
- continue;
- Real r=dat_->elem (i,j);
- dat_->elem (i,j) = dat_->elem (j,i);
- dat_->elem (j,i)=r;
- }
-}
-
-Matrix
-Matrix::operator-() const
-{
- OK();
- Matrix m (*this);
- m*=-1.0;
- return m;
-}
-
-Matrix
-Matrix::transposed() const
-{
- Matrix m (*this);
- m.transpose();
- return m;
-}
-
-Matrix
-operator *(Matrix const &m1, Matrix const &m2)
-{
- Matrix result (m1.rows (), m2.cols ());
- result.set_product (m1,m2);
- return result;
-}
-
-
-
-void
-Matrix::set_product (Matrix const &m1, Matrix const &m2)
-{
- assert (m1.cols()==m2.rows ());
- assert (cols()==m2.cols () && rows ()==m1.rows ());
-
- for (int i=0; i < rows (); i++)
- for (int j=0; j < cols (); j++)
- {
- Real r=0.0;
- for (int k = 0 >? i - m1.band_i () >? j - m2.band_i ();
- k < m1.cols(); k++)
- {
- r += m1(i,k)*m2(k,j);
- }
- dat_->elem (i,j)=r;
- }
-}
-
-void
-Matrix::insert_row (Vector v, int k)
-{
- int c = cols();
- assert (v.dim()==cols ());
- dat_->insert_row (k);
- for (int j=0; j < c; j++)
- dat_->elem (k,j)=v (j);
-}
-
-
-void
-Matrix::swap_columns (int c1, int c2)
-{
- assert (c1>=0&& c1 < cols()&&c2 < cols () && c2 >=0);
- int r = rows();
- for (int i=0; i< r; i++)
- {
- Real r=dat_->elem (i,c1);
- dat_->elem (i,c1) = dat_->elem (i,c2);
- dat_->elem (i,c2)=r;
- }
-}
-
-void
-Matrix::swap_rows (int c1, int c2)
-{
- assert (c1>=0&& c1 < rows()&&c2 < rows () && c2 >=0);
- int c = cols();
- for (int i=0; i< c; i++)
- {
- Real r=dat_->elem (c1,i);
- dat_->elem (c1,i) = dat_->elem (c2,i);
- dat_->elem (c2,i)=r;
- }
-}
-
-
-int
-Matrix::dim() const
-{
- assert (cols() == rows ());
- return rows();
-}
-
-Matrix::~Matrix ()
-{
- delete dat_;
-}
+++ /dev/null
-/*
- vector.cc -- implement Vector
-
- source file of the Flower Library
-
- (c) 1996-98 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-
- */
-
-#include "vector.hh"
-
-Vector::Vector (Array<Real> d)
- : dat (d)
-{
-
-}
-Vector::Vector (Vector const &n)
- : dat (n.dat)
-{
-}
-
-Vector
-Vector::operator-() const
-{
- Vector v (*this);
- v*=-1;
- return v;
-}
-
-void
-Vector::set_unit (int j)
-{
- fill (0.0);
- dat[j] = 1.0;
-}
+++ /dev/null
-\score{
- \addlyrics
- \context Staff {
- \notes\relative c''{ c () d r e f }
- \notes\relative c''{ c ~ c r c c }
- }
- {
- \context Lyrics \lyrics { foo __ bar baz }
- \context Lyrics \lyrics { foo __ bar baz }
- }
- \paper{
- \translator{
- \VoiceContext
- automaticMelismas=1;
- }
- }
-}
+++ /dev/null
-
-
-The following example shows 2 bugs with HaraKiriStaffContext.
-
-The first bug is that when the second staff disappears, the instr
-string in the margin "2" is still printed, and surprisingly appears
-above staff 1.
-
-The second and more serious bug is that when staff 2 reappears,
-the bar separators are no longer being printed for that staff.
-
-Regards,
-Roy Rankin
-
-\header {
- composer = "Music: Trad.";
- crossRefNumber = "4";
- title = "HaraKiri Bugs";
-}
-voice1 = \notes {
-\property Staff.instrument = "Part1"
-\property Staff.instr = "1"
-\property Staff.timeSignatureStyle="C"
- \time 4/4; \key c;
- e''4. e''8 e''4 r4 | e''4 d''4
- c''2 |
- e''4 d''4 c''2 ( |
- e''2 c''2 ~ |
- ) c''1 \bar "|."; \break
- e''4. e''8 e''4 r4 | e''4 d''4
- c''2 |
- e''4 d''4 c''2 ( |
- e''2 c''2 ~ |
- ) c''1 \bar "|."; \break
- e''4. e''8 e''4 r4 | e''4 d''4
- c''2 |
- e''4 d''4 c''2 ( |
- e''2 c''2 ~ |
- ) c''1 \bar "|."; \break
- e''4. e''8 e''4 r4 | e''4 d''4
- c''2 |
- e''4 d''4 c''2 ( |
- e''2 c''2 ~ |
- ) c''1 \bar "|.";
-}
-voice2 = \notes {
-\property Staff.instrument = "Part2"
-\property Staff.instr = "2"
-\property Staff.timeSignatureStyle="C"
- \time 4/4; \key c;
- e''4. e''8 e''4 r4 | e''4 d''4
- c''2 |
- e''4 d''4 c''2 ( |
- e''2 c''2 ~ |
- ) c''1 \bar "|.";
-s1 * 10
- e''4. e''8 e''4 r4 | e''4 d''4
- c''2 |
- e''4 d''4 c''2 ( |
- e''2 c''2 ~ |
- ) c''1 \bar "|.";
-}
-\score{
- \notes <
-
- \context Staff="1"
- {
- \$voice1
- }
- \context Staff="2"
- {
- \$voice2
- }
-
- >
- \paper
- {
- \translator
- {
- \HaraKiriStaffContext
- \consists Staff_margin_engraver;
- }
- }
- \midi {}
-}
-
-
-
+++ /dev/null
-\score {
- \notes {
- \repeat volta 5 {
- c' c' c' c'
- }
- \alternative {
- {d' d' d' d'}
- {e' e' e' e'}
- }
- f'
- }
- }
-
+++ /dev/null
-\header {\r
-filename = "utremi.ly";\r
-enteredby = "Christian Mondrup";\r
-%composer = "Thomas Ravenscroft";\r
-opus = "Thomas Ravenscroft";\r
-%opus = "Pammelia, 1609, no. 31";\r
-arranger = "Pammelia, 1609, no. 31";\r
-title = "Vt, re, me, fa, sol, la, sol, fa, re, vt";\r
-tagline = "Typeset with GNU LilyPond by Christian Mondrup. Non-commercial copying welcome.";\r
-}\r
-\r
-papersize = "letter"\r
-\r
-\include "paper16.ly"\r
-%\include "a4.ly"\r
-\r
-lyricsOne = \lyrics {\r
- ""\breve ""8\r
- "Vt,"\breve re, me, fa, sol, la, ""\r
- sol, la, fa, me, re, vt. \r
-}\r
-\r
-lyricsTwo = \lyrics {\r
- "" ""\r
- Hey downe downe hey downe downe downe downe hey down\r
- hey down "" down down a.\r
- My heart of gold as true as steele\r
- as I me leant "" vn -- to the bowres.\r
- but if my La -- dy loue me well,\r
- Lord so Ro -- bin bowres,\r
-}\r
-\r
-lyricsThree = \lyrics {\r
- "" "" "" "" "" "" "" "" "" "" "" ""\r
- heaue and hoe Rum -- be -- lo, hey tro -- lo tro -- ly lo,\r
- hey tro -- ly trol -- ly hey tro -- ly trol -- ly \r
- hey tro -- ly trol -- ly hey trol -- "" ly trol -- ly lo. \r
- My La -- dies gone to Can -- ter -- bu -- ry\r
- S. Tho -- mas be her boote.\r
- Shee met with Kate of Malms -- "" bu -- ry, \r
- why weepst thou ma -- ple \r
-}\r
-\r
-lyricsFour = \lyrics {\r
- "" \r
- roote:\r
- O sleepst thou or wakst thon Ies -- se -- ry, Cooke,\r
- the rost it burnes, turne round turne round "" a -- bout ""\r
- turne round "" a -- bout, ""\r
- turne round "" a -- bout, "" turne round. ""\r
- O Fri -- er how fares thy ban -- de -- low ban -- de -- low\r
- Fri -- er, how fares thy San -- de -- low, San -- de -- low.\r
-}\r
-\r
-global = \notes {\r
- \key c;\r
- %\time 4/2;\r
- %\time 17/8;\r
- \time 33/16;\r
- \clef "tenor";\r
-}\r
-\r
-dummyBeat = \notes {\r
- \property Staff.defaultBarType = "||"\r
- \property Staff.barAlways = "1" \r
- s16\r
- \property Staff.barAlways = "0"\r
- \property Staff.defaultBarType = "|"\r
-}\r
-\r
-incipOne = \notes\relative c' {\r
- g\breve \r
- \dummyBeat\r
- \bar "";\r
-}\r
-\r
-incipTwo = \notes\relative c'' {\r
- \property Voice.noteHeadStyle = "diamond"\r
- g1. \r
- \property Voice.noteHeadStyle = ""\r
- f2 \r
- \dummyBeat\r
- \bar "";\r
-}\r
-\r
-incipThree = \notes\relative c' {\r
- \property Voice.noteHeadStyle = "harmonic"\r
- \property Voice.tupletVisibility = 0\r
- \times 2/3 { d4. e8 d4 } \times 2/3 { b4. a8 g4 }\r
- \times 2/3 { b4. c8 b4 } \times 2/3 { b4. a8 g4 }\r
- \dummyBeat\r
- \property Voice.tupletVisibility = 3\r
- \bar "";\r
-}\r
-\r
-incipFour = \notes\relative c' {\r
- \property Voice.noteHeadStyle = "diamond"\r
- %\property Voice.restStyle = "mensural"\r
- b1 r1\r
- \dummyBeat\r
- \bar "";\r
-}\r
-\r
-partOne = \notes\relative c' {\r
- %\property Score.currentBarNumber = "1"\r
- \property Voice.noteHeadStyle = ""\r
- %\property Staff.barNumberScriptPadding = 1\r
- g1 ~ | g | a ~ | a | bes ~ | bes | c ~ | c | d ~ | d | e ~ | e |\r
- r | r | \r
- e ~ | e | d ~ | d | c ~ | c | bes ~ | bes | a ~ | a | g ~ | g\r
- \bar "|.";\r
-}\r
-\r
-partTwo = \notes\relative c'' {\r
- \property Voice.noteHeadStyle = ""\r
- g2. f4 | e2 d | c4. b8 a4 g | f2 e |\r
- d g ~ | g4 a g2 | c,1 | r2 c |\r
- g'2. a4 | b2 g | c c | c c |\r
- b4 c d2 ~ | d4 d b g | c1 | r2 c |\r
- g g | d2. e4 | f2 e | c1 |\r
- g'2. g4 | g2 g | c,1 ~ | c | r | r \r
- \bar "|.";\r
-}\r
-\r
-partThree = \notes\relative c' {\r
- \property Voice.noteHeadStyle = ""\r
- \times 2/3 { d4. e8 d4 } \times 2/3 {b4. a8 g4 } | \r
- \times 2/3 { b4. c8 b4 } \times 2/3 { b4. a8 g4 }\r
- c2 c4 c | c2 c |\r
- d2 d4 d | d2 d | e2. f4 | e2 d4 c |\r
- b2 g ~ | g4 f e d | c1 | r2 c |\r
- g'2. g4 | g2 g | c,2. c4 | c c c'2 |\r
- bes g | bes bes | a1 | r2 e' |\r
- d d | g g | a2. g4 | fis e d2 |\r
- g g | d2. c4 \r
- \bar "|.";\r
-}\r
-\r
-partFour = \notes\relative c' {\r
- \property Voice.noteHeadStyle = ""\r
- b1^\fermata | r2 g' | e e4 e | a2. g4 | \r
- fis2 e4 d | g1 | g2 g | g g |\r
- g g | g g ~ | g4 f e d | c2 d ~ |\r
- d4 c b a | g2 g' ~ | g4 f e d | c2 g' ~ |\r
- g r | \times 2/3 { r1 d2 } \times 2/3 { a'2 a g } |\r
- \times 2/3 { e1 f2 } \times 2/3 { g2. g4 g2 }\r
- \times 2/3 { d2. d4 e2 } \times 2/3 { f2. f4 e2 }\r
- \times 2/3 { d1 c2 } \times 2/3 { b2. a4 g2 }\r
- \times 2/3 { d'2. d4 d2 }\r
- \bar "|.";\r
-}\r
-\r
-partOneStaff = <\r
- \context Staff = vocal1 <\r
- \property Staff.clefStyle = "fullSizeChanges"\r
- \property Staff.instrument = "\large{1.}"\r
- %\property Staff.instr = ""\r
- \notes { \r
- \global\r
- \property Staff.timeSignatureStyle = "old4/4"\r
- \incipOne\r
- \clef "G_8";\r
- \property Staff.timeSignatureStyle = "C2/2"\r
- \time 4/4;\r
- \partOne \r
- }\r
- {\context Lyrics = lyrOne \lyricsOne }\r
- >\r
->\r
-\r
-partTwoStaff = <\r
- \context Staff = vocal2 <\r
- \property Staff.clefStyle = "fullSizeChanges"\r
- \property Staff.instrument = "\large{2.}"\r
- %\property Staff.instr = ""\r
- %\property Voice.automaticMelismata = "1"\r
- \addlyrics\r
- \notes { \r
- \global\r
- \property Staff.timeSignatureStyle = "old4/4"\r
- \incipTwo\r
- \clef "G_8";\r
- \property Staff.timeSignatureStyle = "C2/2"\r
- \time 4/4;\r
- \partTwo\r
- }\r
- {\context Lyrics = lyrTwo \lyricsTwo }\r
- >\r
->\r
-\r
-\r
-partThreeStaff = <\r
- \context Staff = vocal3 <\r
- \property Staff.clefStyle = "fullSizeChanges"\r
- \property Staff.instrument = "\large{3.}"\r
- %\property Staff.instr = ""\r
- %\property Voice.automaticMelismata = "1"\r
- \addlyrics\r
- \notes { \r
- \context Voice = vthree \r
- \property Voice.tupletDirection = \up\r
- \global\r
- \property Staff.timeSignatureStyle = "old4/4"\r
- \incipThree\r
- \clef "G_8";\r
- \property Staff.timeSignatureStyle = "C2/2"\r
- \time 4/4;\r
- \partThree\r
- }\r
- {\context Lyrics = lyrThree \lyricsThree }\r
- >\r
->\r
-\r
-partFourStaff = <\r
- \context Staff = vocal4 <\r
- \property Staff.clefStyle = "fullSizeChanges"\r
- \property Staff.timeSignatureStyle = "C2/2"\r
- \property Staff.barNumberDirection = \up\r
- \property Staff.instrument = "\large{4.}"\r
- \property Staff.instr = ""\r
- %\property Voice.automaticMelismata = "1"\r
- \addlyrics\r
- \notes { \r
- \context Voice = vfour \r
- \property Voice.tupletDirection = \up\r
- \global\r
- %\property Staff.instrument = "\large{4.}"\r
- %\property Staff.instr = ""\r
- \property Staff.timeSignatureStyle = "old4/4"\r
- \incipFour\r
- \clef "G_8";\r
- \property Staff.timeSignatureStyle = "C2/2"\r
- \time 4/4;\r
- \partFour \r
- }\r
- {\context Lyrics = lyrFour \lyricsFour} \r
- >\r
->\r
-\r
-\score {\r
- \context StaffGroup <\r
- \partFourStaff\r
- \partThreeStaff\r
- \partTwoStaff\r
- \partOneStaff\r
- >\r
- %\paper {\translator {\BarNumberingStaffContext } } \r
- \paper {\r
- indent = 0.\mm;\r
- linewidth = 17.0\cm;\r
- textheight = 27.0\cm;\r
- gourlay_maxmeasures=6.0;\r
- \translator { \StaffContext\r
- \consists "Staff_margin_engraver"; } \r
- %\translator { \ScoreContext\r
- % minVerticalAlign = 1.5*\staffheight; }\r
- }\r
- %\midi { \r
- % output = "utremi.mid";\r
- % \tempo 2 = 80;\r
- %}\r
-}\r
+++ /dev/null
-/*
- col-info.cc -- implement Column_info
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-#include "paper-column.hh"
-#include "col-info.hh"
-#include "debug.hh"
-
-void
-Column_info::print() const
-{
-#ifndef NPRINT
- DOUT << "column { ";
- if (fixed_b())
- DOUT << "fixed at " << fixed_position() << ", ";
- assert (pcol_l_);
- DOUT << width_.str();
- Direction d = LEFT;
- do {
- for (int i=0; i < rods_[d].size (); i++)
- rods_[d][i].print ();
- } while (flip (&d) != LEFT);
-
- DOUT <<"}\n";
-#endif
-}
-
-Column_info::Column_info (Paper_column *col_l, Real const *fixed_C)
-{
- if (fixed_C)
- fixpos_p_.set_l (fixed_C);
- pcol_l_ = col_l;
- width_ = pcol_l_->extent(X_AXIS);
- if (width_.empty_b())
- width_ = Interval(0,0);
-}
-
-
-Column_info::Column_info()
-{
- pcol_l_ =0;
-}
-
-bool
-Column_info::fixed_b () const
-{
- return fixpos_p_.get_C();
-}
-
-Real
-Column_info::fixed_position () const
-{
- return *fixpos_p_;
-}
-
-int
-Column_info::rank_i () const
-{
- return pcol_l_->rank_i ();
-}
-
-void
-Spacer_rod::print ()const
-{
- DOUT << "Other " << other_idx_ << "dist = " << distance_f_ << '\n';
-}
+++ /dev/null
-/*
- ineq-constrained-qp.hh -- declare Ineq_constrained_qp
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#ifndef INEQ_CONSTRAINED_QP_HH
-#define INEQ_CONSTRAINED_QP_HH
-
-
-
-#include "matrix.hh"
-
-/** inequality constrained quadratic program
-
- It takes the form of
-
- optimise for x : x*quad_ *x + lin_* x + const_term_
-
- subject to for all i: cons_[i] * x >= consrhs_[i]
-
-
- @usage:
- instantiate Ineq_constrained_qp.
-
- Modify quad_, lin_ and const_term_ directly. Use
- add_inequality_cons () to add constraints. Call solve () with a
- feasible solution
-
-
- */
-class Ineq_constrained_qp {
- friend class Active_constraints;
-
- Array<Vector> cons_;
- Array<Real> consrhs_;
-public:
- Matrix quad_;
- Vector lin_;
- Real const_term_;
-
-
- /**
- use a KKT method to assert optimality of sol
- */
- void assert_solution (Vector sol) const;
- /// solve the problem using a projected gradient method
- Vector constraint_solve (Vector) const;
- /**
- Solve it. First try it the easy way.
- */
- Vector solve (Vector start) const;
-
- /**
- @return the number of variables in the problem
- */
- int dim() const;
-
- /**
- add a constraint
-
-
- c*vars >= r
-
- PRE
- c.dim() == dim ();
-
- */
- void add_inequality_cons (Vector c, double r);
-
- /** set up matrices to go with the problem. */
- Ineq_constrained_qp (int novars);
-
- /**
- evaluate the quadratic function for input #v#
- */
- Real eval (Vector v);
-
- void eliminate_var (int idx, Real value);
- void OK() const;
- void print() const;
-
-};
-
-// ugh
-const Real EPS=1e-7;
-
-#endif // INEQ_CONSTRAINED_QP_HH
+++ /dev/null
-/*
- qlp.hh -- declare Ineq_constrained_qp, Mixed_qp
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-#ifndef QLP_HH
-#define QLP_HH
-
-#include "ineq-constrained-qp.hh"
-
-/**
- Quadratic programming with mixed linear constraints.
- problem definition of a quadratic optimisation problem with linear
- inequality and equality constraints
-
-
- x^T QUAD x /2 + b^T x
-*/
-class Mixed_qp :public Ineq_constrained_qp {
- Array<int> eq_cons;
- Array<Real> eq_consrhs;
-public:
- Mixed_qp (int n);
- void OK() const;
- void print() const;
-
- Vector solve (Vector start) const;
- void add_fixed_var (int i , Real value);
-
-
- /**
- add a constraint,
-
- c*vars == r
-
- PRE
- c.dim()==dim ();
- */
- void add_equality_cons (Vector c, double r);
-};
-#endif
+++ /dev/null
-/*
- qlpsolve.hh -- declare Active_constraints, Inactive_iter
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#ifndef QLPSOLVE_HH
-#define QLPSOLVE_HH
-
-#include "matrix.hh"
-
-
-/**
- This class represents the set of active (binding) constraints
- which can be active while the QLP algorithm is in a feasible
- point. The active constraints are numbered.
- If the constraints are of the form
-
- A^T*x >= b
-
- then the binding constraints are those where the >= is equality.
-
- */
-
-class Active_constraints {
- friend class Inactive_iter;
-
-
- Matrix A,H;
- Array<int> active;
- Array<int> inactive; // actually this is a set, not an array.
- Ineq_constrained_qp const *opt;
-
-public:
-
- /** This counts the number of errors the algorithms makes. The
- optimum search should be abandoned if it becomes too high. */
- int degenerate_count_i_;
- String status() const;
-
- Vector vec (int k) const { return opt->cons_[k]; }
- Real rhs (int k) const { return opt->consrhs_[k]; }
-
-
- /** drop constraint. drop constraint k from the active set. k is the index of the
- constraint in #active#
-
- */
- void drop_constraint (int k);
-
-
- /** add constraint j.
- add constraint j to the active set j is the index of the
- constraint in #inactive#
- */
- void add_constraint (int j);
-
- /// exchange in and out.
- void exchange (int in, int out) { add_constraint (in); drop_constraint (out); }
-
-
- Vector find_active_optimum (Vector g);
-
- /// get lagrange multipliers.
- Vector get_lagrange (Vector v);
-
- Active_constraints (Ineq_constrained_qp const *op);
- /** construct: no constraints active, n vars. Put the equalities
- into the constraints. */
-
- /// check invariants
- void OK();
-};
-
-
-/**
- loop through the inactive constraints.
- */
-class Inactive_iter {
- int j;
- Active_constraints const* ac;
-public:
- Inactive_iter (Active_constraints const &c) { ac=&c; j=0; }
- int idx() const { return j; }
- void operator ++(int) { j++; }
- int constraint_id() const { return ac->inactive[j]; }
- Vector vec() const { return ac->vec (constraint_id ()); }
- Real rhs() const { return ac->rhs (constraint_id ()); }
- bool ok() const { return j < ac->inactive.size (); }
-};
-
-#endif // QLPSOLVE_HH
+++ /dev/null
-/*
- spring-spacer.hh -- declare Spring_spacer
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#ifndef SPRING_SPACER_HH
-#define SPRING_SPACER_HH
-
-#include "line-spacer.hh"
-#include "cons.hh"
-#include "col-info.hh"
-#include "column-x-positions.hh"
-#include "moment.hh"
-
-
-
-/**
- Determine positions of columns connected by springs and held apart by rods.
-
-
- Generate a spacing which
- \begin{itemize}
- \item
- Satisfies spacing constraints (notes can't be printed through each other)
- \item
- Looks good, ie tries to conform to an ideal spacing as much as possible.
- \end{itemize}
-
- This is converted by regarding idealspacing as "springs" attached
- to columns. The equilibrium of one spring is the ideal
- distance. The columns have a size, this imposes "hard" constraints
- on the distances. This transforms the problem into a quadratic
- programming problem with linear constraints.
-
- The quality is given by the total potential energy in the
- springs. The lower the energy, the better the configuration.
-
-
- TODO: too complicated. Revise.
- Use force iso. energy. Also optimise for uniform density
-*/
-
-class Spring_spacer : public Line_spacer {
- friend class Durations_iter;
-private:
- // can't copy me.
- Spring_spacer (Spring_spacer const&s);
- Cons<Idealspacing> *ideal_p_list_;
- Array<Column_info> cols_;
- Real indent_f_;
-
- /// the index of #c# in #cols#
- int col_id (Paper_column const *c) const;
-
- /// generate an (nonoptimal) solution
- Vector find_initial_solution() const;
-
- /// does #this# contain the column #w#?
- bool contains_b (Paper_column const *w);
-
- /// make the energy function
- void make_matrices (Matrix &quad, Vector &lin,Real&) const;
-
-
- /// generate the LP constraints
- void make_constraints (Mixed_qp& lp) const;
-
-
- void handle_loose_cols();
- bool try_initial_solution_and_tell (Vector&)const;
- Vector try_initial_solution() const;
-
- void set_fixed_cols (Mixed_qp&) const;
-
- Score_column* scol_l (int);
- void connect (int i,int j, Real,Real);
- Real calculate_energy_f (Vector) const;
-public:
- static Line_spacer *constructor();
- Real energy_normalisation_f_;
- Spring_spacer ();
- virtual ~Spring_spacer ();
- virtual void solve (Column_x_positions*) const;
- virtual void lower_bound_solution (Column_x_positions*) const;
- virtual void add_columns (Link_array<Paper_column>);
- void add_column (Paper_column *, bool, Real);
-
- virtual Vector default_solution() const;
- virtual bool check_constraints (Vector v) const;
- virtual void OK() const;
- virtual void print() const;
- virtual void prepare();
-};
-
-#endif // SPRING_SPACER_HH
+++ /dev/null
-/*
- ineq-constrained-qp.cc -- implement Ineq_constrained_qp
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1996, 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-#include "ineq-constrained-qp.hh"
-#include "qlpsolve.hh"
-#include "debug.hh"
-#include "choleski.hh"
-#include "main.hh"
-
-/*
- May be this also should go into a library
- */
-
-const int MAXITER=100; // qlpsolve.hh
-
-const int MAXDEGEN=5;
-
-/*
- 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=infinity_f;
- int idx=-1;
- for (int i = 0; i < v.dim(); i++)
- {
- if (v (i) < m)
- {
- idx = i;
- m = v (i);
- }
- assert (v (i) <= infinity_f);
- }
- 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::constraint_solve (Vector start) const
-{
- if (!dim())
- return Vector (0);
-
- 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 && act.degenerate_count_i_ < MAXDEGEN)
- {
- if (experimental_features_global_b)
- assert_solution (x);
-
- Vector direction= - act.find_active_optimum (gradient);
-
- DOUT << "gradient "<< gradient<< "\ndirection " << direction<< '\n';
-
- if (direction.norm() > EPS)
- {
- DOUT << act.status() << '\n';
-
- Real unbounded_alfa = 1.0;
- Real minalf = unbounded_alfa;
-
- 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++)
- {
- Real dot = ia.vec() * direction;
- Real mindot = (experimental_features_global_b)
- ? -EPS
- : 0;
-
- if (dot >= mindot)
- continue;
-
-
- Real numerator = ia.rhs () - ia.vec()*x;
- if (numerator >= 0)
- {
- if (numerator > EPS)
- {
- warning (_f ("Ineq_constrained_qp::solve (): Constraint off by %f", numerator));
- act.degenerate_count_i_ ++;
- }
- minalf = -numerator;
- minidx = ia;
- break;
- }
-
- Real alfa = numerator / dot;
-
-
- if (minalf > alfa)
- {
- minidx = ia;
- minalf = alfa;
- }
- }
-
- Real optimal_step = minalf;
-
- Vector deltax = direction * optimal_step;
- x += deltax;
- gradient += optimal_step * (quad_ * deltax);
-
- DOUT << "step = " << optimal_step << " (|dx| = " <<
- to_str (deltax.norm()) << ")\n";
-
- if (minalf < unbounded_alfa)
- {
- /* bumped into an edge. try again, in smaller space. */
- act.add_constraint (minidx.idx());
- DOUT << "adding cons "<< minidx.idx () << '\n';
- continue;
- }
- /*ASSERT: we are at the 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)
- {
- Real n =gradient.norm();
- if (n >= EPS)
- {
- programming_error ("Huh? Gradient should be zero ");
- act.degenerate_count_i_ ++;
- }
- else if (n <0)
- {
- programming_error ("Huh? Norm should be positive");
- act.degenerate_count_i_++;
- }
- break;
- }
-
- DOUT << "dropping cons " << m << '\n';
- act.drop_constraint (m);
- }
- if (iterations >= MAXITER)
- WARN << _ ("Didn't converge!") << '\n';
- if (act.degenerate_count_i_ >= MAXDEGEN)
- WARN << _ ("too much degeneracy") << '\n';
- DOUT << ": found " << x << " in " << iterations <<" iterations\n";
- assert_solution (x);
- return x;
-}
-
-
-Vector
-Ineq_constrained_qp::solve (Vector start) const
-{
- /* no hassle if no constraints*/
- if (! cons_.size())
- {
- Choleski_decomposition chol (quad_);
- return - chol.solve (lin_);
- }
- else
- {
- return constraint_solve (start);
- }
-}
-
-int
-Ineq_constrained_qp::dim () const
-{
- return lin_.dim();
-}
-
-
-
-
-void
-Ineq_constrained_qp::assert_solution (Vector sol) const
-{
- bool sol_b=true;
-
- for (int i=0; sol_b && i < cons_.size(); i++)
- {
- Real R=cons_[i] * sol- consrhs_[i];
- if (R> -EPS)
- sol_b = false;
- }
-}
-
-void
-Ineq_constrained_qp::print() const
-{
-#ifndef NPRINT
- DOUT << "Ineq_constrained_qp { " << '\n';
- DOUT << "Quad " << quad_;
- DOUT << "lin " << lin_ << '\n'
- << "const " << const_term_<< '\n';
- for (int i=0; i < cons_.size(); i++)
- {
- DOUT << "constraint["<<i<<"]: " << cons_[i] << " >= " << consrhs_[i];
- DOUT << '\n';
- }
- DOUT << "}\n";
-#endif
-}
+++ /dev/null
-/*
- qlp.cc -- implement Mixed_qp
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-#include "debug.hh"
-#include "qlp.hh"
-
-
-void
-Mixed_qp::add_equality_cons (Vector , double)
-{
- assert (false);
-}
-
-void
-Mixed_qp::add_fixed_var (int i, Real r)
-{
- eq_cons.push (i);
- eq_consrhs.push (r);
-}
-
-
-/**
- eliminate appropriate variables, until we have a Ineq_constrained_qp
- then solve that.
-
- PRE
- cons should be ascending
- */
-Vector
-Mixed_qp::solve (Vector start) const
-{
- if (!dim())
- return Vector (0);
-
- print();
- Ineq_constrained_qp pure (*this);
-
- for (int i= eq_cons.size()-1; i>=0; i--)
- {
- pure.eliminate_var (eq_cons[i], eq_consrhs[i]);
- start.del (eq_cons[i]);
- }
- Vector sol = pure.solve (start);
- for (int i= 0; i < eq_cons.size(); i++)
- {
- sol.insert (eq_consrhs[i],eq_cons[i]);
- }
- return sol;
-}
-
-
-Mixed_qp::Mixed_qp (int n)
- : Ineq_constrained_qp (n)
-{
-}
-
-void
-Mixed_qp::OK() const
-{
-#ifndef NDEBUG
- Ineq_constrained_qp::OK();
- assert (eq_consrhs.size() == eq_cons.size ());
-#endif
-}
-
-void
-Mixed_qp::print() const
-{
-#ifndef NPRINT
- Ineq_constrained_qp::print();
- for (int i=0; i < eq_cons.size(); i++)
- {
- DOUT << "eq cons "<<i<<": x["<<eq_cons[i]<<"] == " << eq_consrhs[i]<< '\n';
- }
-#endif
-}
-
+++ /dev/null
-/*
- paper-column.cc -- implement Paper_column
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-#include "paper-column.hh"
-#include "paper-score.hh"
-#include "debug.hh"
-
-void
-Paper_column::add_rod (Paper_column * p, Real d)
-{
- Direction dir = Direction (sign (p->rank_i () - rank_i ()));
-
- if (!dir)
- {
- programming_error ("Must set minimum distance between differing columns.");
- return;
- }
-
- for (int i=0; i < minimal_dists_arr_drul_[dir].size (); i++)
- {
- Column_rod &rod = minimal_dists_arr_drul_[dir][i];
- if (rod.other_l_ == p)
- {
- rod.distance_f_ = rod.distance_f_ >? d;
- return ;
- }
- }
-
- Column_rod cr;
- cr.distance_f_ = d;
- cr.other_l_ = p;
-
- minimal_dists_arr_drul_[dir].push (cr);
-}
-
-void
-Paper_column::add_spring (Paper_column * p, Real d, Real s)
-{
- Direction dir = Direction (sign (p->rank_i () - rank_i ()));
-
- if (!dir)
- {
- warning (_ ("Must set spring between differing columns"));
- return;
- }
-
- for (int i=0; i < spring_arr_drul_[dir].size (); i++)
- {
- Column_spring &spring = spring_arr_drul_[dir][i];
- if (spring.other_l_ == p)
- {
- spring.distance_f_ = spring.distance_f_ >? d;
- return ;
- }
- }
-
- Column_spring cr;
- cr.distance_f_ = d;
- cr.strength_f_ = s;
- cr.other_l_ = p;
-
- spring_arr_drul_[dir].push (cr);
-}
-
-int
-Paper_column::rank_i() const
-{
- return rank_i_;
-}
-
-void
-Paper_column::set_rank (int i)
-{
- rank_i_ = i;
-}
-
-void
-Paper_column::do_print() const
-{
-#ifndef NPRINT
- DOUT << "rank: " << rank_i_ << '\n';
- Direction d = LEFT;
- do
- {
- for (int i=0; i < minimal_dists_arr_drul_[d].size (); i++)
- {
- minimal_dists_arr_drul_[d][i].print ();
- }
- for (int i=0; i < spring_arr_drul_[d].size (); i++)
- {
- spring_arr_drul_[d][i].print ();
- }
-
- }
- while ((flip (&d))!=LEFT);
- Item::do_print ();
-#endif
-}
-
-bool
-Paper_column::breakpoint_b() const
-{
- return !line_l_;
-}
-
-Paper_column::Paper_column()
-{
- set_axes (X_AXIS, X_AXIS);
-
- line_l_=0;
- rank_i_ = -1;
-}
-
-Line_of_score*
-Paper_column::line_l() const
-{
- return line_l_;
-}
-
-
-
-
-Paper_column*
-Paper_column::column_l () const
-{
- return (Paper_column*)(this);
-}
-
-/*
- ugh.
- */
-void
-Paper_column::preprocess ()
-{
- minimal_dists_arr_drul_[LEFT].sort (Column_rod::compare);
- minimal_dists_arr_drul_[RIGHT].sort (Column_rod::compare);
- spring_arr_drul_[LEFT].sort (Column_spring::compare);
- spring_arr_drul_[RIGHT].sort (Column_spring::compare);
-}
+++ /dev/null
-/*
- qlpsolve.cc -- implement Active_constraints, Inactive_iter
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1996, 1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-
- TODO:
- try fixed point arithmetic, to speed up lily.
- */
-
-#include "ineq-constrained-qp.hh"
-#include "qlpsolve.hh"
-#include "debug.hh"
-#include "choleski.hh"
-
-const Real TOL=1e-1; // roughly 1/30 mm
-
-String
-Active_constraints::status() const
-{
- String s ("Active|Inactive [");
- for (int i=0; i< active.size(); i++)
- {
- s += to_str (active[i]) + " ";
- }
-
- s+="| ";
- for (int i=0; i< inactive.size(); i++)
- {
- s += to_str (inactive[i]) + " ";
- }
- s+="]";
-
- return s;
-}
-
-void
-Active_constraints::OK()
-{
-#ifndef NDEBUG
- H.OK();
- A.OK();
- assert (active.size() +inactive.size () == opt->cons_.size ());
- assert (H.dim() == opt->dim ());
- assert (active.size() == A.rows ());
- Array<int> allcons;
-
- for (int i=0; i < opt->cons_.size(); i++)
- allcons.push (0);
- for (int i=0; i < active.size(); i++)
- {
- int j = active[i];
- allcons[j]++;
- }
- for (int i=0; i < inactive.size(); i++)
- {
- int j = inactive[i];
- allcons[j]++;
- }
- for (int i=0; i < allcons.size(); i++)
- assert (allcons[i] == 1);
-#endif
-}
-
-Vector
-Active_constraints::get_lagrange (Vector gradient)
-{
- return (A*gradient);
-}
-
-void
-Active_constraints::add_constraint (int k)
-{
- // add indices
- int cidx=inactive[k];
-
- Vector a (opt->cons_[cidx]);
- // update of matrices
- Vector Ha = H*a;
- Real aHa = a*Ha;
- Vector addrow (Ha.dim());
- bool degenerate = (abs (aHa) < EPS);
-
- if (degenerate)
- {
- warning (String ("Active_constraints::add ():")
- + _("degenerate constraints"));
- DOUT << "Ha = " << Ha.str () << '\n';
- /*
- a != 0, so if Ha = O(EPS), then
- Ha * aH / aHa = O(EPS^2/EPS)
-
- if H*a == 0, the constraints are dependent.
- */
- degenerate_count_i_ ++;
- }
- if (!degenerate)
- {
- active.push (cidx);
- inactive.swap (k,inactive.size()-1);
- inactive.pop();
-
- H -= Matrix (Ha/aHa , Ha);
-
- addrow=Ha;
- addrow /= aHa;
- A -= Matrix (A*a, addrow);
- A.insert_row (addrow,A.rows());
- }
-}
-
-void
-Active_constraints::drop_constraint (int k)
-{
- int q=active.size()-1;
-
-
- Vector a (A.row (q));
- if (a.norm() > EPS)
- {
- // drop indices
- inactive.push (active[k]);
- active.swap (k,q);
- A.swap_rows (k,q);
- active.pop();
- /*
-
- */
- Real aqa = a*opt->quad_*a;
- Matrix aaq (a,a/aqa);
- H += aaq;
- A -= A*opt->quad_*aaq;
-
- A.delete_row (q);
- }else {
- degenerate_count_i_ ++;
- warning (String ("Active_constraints::drop ():")
- + _("degenerate constraints"));
- }
-}
-
-
-Active_constraints::Active_constraints (Ineq_constrained_qp const *op)
- : A(0,op->dim()),
- H(op->dim()),
- opt (op)
-{
- for (int i=0; i < op->cons_.size(); i++)
- inactive.push (i);
- Choleski_decomposition chol (op->quad_);
-
- /*
- ugh.
- */
- H=chol.inverse();
- OK();
-
- degenerate_count_i_ = 0;
-}
-
-/** Find the optimum which is in the planes generated by the active
- constraints.
- */
-Vector
-Active_constraints::find_active_optimum (Vector g)
-{
- return H*g;
-}
+++ /dev/null
-/*
- spring-spacer.cc -- implement Spring_spacer
-
- source file of the GNU LilyPond music typesetter
-
- (c) 1996--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-*/
-
-
-#include <math.h>
-#include <limits.h>
-#include "killing-cons.tcc"
-#include "spring-spacer.hh"
-#include "paper-column.hh"
-#include "debug.hh"
-#include "dimensions.hh"
-#include "qlp.hh"
-#include "unionfind.hh"
-#include "idealspacing.hh"
-#include "pointer.tcc"
-#include "score-column.hh"
-#include "paper-def.hh"
-#include "column-x-positions.hh"
-#include "main.hh"
-
-Vector
-Spring_spacer::default_solution() const
-{
- return try_initial_solution();
-}
-
-Score_column*
-Spring_spacer::scol_l (int i)
-{
- return dynamic_cast<Score_column*>(cols_[i].pcol_l_);
-}
-
-const Real COLFUDGE=1e-3;
-template class P<Real>; // ugh.
-
-bool
-Spring_spacer::contains_b (Paper_column const *w)
-{
- for (int i=0; i< cols_.size(); i++)
- if (cols_[i].pcol_l_ == w)
- return true;
- return false;
-}
-
-
-void
-Spring_spacer::OK() const
-{
-#ifndef NDEBUG
- for (int i = 1; i < cols_.size(); i++)
- assert (cols_[i].rank_i_ > cols_[i-1].rank_i_);
-#endif
-}
-
-/**
- Make sure no unconnected columns happen.
- */
-void
-Spring_spacer::handle_loose_cols()
-{
- Union_find connected (cols_.size());
- Array<int> fixed;
-
- for (Cons<Idealspacing> *i = ideal_p_list_; i; i = i->next_)
- {
- connected.connect (i->car_->cols_drul_[LEFT],i->car_->cols_drul_[RIGHT]);
- }
- for (int i = 0; i < cols_.size(); i++)
- if (cols_[i].fixed_b())
- fixed.push (i);
- for (int i=1; i < fixed.size(); i++)
- connected.connect (fixed[i-1], fixed[i]);
-
- /*
- If columns do not have spacing information set, we need to supply our own.
- */
- Real d = default_space_f_;
- for (int i = cols_.size(); i--;)
- {
- if (! connected.equiv (fixed[0], i))
- {
- connected.connect (i-1, i);
- connect (i-1, i, d, 1.0);
- }
- }
-}
-
-bool
-Spring_spacer::check_constraints (Vector v) const
-{
- int dim=v.dim();
- assert (dim == cols_.size());
- DOUT << "checking " << v;
- for (int i=0; i < dim; i++)
- {
- if (cols_[i].fixed_b() &&
- abs (cols_[i].fixed_position() - v (i)) > COLFUDGE)
- {
- DOUT << "Fixpos broken\n";
- return false;
- }
- Array<Spacer_rod> const &rods (cols_[i].rods_[RIGHT]);
- for (int j =0; j < rods.size (); j++)
- {
- int other =rods[j].other_idx_;
- Real diff =v (other) - v (i) ;
- if (COLFUDGE +diff < rods[j].distance_f_)
- {
- DOUT << "i, other_i: " << i << " " << other << '\n';
- DOUT << "dist, minimal = " << diff << " "
- << rods[j].distance_f_ << '\n';
- return false;
- }
- }
-
- }
- return true;
-}
-
-/** try to generate a solution which obeys the min
- distances and fixed positions
- */
-Vector
-Spring_spacer::try_initial_solution() const
-{
- Vector v;
- if (!try_initial_solution_and_tell (v))
- {
- warning (_ ("I'm too fat; call Oprah"));
- }
- return v;
-
-}
-
-bool
-Spring_spacer::try_initial_solution_and_tell (Vector &v) const
-{
- int dim=cols_.size();
- bool succeeded = true;
- Vector initsol (dim);
-
- assert (cols_[0].fixed_b ());
- DOUT << "fixpos 0 " << cols_[0].fixed_position ();
- for (int i=0; i < dim; i++)
- {
- Real min_x = i ? initsol (i-1) : cols_[0].fixed_position ();
- Array<Spacer_rod> const &sr_arr(cols_[i].rods_[LEFT]);
- for (int j=0; j < sr_arr.size (); j++)
- {
- min_x = min_x >? (initsol (sr_arr[j].other_idx_) + sr_arr[j].distance_f_);
- }
- initsol (i) = min_x;
-
- if (cols_[i].fixed_b())
- {
- initsol (i)=cols_[i].fixed_position();
- if (initsol (i) < min_x )
- {
- DOUT << "failing: init, min : " << initsol (i) << " " << min_x << '\n';
- initsol (i) = min_x;
- succeeded = false;
- }
- }
- }
- v = initsol;
-
- DOUT << "tried and told solution: " << v;
- if (!succeeded)
- DOUT << "(failed)\n";
- return succeeded;
-}
-
-
-
-// generate the matrices
-void
-Spring_spacer::make_matrices (Matrix &quad, Vector &lin, Real &c) const
-{
- quad.fill (0);
- lin.fill (0);
- c = 0;
-
- for (Cons<Idealspacing> *p =ideal_p_list_; p; p = p->next_)
- {
- Idealspacing *i = p->car_;
- int l = i->cols_drul_[LEFT];
- int r = i->cols_drul_[RIGHT];
-
- quad (r,r) += i->hooke_f_;
- quad (r,l) -= i->hooke_f_;
- quad (l,r) -= i->hooke_f_;
- quad (l,l) += i->hooke_f_;
-
- lin (r) -= i->space_f_*i->hooke_f_;
- lin (l) += i->space_f_*i->hooke_f_;
-
- c += sqr (i->space_f_);
- }
-
- if (quad.dim() > 10)
- quad.set_band();
-
-
-}
-
-void
-Spring_spacer::set_fixed_cols (Mixed_qp &qp) const
-{
- for (int j=0; j < cols_.size(); j++)
- if (cols_[j].fixed_b())
- qp.add_fixed_var (j,cols_[j].fixed_position());
-}
-
-// put the constraints into the LP problem
-void
-Spring_spacer::make_constraints (Mixed_qp& lp) const
-{
- int dim=cols_.size();
-
- for (int j=0; j < dim -1; j++)
- {
- Array<Spacer_rod> const&rod_arr (cols_[j].rods_[RIGHT]);
- for (int i = 0; i < rod_arr.size (); i++)
- {
- Vector c1(dim);
- c1(rod_arr[i].other_idx_)=1.0 ;
- c1(j)=-1.0 ;
-
- lp.add_inequality_cons (c1, rod_arr[i].distance_f_);
- }
- }
-}
-
-
-Real
-Spring_spacer::calculate_energy_f (Vector solution) const
-{
- Real e = 0.0;
- for (Cons<Idealspacing>*p =ideal_p_list_; p; p = p->next_)
- {
- Idealspacing * i = p->car_;
- e += i->energy_f(solution(i->cols_drul_[RIGHT]) - solution(i->cols_drul_[LEFT]));
- }
-
- return e;
-}
-void
-Spring_spacer::lower_bound_solution (Column_x_positions*positions) const
-{
- Mixed_qp lp (cols_.size());
- make_matrices (lp.quad_,lp.lin_, lp.const_term_);
- set_fixed_cols (lp);
-
- Vector start (cols_.size());
- start.fill (0.0);
- Vector solution_vec (lp.solve (start));
- for (int i=0; i < solution_vec.dim (); i++)
- solution_vec(i) += indent_f_;
-
- DOUT << "Lower bound sol: " << solution_vec;
- positions->energy_f_ = calculate_energy_f (solution_vec);
- positions->config_ = solution_vec;
- positions->satisfies_constraints_b_ = check_constraints (solution_vec);
-}
-
-Spring_spacer::Spring_spacer ()
-{
- ideal_p_list_ =0;
- energy_normalisation_f_ = 1.0;
-}
-
-void
-Spring_spacer::solve (Column_x_positions*positions) const
-{
- Vector solution_try;
-
- bool constraint_satisfaction = try_initial_solution_and_tell (solution_try);
- if (constraint_satisfaction)
- {
- Mixed_qp lp (cols_.size());
- make_matrices (lp.quad_,lp.lin_, lp.const_term_);
- make_constraints (lp);
- set_fixed_cols (lp);
-
-
- Vector solution_vec (lp.solve (solution_try));
- for (int i=0; i < solution_vec.dim (); i++)
- solution_vec(i) += indent_f_;
-
-
- positions->satisfies_constraints_b_ = check_constraints (solution_vec);
- if (!positions->satisfies_constraints_b_)
- {
- warning (_("Solution doesn't satisfy constraints"));
- }
-
- positions->energy_f_ = calculate_energy_f (solution_vec);
- positions->config_ = solution_vec;
- }
- else
- {
- positions->set_stupid_solution (solution_try);
- }
-
-}
-
-/**
- add one column to the problem.
-
- TODO: ugh merge with add_columns.
-*/
-void
-Spring_spacer::add_column (Paper_column *col, bool fixed, Real fixpos)
-{
- Column_info c (col,(fixed)? &fixpos : 0);
- int this_rank = cols_.size();
- c.rank_i_ = this_rank;
-
- for (int i=0; i < col->minimal_dists_arr_drul_[LEFT].size (); i++)
- {
- Column_rod &cr = col->minimal_dists_arr_drul_[LEFT][i];
- int left_idx = cr.other_l_->rank_i () - cols_[0].pcol_l_->rank_i ();
- if (left_idx < 0)
- continue;
-
- if (cols_[left_idx].pcol_l_ != cr.other_l_)
- continue;
-
- Spacer_rod l_rod;
- l_rod.distance_f_ = cr.distance_f_;
- l_rod.other_idx_ = left_idx;
- c.rods_[LEFT].push (l_rod);
-
- Spacer_rod r_rod;
- r_rod.distance_f_ = cr.distance_f_;
- r_rod.other_idx_ = this_rank;
- cols_[left_idx].rods_[RIGHT].push (r_rod);
- }
-
- for (int i=0; i < col->spring_arr_drul_[LEFT].size (); i++)
- {
- Column_spring &cr = col->spring_arr_drul_[LEFT][i];
- int idx = cr.other_l_->rank_i () - cols_[0].pcol_l_->rank_i ();
- if (idx < 0)
- continue;
-
- if (cols_[idx].pcol_l_ != cr.other_l_)
- continue;
- connect (idx, this_rank, cr.distance_f_, cr.strength_f_);
- }
-
- cols_.push (c);
-}
-
-
-void
-Spring_spacer::add_columns (Link_array<Paper_column> cols)
-{
- energy_normalisation_f_ = sqr (line_len_f_);
- add_column (cols[0], true, 0.0);
- for (int i=1; i< cols.size ()-1; i++)
- add_column (cols[i],false,0.0);
-
- if (line_len_f_ > 0)
- add_column (cols.top (), true, line_len_f_);
- else
- add_column (cols.top (), false, 0);
-}
-
-
-
-void
-Spring_spacer::print() const
-{
-#ifndef NPRINT
- for (int i=0; i < cols_.size(); i++)
- {
- DOUT << "col " << i << " ";
- cols_[i].print();
- }
-
- for (Cons<Idealspacing> *p =ideal_p_list_; p; p = p->next_)
- {
- p->car_->print();
- }
-#endif
-}
-
-
-void
-Spring_spacer::connect (int i, int j, Real d, Real h)
-{
- if (d > 100 CM)
- {
- programming_error( _f ("Improbable distance: %f point, setting to 10 mm", d));
- d = 1 CM;
- }
- if(d < 0)
- {
- programming_error (_ ("Negative distance, setting to 10 mm"));
- d = 10 MM;
- }
-
- assert(h >=0);
-
- Idealspacing * s = new Idealspacing;
-
- s->cols_drul_[LEFT] = i ;
- s->cols_drul_[RIGHT] = j;
- s->space_f_ = d;
- s->hooke_f_ = h;
-
- ideal_p_list_ = new Killing_cons<Idealspacing> (s, ideal_p_list_);
-}
-
-void
-Spring_spacer::prepare()
-{
- handle_loose_cols();
- print();
-}
-
-
-Spring_spacer::~Spring_spacer()
-{
- delete ideal_p_list_;
-}
projects / lilypond.git / commitdiff