Index: bdm/stat/libEF.h
===================================================================
--- bdm/stat/libEF.h (revision 176)
+++ bdm/stat/libEF.h (revision 178)
@@ -44,7 +44,7 @@
 	virtual double lognc() const =0;
 	//!TODO decide if it is really needed
-	virtual void dupdate ( mat &v ) {it_error ( "Not implemneted" );};
+	virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
 	//!Evaluate normalized log-probability
-	virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemneted" );return 0.0;};
+	virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
 	//!Evaluate normalized log-probability
 	virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();}
@@ -90,7 +90,12 @@
 	//original Bayes
 	void bayes ( const vec &dt );
-	//!Flatten the posterior
-	virtual void flatten ( BMEF * B ) {it_error ( "Not implemented" );}
-};
+	//!Flatten the posterior according to the given BMEF (of the same type!)
+	virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
+	//!Flatten the posterior as if to keep nu0 data
+	virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
+};
+
+template<class sq_T>
+class mlnorm;
 
 /*!
@@ -100,5 +105,4 @@
 */
 template<class sq_T>
-
 class enorm : public eEF {
 protected:
@@ -110,11 +114,10 @@
 	int dim;
 public:
-//	enorm() :eEF() {};
 	//!Default constructor
-	enorm ( RV &rv );
+	enorm ( const RV &rv );
 	//! Set mean value \c mu and covariance \c R
 	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 );
+	////! tupdate in exponential form (not really handy)
+	//void tupdate ( double phi, mat &vbar, double nubar );
 	//! dupdate in exponential form (not really handy)
 	void dupdate ( mat &v,double nu=1.0 );
@@ -124,8 +127,9 @@
 	mat sample ( int N ) const;
 	double eval ( const vec &val ) const ;
-	double evalpdflog ( const vec &val ) const;
+	double evalpdflog_nn ( const vec &val ) const;
 	double lognc () const;
 	vec mean() const {return mu;}
-
+	mlnorm<sq_T>* condition ( const RV &rvn );
+	enorm<sq_T>* marginal ( const RV &rv );
 //Access methods
 	//! returns a pointer to the internal mean value. Use with Care!
@@ -139,5 +143,5 @@
 
 	//! access method
-	mat getR () {return R.to_mat();}
+//	mat getR () {return R.to_mat();}
 };
 
@@ -215,10 +219,10 @@
 	};
 	//!access function
-	vec& _beta() {return beta;}
+	vec& _beta()  {return beta;}
 	//!Set internal parameters
-	void set_parameters(const vec &beta0){
-		if(beta0.length()!=beta.length()){
-			it_assert_debug(rv.length()==1,"Undefined");
-			rv.set_size(0,beta0.length());
+	void set_parameters ( const vec &beta0 ) {
+		if ( beta0.length() !=beta.length() ) {
+			it_assert_debug ( rv.length() ==1,"Undefined" );
+			rv.set_size ( 0,beta0.length() );
 		}
 		beta= beta0;
@@ -258,8 +262,8 @@
 		return pred.lognc()-lll;
 	}
-	void flatten ( BMEF* B ) {
-		eDirich* E=dynamic_cast<eDirich*> ( B );
+	void flatten ( const BMEF* B ) {
+		const eDirich* E=dynamic_cast<const eDirich*> ( B );
 		// sum(beta) should be equal to sum(B.beta)
-		const vec &Eb=E->_beta();
+		const vec &Eb=const_cast<eDirich*> ( E )->_beta();
 		est.pow ( sum ( beta ) /sum ( Eb ) );
 		if ( evalll ) {last_lognc=est.lognc();}
@@ -267,7 +271,7 @@
 	const epdf& _epdf() const {return est;};
 	void set_parameters ( const vec &beta0 ) {
-		est.set_parameters(beta0);
+		est.set_parameters ( beta0 );
 		rv = est._rv();
-		if(evalll){last_lognc=est.lognc();}
+		if ( evalll ) {last_lognc=est.lognc();}
 	}
 };
@@ -359,5 +363,5 @@
  \brief Normal distributed linear function with linear function of mean value;
 
- Mean value \f$mu=A*rvc\f$.
+ Mean value \f$mu=A*rvc+mu_0\f$.
 */
 template<class sq_T>
@@ -366,16 +370,17 @@
 	enorm<sq_T> epdf;
 	mat A;
+	vec mu_const;
 	vec& _mu; //cached epdf.mu;
 public:
 	//! Constructor
-	mlnorm ( RV &rv,RV &rvc );
+	mlnorm (const RV &rv, const RV &rvc );
 	//! Set \c A and \c R
-	void set_parameters ( const  mat &A, const sq_T &R );
+	void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
 	//!Generate one sample of the posterior
-	vec samplecond ( vec &cond, double &lik );
+	vec samplecond (const vec &cond, double &lik );
 	//!Generate matrix of samples of the posterior
-	mat samplecond ( vec &cond, vec &lik, int n );
+	mat samplecond (const vec &cond, vec &lik, int n );
 	//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
-	void condition ( vec &cond );
+	void condition (const vec &cond );
 };
 
@@ -453,6 +458,8 @@
 	//! Default constructor
 	eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
+	//! Set samples and weights
+	void set_parameters ( const vec &w0, const epdf* pdf0 );
 	//! Set sample
-	void set_parameters ( const vec &w0, epdf* pdf0 );
+	void set_samples ( const epdf* pdf0 );
 	//! Potentially dangerous, use with care.
 	vec& _w()  {return w;};
@@ -476,5 +483,5 @@
 
 template<class sq_T>
-enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
+enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
 
 template<class sq_T>
@@ -490,8 +497,8 @@
 };
 
-template<class sq_T>
-void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
-	//
-};
+// template<class sq_T>
+// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
+// 	//
+// };
 
 template<class sq_T>
@@ -531,7 +538,7 @@
 
 template<class sq_T>
-double enorm<sq_T>::evalpdflog ( const vec &val ) const {
+double enorm<sq_T>::evalpdflog_nn ( const vec &val ) const {
 	// 1.83787706640935 = log(2pi)
-	return  -0.5* ( +R.invqform ( mu-val ) ) - lognc();
+	return  -0.5* ( R.invqform ( mu-val ) );// - lognc();
 };
 
@@ -539,20 +546,21 @@
 inline double enorm<sq_T>::lognc () const {
 	// 1.83787706640935 = log(2pi)
-	return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
-};
-
-template<class sq_T>
-mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
+	return 0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
+};
+
+template<class sq_T>
+mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
 	ep =&epdf;
 }
 
 template<class sq_T>
-void mlnorm<sq_T>::set_parameters ( const mat &A0, const sq_T &R0 ) {
+void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
 	epdf.set_parameters ( zeros ( rv.count() ),R0 );
 	A = A0;
+	mu_const = mu0;
 }
 
 template<class sq_T>
-vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
+vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
 	this->condition ( cond );
 	vec smp = epdf.sample();
@@ -562,5 +570,5 @@
 
 template<class sq_T>
-mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
+mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
 	int i;
 	int dim = rv.count();
@@ -579,9 +587,50 @@
 
 template<class sq_T>
-void mlnorm<sq_T>::condition ( vec &cond ) {
-	_mu = A*cond;
+void mlnorm<sq_T>::condition (const vec &cond ) {
+	_mu = A*cond + mu_const;
 //R is already assigned;
 }
 
+template<class sq_T>
+enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) {
+	ivec irvn = rvn.dataind ( rv );
+
+	sq_T Rn ( R,irvn );
+	enorm<sq_T>* tmp = new enorm<sq_T>( rvn );
+	tmp->set_parameters ( mu ( irvn ), Rn );
+	return tmp;
+}
+
+template<class sq_T>
+mlnorm<sq_T>* enorm<sq_T>::condition ( const RV &rvn ) {
+
+	RV rvc = rv.subt ( rvn );
+	it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" );
+	//Permutation vector of the new R
+	ivec irvn = rvn.dataind ( rv );
+	ivec irvc = rvc.dataind ( rv );
+	ivec perm=concat ( irvn , irvc );
+	sq_T Rn ( R,perm );
+
+	//fixme - could this be done in general for all sq_T?
+	mat S=R.to_mat();
+	//fixme
+	int n=rvn.count()-1;
+	int end=R.rows()-1;
+	mat S11 = S.get ( 0,n, 0, n );
+	mat S12 = S.get ( rvn.count(), end, 0, n );
+	mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
+
+	vec mu1 = mu ( irvn );
+	vec mu2 = mu ( irvc );
+	mat A=S12*inv ( S22 );
+	sq_T R_n ( S11 - A *S12.T() );
+
+	mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc );
+
+	tmp->set_parameters ( A,mu1-A*mu2,R_n );
+	return tmp;
+}
+
 ///////////
 
