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

[Ethan Eade]

CVSROOT:        /sources/toon
Module name:    TooN
Changes by:     Ethan Eade <ethaneade>  06/07/22 12:29:28

Modified files:
        .              : Cholesky.h helpers.h 

Log message:
        Fixed some methods for dynamic Cholesky; added dynamic 
transformCovariance.   

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.10&r2=1.11
http://cvs.savannah.gnu.org/viewcvs/TooN/helpers.h?cvsroot=toon&r1=1.14&r2=1.15

Patches:
Index: Cholesky.h
===================================================================
RCS file: /sources/toon/TooN/Cholesky.h,v
retrieving revision 1.10
retrieving revision 1.11
diff -u -b -r1.10 -r1.11
--- Cholesky.h  16 Jul 2006 19:45:01 -0000      1.10
+++ Cholesky.h  22 Jul 2006 12:29:28 -0000      1.11
@@ -102,6 +102,36 @@
            }
        }
        
+       template <class A1, class A2, class A3> inline void 
cholesky_inverse(const DynamicMatrix<A1>& L, const DynamicVector<A2>& invdiag, 
DynamicMatrix<A3>& I)
+       {
+           assert(L.num_rows() == L.num_cols() && 
+                  L.num_rows() == invdiag.size() && 
+                  L.num_rows() == I.num_rows() && 
+                  I.num_rows() == I.num_cols());
+           int S = L.num_rows();
+           for (int col = 0; col<S; col++) {
+               Vector<> t(S),x(S);
+               t[col] = invdiag[col];
+               for (int i=col+1; i<S; i++) {
+                   double psum = 0;
+                   for (int j=col; j<i; j++)
+                       psum -= L[i][j]*t[j];               
+                   t[i] = invdiag[i] * psum;
+               }
+               for (int i=S-1; i>col; i--) {
+                   double psum = t[i];
+                   for (int j=i+1; j<S; j++)
+                       psum -= L[j][i]*x[j];
+                   I[i][col] = I[col][i] = x[i] = invdiag[i] * psum;
+               }
+               double psum = t[col];
+               for (int j=col+1; j<S; j++)
+                   psum -= L[j][col]*x[j];             
+               I[col][col] = invdiag[col]*psum;
+           }
+       }
+
+       
     }
 
     template <int Size=-1>
@@ -220,32 +250,14 @@
            return det*det;
        }
 
-       Matrix<> get_inverse() const {
-           int Size = L.num_rows();
-           Matrix<> M(Size, Size);
-           // compute inverse(L)
-           for (int r=0; r<Size; r++) {
-               M[r][r] = invdiag[r];
-               for (int i=r+1; i<Size; i++) {
-                   double b = -L[i][r]*M[r][r];
-                   for (int j=r+1; j<i;j++)
-                       b -= L[i][j]*M[r][j];
-                   M[r][i] = b*invdiag[i];
-               }      
-           }
-           // multiply (inverse(L))^T * inverse(L)
-           for (int i=0;i<Size; i++) {
-               double s = 0;
-               for (int k=i; k<Size; k++)
-                   s += M[i][k]*M[i][k];
-               M[i][i] = s;
-               for (int j=i+1; j<Size; j++) {
-                   double s = 0;
-                   for (int k=j; k<Size; k++)
-                       s += M[j][k]*M[i][k];
-                   M[i][j] = M[j][i] = s;
-               }
+       template <class Mat> void get_inverse(Mat& M) const {
+           assert(M.num_rows() == M.num_cols() && M.num_rows() == 
L.num_rows());
+           util::cholesky_inverse(L, invdiag, M);
            }
+
+       Matrix<> get_inverse() const {
+           Matrix<> M(L.num_rows(), L.num_rows());
+           get_inverse(M);
            return M;
        }
 

Index: helpers.h
===================================================================
RCS file: /sources/toon/TooN/helpers.h,v
retrieving revision 1.14
retrieving revision 1.15
diff -u -b -r1.14 -r1.15
--- helpers.h   3 Jul 2006 10:03:31 -0000       1.14
+++ helpers.h   22 Jul 2006 12:29:28 -0000      1.15
@@ -240,9 +240,14 @@
 }
 
  namespace util {
-     template <int R, int N, class Accessor1, class Accessor2> Matrix<R,R> 
transformCovariance(const FixedMatrix<R,N,Accessor1>& A, const 
FixedMatrix<N,N,Accessor2>& B)
+     template <class MatA, class MatB, class MatM> void 
transformCovariance(const MatA& A, const MatB& B, MatM& M)
      {
-        Matrix<R> M;
+        assert(M.num_rows() == M.num_cols() && 
+               M.num_rows() == A.num_rows() && 
+               A.num_cols() == B.num_rows() && 
+               B.num_rows() == B.num_cols());
+        const int R = A.num_rows();
+        const int N = A.num_cols();
         for (int i=0; i<R; i++) {
             double sum = 0;
             for (int k=0; k<N; k++) {
@@ -259,13 +264,33 @@
                 M[i][j] = M[j][i] = sum;
             }
         }
+     }
+
+     template <int R, int N, class Accessor1, class Accessor2> Matrix<R,R> 
transformCovariance(const FixedMatrix<R,N,Accessor1>& A, const 
FixedMatrix<N,N,Accessor2>& B)
+     {
+        Matrix<R> M;
+        transformCovariance(A,B,M);
+        return M;
+     }
+ }
+ 
+ template <class A1, class A2> Matrix<> transformCovariance(const 
DynamicMatrix<A1>& A, const DynamicMatrix<A2>& B)
+ {
+     Matrix<> M(A.num_rows(), A.num_rows());
+     util::transformCovariance(A,B,M);
         return M;
      }
+
+ template <int R, int N, class Accessor1, class Accessor2> inline Matrix<R> 
transformCovariance(const FixedMatrix<R,N,Accessor1>& A, const 
FixedMatrix<N,N,Accessor2>& B)
+ {
+     Matrix<R> M;
+     util::transformCovariance(A,B,M);
+     return M;
  }
 
- template <int R, int N, class Accessor1, class Accessor2> inline Matrix<R,R> 
transformCovariance(const FixedMatrix<R,N,Accessor1>& A, const 
FixedMatrix<N,N,Accessor2>& B)
+ template <class MatA, class MatB, class MatM> void transformCovariance(const 
MatA& A, const MatB& B, MatM& M)
  {
-     return util::transformCovariance(A,B);
+     util::transformCovariance(A,B,M);
  }
  
 #ifndef TOON_NO_NAMESPACE