[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN Cholesky.h helpers.h
From: |
Ethan Eade |
Subject: |
[Toon-members] TooN Cholesky.h helpers.h |
Date: |
Sat, 22 Jul 2006 12:29:28 +0000 |
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN Cholesky.h helpers.h,
Ethan Eade <=