flower-1.0.2

diff --git a/flower/Sources.make b/flower/Sources.make
new file mode 100644 (file)
index 0000000..18fb18c
--- /dev/null
+++ b/flower/Sources.make
@@ -0,0 +1,12 @@
+
+cc=lgetopt.cc   string.cc dataf.cc textdb.cc unionfind.cc  \
+       smat.cc matrix.cc choleski.cc vector.cc dstream.cc\
+       matdebug.cc
+
+templatecc=cursor.cc list.cc tsmat.cc 
+inl=findcurs.inl link.inl list.inl
+hh=cursor.hh cursor.inl lgetopt.hh link.hh list.hh dstream.hh \
+       string.hh stringutil.hh vray.hh textdb.hh textstr.hh  assoc.hh\
+       findcurs.hh unionfind.hh compare.hh handle.hh matrix.hh\
+       smat.hh vsmat.hh  vector.hh  real.hh choleski.hh\
+       tsmat.hh tvsmat.hh

diff --git a/flower/choleski.cc b/flower/choleski.cc
new file mode 100644 (file)
index 0000000..82740a2
--- /dev/null
+++ b/flower/choleski.cc
@@ -0,0 +1,97 @@
+#include "choleski.hh"
+const Real EPS = 1e-7;         // so sue me. Hard coded
+
+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);
+    }
+
+    // backward subst
+    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;
+}
+
+/*
+  Standard matrix algorithm.
+  */
+
+Choleski_decomposition::Choleski_decomposition(Matrix P)
+    : L(P.dim()), D(P.dim())
+{
+    int n = P.dim();
+    assert((P-P.transposed()).norm()/P.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() / 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()/original.norm() < EPS);
+#endif
+    
+    return invm;
+}
+
+
+
+

diff --git a/flower/matrix.cc b/flower/matrix.cc
new file mode 100644 (file)
index 0000000..1afe595
--- /dev/null
+++ b/flower/matrix.cc
@@ -0,0 +1,260 @@
+#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();
+}
+

diff --git a/flower/vector.cc b/flower/vector.cc
new file mode 100644 (file)
index 0000000..2cbfdfc
--- /dev/null
+++ b/flower/vector.cc
@@ -0,0 +1,20 @@
+#include "vector.hh"
+#include "string.hh"
+
+Vector::Vector(const Vector&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;
+}

diff --git a/flower/vector.hh b/flower/vector.hh
new file mode 100644 (file)
index 0000000..1929c67
--- /dev/null
+++ b/flower/vector.hh
@@ -0,0 +1,111 @@
+#ifndef VECTOR_HH
+#define VECTOR_HH
+
+#include <math.h>
+#include "real.hh"
+#include "vray.hh"
+
+class Dstream;
+class String;
+void set_matrix_debug(Dstream&ds);
+
+/// a row of numbers
+class Vector  {
+    svec<Real> dat;
+public:
+    void OK() const { dat.OK();}
+    int dim() const { return dat.sz(); }
+    Vector() { }
+    Vector(const Vector&n);
+    Vector(int n) {
+       dat.set_size(n);
+       fill(0);
+    }
+    void insert(Real v, int i) {
+       dat.insert(v,i);
+    }
+    void del(int i) { dat.del(i); }
+    operator String() const;
+    void fill(Real r) {
+       for (int i=0; i < dim(); i++)
+           dat[i] =r;
+    }
+
+    void operator +=(Vector v) {
+       assert(v.dim() == dim());
+       for (int i=0; i < dim(); i++)
+           dat[i] += v.dat[i];
+    }
+    
+    void operator /=(Real a) {
+       (*this) *= 1/a;
+    }
+
+    void operator *=(Real a) {
+       for (int i=0; i < dim(); i++)
+           dat[i] *= a;
+    }
+
+    void operator -=(Vector v) {
+       assert(v.dim() == dim());
+       for (int i=0; i < dim(); i++)
+           dat[i] -= v(i);     
+    }
+
+    Real &operator()(int i) { return dat[i]; }
+    Real operator()(int i) const { return dat[i]; }
+    Real elem(int i) { return dat[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;
+    }
+    Vector operator-() const;
+    Real norm() {
+       return sqrt(norm_sq() );
+    }
+    Real norm_sq() {
+       return ((*this) * (*this));
+    }
+    operator svec<Real> () { return dat; }
+    void print() const;
+    /// set to j-th element of unit-base
+    void set_unit(int j) ;
+};
+/**
+    a vector. Storage is handled in svec, Vector only does the mathematics.
+ */
+
+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