Index: bdm/stat/libEF.h
===================================================================
--- bdm/stat/libEF.h (revision 27)
+++ bdm/stat/libEF.h (revision 28)
@@ -26,9 +26,14 @@
 * More?...
 */
+
 class eEF : public epdf {
 
 public:
-	virtual void tupdate( double phi, mat &vbar, double nubar ) {};
-	virtual void dupdate( mat &v,double nu=1.0 ) {};
+	//! default constructor
+	eEF() :epdf() {};
+
+	virtual void tupdate ( double phi, mat &vbar, double nubar ) {};
+
+	virtual void dupdate ( mat &v,double nu=1.0 ) {};
 };
 
@@ -45,23 +50,45 @@
 */
 template<class sq_T>
+
 class enorm : public eEF {
+protected:
+	//! mean value
+	vec mu;
+	//! Covariance matrix in decomposed form
+	sq_T R;
+	//! Cache:  _iR = inv(R);
+	sq_T _iR;
+	//! indicator if \c _iR is chached
+	bool cached;
+	//! Random number generator for sampling; Fixme: what about *static* correct ?
+	Normal_RNG RNG;
+	//! dimension (redundant from rv.count() for easier coding )
 	int dim;
-	vec mu;
-	sq_T R;
-public:
-	enorm( RV &rv, vec &mu, sq_T &R );
-	enorm();
+public:
+	enorm() :eEF() {};
+
+	enorm ( RV &rv );
+	void set_parameters ( const vec &mu,const sq_T &R );
 	//! tupdate in exponential form (not really handy)
-	void tupdate( double phi, mat &vbar, double nubar );
-	void dupdate( mat &v,double nu=1.0 );
-	//!tupdate used in KF
-	void tupdate();
-	//!dupdate used in KF
-	double dupdate();
-	 
+	void tupdate ( double phi, mat &vbar, double nubar );
+	void dupdate ( mat &v,double nu=1.0 );
+
 	vec sample();
-	mat sample(int N);
-	double eval( const vec &val );
-	Normal_RNG RNG;
+	mat sample ( int N );
+	double eval ( const vec &val );
+	double evalpdflog ( const vec &val );
+
+//Access methods
+	//! returns a pointer to the internal mean value. Use with Care!
+	vec* _mu() {return &mu;}
+
+	//! returns pointers to the internal variance and its inverse. Use with Care!
+	void _R ( sq_T* &pR, sq_T* &piR ) {
+		pR=&R;
+		piR=&_iR;
+	}
+
+	//! set cache as inconsistent
+	void _cached ( bool what ) {cached=what;}
 };
 
@@ -70,4 +97,5 @@
 */
 template<class sq_T>
+
 class mlnorm : public mEF {
 	enorm<sq_T> epdf;
@@ -75,9 +103,9 @@
 public:
 	//! Constructor
-	mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R );
+	mlnorm ( RV &rv,RV &rvc, mat &A, sq_T &R );
 	//!Generate one sample of the posterior
-	vec samplecond( vec &cond, double &lik );
-	mat samplecond( vec &cond, vec &lik, int n );
-	void condition( vec &cond );
+	vec samplecond ( vec &cond, double &lik );
+	mat samplecond ( vec &cond, vec &lik, int n );
+	void condition ( vec &cond );
 };
 
@@ -85,17 +113,23 @@
 
 template<class sq_T>
-enorm<sq_T>::enorm( RV &rv, vec &mu0, sq_T &R0 ) {
-	dim = rv.count();
+enorm<sq_T>::enorm ( RV &rv ) :eEF(), mu ( rv.count() ),R(rv.count()),_iR(rv.count()),cached ( false ),dim ( rv.count() ) {};
+
+template<class sq_T>
+void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
+//Fixme test dimensions of mu0 and R0;
 	mu = mu0;
 	R = R0;
-};
-
-template<class sq_T>
-void enorm<sq_T>::dupdate( mat &v, double nu ) {
+	if(_iR.rows()!=R.rows()) _iR=R; // memory allocation!
+	R.inv ( _iR ); //update cache
+	cached=true;
+};
+
+template<class sq_T>
+void enorm<sq_T>::dupdate ( mat &v, double nu ) {
 	//
 };
 
 template<class sq_T>
-void enorm<sq_T>::tupdate( double phi, mat &vbar, double nubar ) {
+void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
 	//
 };
@@ -103,7 +137,7 @@
 template<class sq_T>
 vec enorm<sq_T>::sample() {
-	vec x( dim );
-	RNG.sample_vector( dim,x );
-	vec smp = R.sqrt_mult( x );
+	vec x ( dim );
+	RNG.sample_vector ( dim,x );
+	vec smp = R.sqrt_mult ( x );
 
 	smp += mu;
@@ -112,63 +146,79 @@
 
 template<class sq_T>
-mat enorm<sq_T>::sample( int N ) {
-	mat X( dim,N );
-	vec x( dim );
+mat enorm<sq_T>::sample ( int N ) {
+	mat X ( dim,N );
+	vec x ( dim );
 	vec pom;
 	int i;
+
 	for ( i=0;i<N;i++ ) {
-		RNG.sample_vector( dim,x );
-		pom = R.sqrt_mult( x );
+		RNG.sample_vector ( dim,x );
+		pom = R.sqrt_mult ( x );
 		pom +=mu;
-		X.set_col( i, pom);
+		X.set_col ( i, pom );
 	}
+
 	return X;
 };
 
 template<class sq_T>
-double enorm<sq_T>::eval( const vec &val ) {
-	return 0.0;	
-};
-
-
-template<class sq_T>
-enorm<sq_T>::enorm() {};
-
-template<class sq_T>
-mlnorm<sq_T>::mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R ) {
+double enorm<sq_T>::eval ( const vec &val ) {
+	return exp ( evalpdflog ( val ) );
+};
+
+template<class sq_T>
+double enorm<sq_T>::evalpdflog ( const vec &val ) {
+	if ( !cached ) {R.inv ( _iR );cached=true;}
+
+	return -0.5* ( R.cols() *0.79817986835811504957 +R.logdet() +_iR.qform ( mu-val ) );
+};
+
+
+template<class sq_T>
+mlnorm<sq_T>::mlnorm ( RV &rv,RV &rvc, mat &A, sq_T &R ) {
 	int dim = rv.count();
-	vec mu( dim );
-
-	epdf = enorm<sq_T>( rv,mu,R );
-}
-
-template<class sq_T>
-vec mlnorm<sq_T>::samplecond( vec &cond, double &lik ) {
-	this->condition( cond );
+	vec mu ( dim );
+
+	epdf = enorm<sq_T> ( rv,mu,R );
+}
+
+template<class sq_T>
+vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
+	this->condition ( cond );
 	vec smp = epdf.sample();
-	lik = epdf.eval( smp );
+	lik = epdf.eval ( smp );
 	return smp;
 }
 
 template<class sq_T>
-mat mlnorm<sq_T>::samplecond( vec &cond, vec &lik, int n ) {
+mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
 	int i;
 	int dim = rv.count();
-	mat Smp( dim,n );
-	vec smp( dim );
-	this->condition( cond );
+	mat Smp ( dim,n );
+	vec smp ( dim );
+	this->condition ( cond );
+
 	for ( i=0; i<dim; i++ ) {
 		smp = epdf.sample();
-		lik( i ) = epdf.eval( smp );
-		Smp.set_col( i ,smp );
+		lik ( i ) = epdf.eval ( smp );
+		Smp.set_col ( i ,smp );
 	}
+
 	return Smp;
 }
 
 template<class sq_T>
-void mlnorm<sq_T>::condition( vec &cond ) {
+void mlnorm<sq_T>::condition ( vec &cond ) {
 	epdf.mu = A*cond;
 //R is already assigned;
 }
 
+///////////
+
+//#ifdef HAVE_EXTERN_TEMPLATE
+
+//  extern template class enorm<fsqmat>;
+//  extern template class enorm<ldmat>;
+//#endif
+
 #endif //EF_H
Index: bdm/stat/libBM.h
===================================================================
--- bdm/stat/libBM.h (revision 27)
+++ bdm/stat/libBM.h (revision 28)
@@ -84,23 +84,9 @@
 	//! access function
 	int _dimy()const{return dimy;}
+
+		//! Destructor for future use;
+		virtual ~fnc(){};
 };
 
-//! Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
-class BM {
-public:
-	//!Logarithm of marginalized data likelihood.
-	double ll;
-
-	//!Default constructor
-	BM(){ll=0;};
-	
-	/*! \brief Incremental Bayes rule
-	@param dt vector of input data
-	@param evall If true, the filter will compute likelihood of the data record and store it in \c ll
-	*/
-	virtual void bayes ( const vec &dt, bool evall=true ) = 0;
-	//! Batch Bayes rule (columns of Dt are observations)
-	void bayes ( mat Dt );
-};
 
 //! Probability density function with numerical statistics, e.g. posterior density.
@@ -108,4 +94,6 @@
 	RV rv;
 public:
+	//!default constructor
+	epdf():rv(ivec(0)){};
 	//! Returns the required moment of the epdf
 //	virtual vec moment ( const int order = 1 );
@@ -113,8 +101,10 @@
 	virtual vec sample ()=0;
 	//! Compute probability of argument \c val
-	virtual double eval(const vec &val)
-	{
-		return 0.0;
-	};
+	virtual double eval(const vec &val){return 0.0;};
+	//! Compute log-probability of argument \c val
+	virtual double evalpdflog(const vec &val){return 0.0;};
+
+		//! Destructor for future use;
+		virtual ~epdf(){};
 };
 
@@ -130,6 +120,9 @@
 //	virtual fnc moment ( const int order = 1 );
 	//! Returns a sample from the density conditioned on \c cond, $x \sim epdf(rv|cond)$
-	virtual vec samplecond (vec &cond, double lik){};
+	virtual vec samplecond (vec &cond, double lik){return vec(0);};
 	virtual void condition (vec &cond){};
+
+		//! Destructor for future use;
+		virtual ~mpdf(){};
 };
 
@@ -164,6 +157,32 @@
 	//! Moves from $t$ to $t+1$, i.e. perfroms the actions and reads response of the system. 
 	void step();
+		
 };
 
+/*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
+	
+*/
+class BM {
+public:
+	//!Logarithm of marginalized data likelihood.
+	double ll;
+	//!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save time.
+	bool evalll;
+
+	//!Default constructor
+	BM():ll(0),evalll(true){};
+	
+	/*! \brief Incremental Bayes rule
+	@param dt vector of input data
+	*/
+	virtual void bayes ( const vec &dt) = 0;
+	//! Batch Bayes rule (columns of Dt are observations)
+	void bayes ( mat Dt );
+	//! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
+	virtual epdf* _epdf()=0;
+
+		//! Destructor for future use;
+		virtual ~BM(){};
+};
 
 #endif // BM_H
