lilypond-0.0.1

diff --git a/choleski.cc b/choleski.cc
new file mode 100644 (file)
index 0000000..db3fe82
--- /dev/null
+++ b/choleski.cc
@@ -0,0 +1,85 @@
+#include "choleski.hh"
+
+Vector
+Choleski_decomposition::solve(Vector rhs)const
+{
+    int n= rhs.dim();
+    assert(n == L.dim());
+    Vector y(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++) {
+       assert(D(i));
+       y(i) /= D(i);
+    }
+    Vector x(n);
+    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);
+    }
+    return x;
+}
+
+Choleski_decomposition::Choleski_decomposition(Matrix P)
+    : L(P.dim()), D(P.dim())
+{
+    int n = P.dim();
+    assert((P-P.transposed()).norm() < EPS);
+    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;
+    }
+
+    #ifdef NDEBUG
+    assert((original()-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);
+    for (int i = 0; i < n; i++) {
+       e_i.set_unit(i);
+       Vector inv(solve(e_i));
+       for (int j = 0 ; j<n; j++)
+           invm(i,j) = inv(j);
+    }
+    
+#ifdef NDEBUG
+    Matrix I1(n), I2(original());
+    I1.unit();
+    assert((I1-original()*invm).norm() < EPS);
+    #endif
+    
+    return invm;
+}

diff --git a/choleski.hh b/choleski.hh
new file mode 100644 (file)
index 0000000..c6cb917
--- /dev/null
+++ b/choleski.hh
@@ -0,0 +1,46 @@
+#ifndef CHOLESKI_HH
+#define CHOLESKI_HH
+
+#include "matrix.hh"
+
+struct Choleski_decomposition {
+
+    /// lower triangle of Choleski decomposition
+    Matrix L;
+
+    /// diagonal 
+    Vector D;
+    ///Create decomposition of P
+    Choleski_decomposition(Matrix P);
+    /**
+    PRE
+    P needs to be symmetric positive definite
+    */
+    
+    Vector solve(Vector rhs) const;
+    Vector operator * (Vector rhs) const { return solve (rhs); }
+    /**
+    solve Px = rhs
+    */
+
+    Matrix inverse() const;
+    /**
+    return the inverse of the matrix P.
+    */
+
+    Matrix original() const;
+    /**
+    return P,  calc'ed from L and D
+    */
+        
+};
+/**
+    structure for using the LU decomposition of a positive definite .
+
+    #P# is split  into
+
+    LD transpose(L)
+    */
+    
+
+#endif

diff --git a/matrix.cc b/matrix.cc
new file mode 100644 (file)
index 0000000..ab59a87
--- /dev/null
+++ b/matrix.cc
@@ -0,0 +1,281 @@
+#include "matrix.hh"
+#include "string.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);
+}
+
+//inline
+Real
+Matrix::operator()(int i,int j) const
+{
+    assert(i >= 0 && j >= 0);
+    assert(i < rows() && j < cols());
+    return dat->elem(i,j);
+}
+
+//inline
+Real &
+Matrix::operator()(int i, int j)
+{
+    assert(i >= 0 && j >= 0);
+    assert(i < rows() && j < cols());
+    return dat->elem(i,j);
+}
+
+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+=(const Matrix&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-=(const Matrix&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=(const Matrix&m)
+{
+    if (&m == this)
+       return ;
+    delete dat;
+    dat = m.dat->clone();
+}
+    
+Matrix::Matrix(const Matrix &m)
+{
+    m.OK();
+    
+    dat = m.dat->clone();
+}
+
+
+Matrix::Matrix(int n, int m)
+{
+    dat = virtual_smat::get_full(n,m);
+    fill(0);
+}
+
+Matrix::Matrix(int n)
+{
+    dat = virtual_smat::get_full(n,n);
+    fill(0);
+}
+
+Matrix::Matrix(Vector v, Vector w)
+{   
+    dat = virtual_smat::get_full(v.dim(), w.dim());
+    for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,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(const Vector& 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 *(const Vector& 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;
+}
+
+void
+Matrix::transpose()            // delegate to storage?
+{
+    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;
+}
+
+
+/* should do something smarter: bandmatrix * bandmatrix is also banded matrix.  */
+Matrix
+operator *(const Matrix &m1, const Matrix &m2)
+{
+    Matrix result(m1.rows(), m2.cols());
+    result.set_product(m1,m2);
+    return result;
+}
+
+void
+Matrix::set_product(const Matrix &m1, const Matrix &m2)
+{
+    assert(m1.cols()==m2.rows());
+    assert(cols()==m2.cols() && rows()==m1.rows());
+    
+    for (int i=0, j=0; dat->mult_ok(i,j);
+        dat->mult_next(i,j)) {
+       Real r=0.0;
+       for (int k = 0; k < m1.cols(); k++)
+           r += m1(i,k)*m2(k,j);
+       dat->elem(i,j)=r;
+    }
+}
+
+void
+Matrix::insert_row(Vector v, int k)
+{
+    assert(v.dim()==cols());
+    dat->insert_row(k);
+    for (int j=0; j < cols(); j++)
+       dat->elem(k,j)=v(j);
+}
+
+
+void
+Matrix::swap_columns(int c1, int c2)
+{
+    assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0);
+    for (int i=0; i< rows(); 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);
+    for (int i=0; i< cols(); 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::operator String() const
+{
+    String s("matrix {\n");
+    for (int i=0; i< rows(); i++){
+       for (int j = 0; j < cols(); j++) {
+           s+= String(dat->elem(i,j), "%6f ");
+       }
+       s+="\n";
+    }
+    s+="}\n";
+    return s;
+}
+
+
+
+#include "debug.hh"
+void
+Matrix::print() const
+{
+    mtor << *this;
+}

diff --git a/vector.cc b/vector.cc
new file mode 100644 (file)
index 0000000..be24232
--- /dev/null
+++ b/vector.cc
@@ -0,0 +1,41 @@
+#include "debug.hh"
+#include "vector.hh"
+#include "string.hh"
+
+Vector::Vector(const Vector&n)
+          :dat(n.dat)
+     // this makes GCC 272 barf
+{
+    //dat = n.dat;
+}    
+
+Vector::operator String() const
+{
+    int i=0;
+    String s("vector [");
+    for (; i < dim(); i++) {
+       s += String(dat[i], "%6f") + ' ';
+    }
+    s+="]";
+    return s;
+}
+
+
+void
+Vector::print() const
+{
+    mtor << *this<<'\n';
+}
+
+Vector
+Vector::operator-() const
+{
+    Vector v(*this); v*=-1; return v;
+}
+
+void
+Vector::set_unit(int j)
+{
+    fill(0.0);
+    dat[j] = 1.0;
+}