[Toon-members] TooN Cholesky.h helpers.h

[Gerhard Reitmayr]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Gerhard Reitmayr <gerhard>      07/05/14 21:04:41

Modified files:
        .              : Cholesky.h helpers.h 

Log message:
        additional versions of various functions with support for DynamicMatrix

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.15&r2=1.16
http://cvs.savannah.gnu.org/viewcvs/TooN/helpers.h?cvsroot=toon&r1=1.18&r2=1.19

Patches:
Index: Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Cholesky.h,v
retrieving revision 1.15
retrieving revision 1.16
diff -u -b -r1.15 -r1.16
--- Cholesky.h  11 Feb 2007 13:20:11 -0000      1.15
+++ Cholesky.h  14 May 2007 21:04:40 -0000      1.16
@@ -48,15 +48,25 @@
                x[I] = v[I] - Dot<0,I-1>::eval(L[I], x);
                Forwardsub_L<N,I+1>::eval(L, v, x);
            }
+           template <class A1, class A2, class Vec> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const DynamicVector<A2>& v, Vec & x) {
+               x[I] = v[I] - Dot<0,I-1>::eval(L[I], x);
+               Forwardsub_L<N,I+1>::eval(L, v, x);
+           }
        };
        template <int N> struct Forwardsub_L<N,0> {
            template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v, 
FixedVector<N,A3>& x) {
                x[0] = v[0];
                Forwardsub_L<N,1>::eval(L, v, x);
            }
+           template <class A1, class A2, class Vec> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const DynamicVector<A2>& v, Vec & x) {
+               assert(v.size() == N && x.size() == N);
+               x[0] = v[0];
+               Forwardsub_L<N,1>::eval(L, v, x);
+           }
        };
        template <int N> struct Forwardsub_L<N,N> {
            template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v, 
FixedVector<N,A3>& x) {}
+           template <class A1, class A2, class Vec> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const DynamicVector<A2>& v, Vec & x) {}
        };
        
        //
@@ -225,6 +235,8 @@
     class Cholesky {
     public:
 
+       Cholesky() {}
+
        template<class Accessor>
        Cholesky(const FixedMatrix<N,N,Accessor>& m){
            compute(m);
@@ -329,15 +341,43 @@
                }
            }
        }
+
+       template <class F, class A1, class M2> void transform_inverse(const 
DynamicMatrix<A1>& J, M2 & T) const {
+           const int M = J.num_rows();
+           assert( T.num_rows() == M && T.num_cols() == M);
+           Matrix<> L_inv_JT(M,N);
+           for (int i=0; i<M; ++i)
+               util::Forwardsub_L<N>::eval(L, J[i], L_inv_JT[i]);
+           for (int i=0; i<M; ++i) {
+               F::eval(T[i][i],util::Dot3<0,N-1>::eval(L_inv_JT[i], 
L_inv_JT[i], invdiag));
+               for (int j=i+1; j<M; ++j) {
+                   const double x = util::Dot3<0,N-1>::eval(L_inv_JT[i], 
L_inv_JT[j], invdiag);
+                   F::eval(T[i][j],x);
+                   F::eval(T[j][i],x);
+               }
+           }
+       }
+
        template <int M, class A1, class A2> void transform_inverse(const 
FixedMatrix<M,N,A1>& J, FixedMatrix<M,M,A2>& T) const {
            transform_inverse<util::Assign>(J,T);
        }
+
+       template <class A1, class M2> void transform_inverse(const 
DynamicMatrix<A1>& J, M2 & T) const {
+           transform_inverse<util::Assign>(J,T);
+       }
+
        template <int M, class A> Matrix<M> transform_inverse(const 
FixedMatrix<M,N,A>& J) const {
            Matrix<M> T;
            transform_inverse(J,T);
            return T;
        }
 
+       template <class A> Matrix<> transform_inverse(const DynamicMatrix<A>& 
J) const {
+           Matrix<> T(J.num_rows(), J.num_rows());
+           transform_inverse(J,T);
+           return T;
+       }
+
        template <class A1, class A2> inline 
        void inverse_times(const FixedVector<N,A1>& v, FixedVector<N,A2>& x) 
const
        {
@@ -397,6 +437,8 @@
     class Cholesky<-1> {
     public:
 
+    Cholesky(){}
+
        template<class Accessor>
        Cholesky(const DynamicMatrix<Accessor>& m) {
            compute(m);

Index: helpers.h
===================================================================
RCS file: /cvsroot/toon/TooN/helpers.h,v
retrieving revision 1.18
retrieving revision 1.19
diff -u -b -r1.18 -r1.19
--- helpers.h   30 Jan 2007 15:57:11 -0000      1.18
+++ helpers.h   14 May 2007 21:04:40 -0000      1.19
@@ -284,7 +284,7 @@
         }
      }
 
-     template <class A1, class A2, class MatM> inline void 
transformCovariance(const DynamicMatrix<A1>& A, const DynamicMatrix<A2>& B, 
MatM& M)
+     template <class F, class A1, class M2, class M3> inline void 
transformCovarianceUpper(const DynamicMatrix<A1>& A, const M2 & B, M3 & M)
      {
         const int R = A.num_rows();
         const int N = A.num_cols();     
@@ -294,17 +294,35 @@
                B.num_cols() == N);
         for (int i=0; i<R; ++i) {
             const Vector<> ABi = B * A[i];
-            M[i][i] = ABi * A[i];
-            for (int j=i+1; j<R; ++j)
-                M[j][i] = M[i][j] = ABi * A[j];
+            for (int j=i; j<R; ++j)
+                F::eval(M[i][j], ABi * A[j]);
+        }
+     }
+
+     template <class F, class A1, class M2, class MatM> inline void 
transformCovariance(const DynamicMatrix<A1> & A, const M2 & B, MatM& M)
+     {
+       const int R = A.num_rows();
+       const int N = A.num_cols();
+       assert(M.num_rows() == R &&
+              M.num_cols() == R &&
+              B.num_rows() == N &&
+              B.num_cols() == N);
+       for (int i=0; i<R; ++i) {
+            const Vector<> ABi = B * A[i];
+            F::eval(M[i][i], ABi * A[i]);
+            for (int j=i+1; j<R; ++j){
+               const double v = ABi * A[j];
+               F::eval(M[j][i], v);
+               F::eval(M[i][j], v);
+            }
         }
      }
  }
  
- template <class A1, class A2> Matrix<> inline transformCovariance(const 
DynamicMatrix<A1>& A, const DynamicMatrix<A2>& B)
+ template <class A1, class A2> Matrix<> inline transformCovariance(const 
DynamicMatrix<A1> & A, const A2 & B)
  {
      Matrix<> M(A.num_rows(), A.num_rows());
-     util::transformCovariance(A,B,M);
+     util::transformCovariance<util::Assign>(A,B,M);
      return M;
  }
 
@@ -349,6 +367,19 @@
         }
  }
 
+ template <class A, class B, class Vec> inline void add_product(const 
DynamicMatrix<A>& Ma, const DynamicVector<B>& Mb, Vec & r)
+ {
+     const int M=Ma.num_rows();
+     const int N=Ma.num_cols();
+     assert(N == Mb.size());
+     for (int i=0; i<M; ++i) {
+       double sum = 0;
+       for (int k=0; k<N; ++k)
+           sum += Ma(i,k)*Mb[k];
+       r[i] += sum;
+     }
+ }
+
  template <int M, int N, int C, class A1, class A2, class Mat> void 
matrix_multiply(const FixedMatrix<M,N,A1>& A, const FixedMatrix<N,C,A2>& B, 
Mat& R)
  {
      util::matrix_multiply<M,N,C>(A,B,R);