Changeset 32 for bdm/math

Timestamp:

03/03/08 13:00:32 (16 years ago)

Author:

smidl

Message:

test KF : estimation of R in KF is not possible! Likelihood of y_t is growing when R -> 0

Location:

bdm/math

Files:

• libDC.cpp (modified) (9 diffs)
• libDC.h (modified) (12 diffs)

bdm/math/libDC.cpp

@@ -17 +17 @@
 void fsqmat::mult_sym ( const mat &C) {M=C *M*C.T();};
 void fsqmat::mult_sym_t ( const mat &C) {M=C.T() *M*C;};
-void fsqmat::mult_sym ( const mat &C, fsqmat &U) { U.M = ( C *(M*C.T()) );};
-void fsqmat::mult_sym_t ( const mat &C, fsqmat &U) { U.M = ( C.T() *(M*C) );};
+void fsqmat::mult_sym ( const mat &C, fsqmat &U) const { U.M = ( C *(M*C.T()) );};
+void fsqmat::mult_sym_t ( const mat &C, fsqmat &U) const { U.M = ( C.T() *(M*C) );};
 void fsqmat::inv ( fsqmat &Inv ) {mat IM = itpp::inv ( M ); Inv=IM;};
 void fsqmat::clear() {M.clear();};
-fsqmat::fsqmat ( const mat &M0 )
+fsqmat::fsqmat ( const mat &M0 ) : sqmat(M0.cols())
 {
         it_assert_debug ( ( M0.cols() ==M0.rows() ),"M0 must be square" );
-        M=M0;dim=M0.cols();
+        M=M0;
 };
 
-fsqmat::fsqmat() {};
-
-fsqmat::fsqmat(const int dim0): M(dim0,dim0) {};
-
-std::ostream &operator<< ( std::ostream &os,  fsqmat &ld ) {
+//fsqmat::fsqmat() {};
+
+fsqmat::fsqmat(const int dim0): sqmat(dim0), M(dim0,dim0) {};
+
+std::ostream &operator<< ( std::ostream &os, const fsqmat &ld ) {
         os << ld.M << endl;
         return os;
@@ -37 +37 @@
 
 
-ldmat::ldmat( const mat &exL, const vec &exD ) {
+ldmat::ldmat( const mat &exL, const vec &exD ) : sqmat(exD.length()) {
         D = exD;
         L = exL;
-        dim = exD.length();
-}
-
-ldmat::ldmat() {
-        dim = 0;
-}
-
-ldmat::ldmat(const int dim0): D(dim0),L(dim0,dim0) {
-        dim = dim0;
-}
-
-ldmat::ldmat(const vec D0) {
+}
+
+ldmat::ldmat() :sqmat(0) {}
+
+ldmat::ldmat(const int dim0): sqmat(dim0), D(dim0),L(dim0,dim0) {}
+
+ldmat::ldmat(const vec D0):sqmat(D0.length()) {
         D = D0;
-        dim = D0.length();
         L = eye(dim);
 }
 
-ldmat::ldmat( const mat &V ) {
+ldmat::ldmat( const mat &V ):sqmat(V.cols()) {
 //TODO check if correct!! Based on heuristic observation of lu()
 
-        dim = V.cols();
         it_assert_debug( dim == V.rows(),"ldmat::ldmat matrix V is not square!" );
                 
@@ -86 +79 @@
 }
 
-std::ostream &operator<< ( std::ostream &os,  ldmat &ld ) {
+std::ostream &operator<< ( std::ostream &os, const ldmat &ld ) {
         os << "L:" << ld.L << endl;
         os << "D:" << ld.D << endl;
@@ -130 +123 @@
 void ldmat::clear(){L.clear(); for ( int i=0;i<L.cols();i++ ){L( i,i )=1;}; D.clear();}
 
-void ldmat::inv( ldmat &Inv ) {
+void ldmat::inv( ldmat &Inv ) const {
         int dim = D.length();
         Inv.clear();   //Inv = zero in LD
@@ -152 +145 @@
 }
 
-void ldmat::mult_sym( const mat &C, ldmat &U) {
+void ldmat::mult_sym( const mat &C, ldmat &U) const {
         mat A=L*C.T(); //could be done more efficiently using BLAS
         U.ldform(A,D);
 }
 
-void ldmat::mult_sym_t( const mat &C, ldmat &U) {
+void ldmat::mult_sym_t( const mat &C, ldmat &U) const {
         mat A=L*C;
 /*      vec nD=zeros(U.rows());
@@ -173 +166 @@
 }
 
-double ldmat::qform( const vec &v ) {
+double ldmat::qform( const vec &v ) const {
         double x = 0.0, sum;
         int i,j;
@@ -179 +172 @@
         for ( i=0; i<D.length(); i++ ) { //rows of L
                 sum = 0.0;
-                for ( j=0; j<=i; j++ ){sum+=L( i,j )*v( j );}
+                for ( j=i; j<=i; j++ ){sum+=L( i,j )*v( i );}
                 x +=D( i )*sum*sum;
         };
@@ -191 +184 @@
 }
 
-vec ldmat::sqrt_mult( const vec &x ) {
+vec ldmat::sqrt_mult( const vec &x ) const {
         int i,j;
         vec res( dim );
@@ -383 +376 @@
         }
 }
-
-

bdm/math/libDC.h

@@ -61 +61 @@
                 Used e.g. in generating normal samples.
                 */
-                virtual vec sqrt_mult (const vec &v ) =0;
+                virtual vec sqrt_mult (const vec &v ) const =0;
 
                 /*!
@@ -67 +67 @@
 
                 */
-                virtual double qform (const vec &v ) =0;
+                virtual double qform (const vec &v ) const =0;
 
 //      //! easy version of the
@@ -83 +83 @@
                 //! Destructor for future use;
                 virtual ~sqmat(){};
+                //! Default constructor
+                sqmat(const int dim0): dim(dim0){};
         protected:
                 int dim;
@@ -101 +103 @@
                 void mult_sym ( const mat &C);
                 void mult_sym_t ( const mat &C);
-                void mult_sym ( const mat &C, fsqmat &U);
-                void mult_sym_t ( const mat &C, fsqmat &U);
+                void mult_sym ( const mat &C, fsqmat &U) const;
+                void mult_sym_t ( const mat &C, fsqmat &U) const;
                 void clear();
 
@@ -124 +126 @@
 
                 double logdet() const {return log ( det ( M ) );};
-                double qform (const  vec &v ) {return ( v* ( M*v ) );};
-                vec sqrt_mult (const vec &v ) {it_error ( "not implemented" );return v;};
+                double qform (const  vec &v ) const {return ( v* ( M*v ) );};
+                vec sqrt_mult (const vec &v ) const {it_error ( "not implemented" );return v;};
 
                 fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;};
@@ -132 +134 @@
 //              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
 
-                friend std::ostream &operator<< ( std::ostream &os, fsqmat &sq );
+                friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );
 
 };
@@ -162 +164 @@
                 void add ( const ldmat &ld2, double w=1.0 );
                 double logdet() const;
-                double qform (const vec &v );
+                double qform (const vec &v ) const;
 //      sqmat& operator -= ( const sqmat & ld2 );
                 void clear();
                 int cols() const;
                 int rows() const;
-                vec sqrt_mult ( const vec &v );
+                vec sqrt_mult ( const vec &v ) const;
 
                 /*! \brief Matrix inversion preserving the chosen form.
@@ -174 +176 @@
 
                 */
-                virtual void inv ( ldmat &Inv );
+                virtual void inv ( ldmat &Inv ) const;
 
                 /*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class.
@@ -181 +183 @@
 
                 */
-                void mult_sym ( const mat &C, ldmat &U);
+                void mult_sym ( const mat &C, ldmat &U) const;
 
                 /*! \brief Symmetric multiplication of $U$ by a transpose of a general matrix $C$, result of which is stored in the current class.
@@ -188 +190 @@
 
                 */
-                void mult_sym_t ( const mat &C, ldmat &U);
+                void mult_sym_t ( const mat &C, ldmat &U) const;
 
 
@@ -200 +202 @@
                 */
                 void ldform (const mat &A,const vec &D0 );
+
+                //! Access functions
+                void setD (const vec &nD){D=nD;}
+                void setL (const vec &nL){L=nL;}
 
                 ldmat& operator += ( const ldmat &ldA );
@@ -205 +211 @@
                 ldmat& operator *= ( double x );
 
-                friend std::ostream &operator<< ( std::ostream &os, ldmat &sq );
+                friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );
+
 
         protected: