Index: /mixpp.kdevelop.filelist
===================================================================
--- /mixpp.kdevelop.filelist (revision 11)
+++ /mixpp.kdevelop.filelist (revision 14)
@@ -1,4 +1,5 @@
 # KDevelop Custom Project File List
 Makefile
+arx.cpp
 doc
 doc/latex
Index: /bdm/estim/libPF.h
===================================================================
--- /bdm/estim/libPF.h (revision 14)
+++ /bdm/estim/libPF.h (revision 14)
@@ -0,0 +1,72 @@
+/*!
+  \file
+  \brief Bayesian Filtering using stochastic sampling (Particle Filters)
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#ifndef PF_H
+#define PF_H
+
+#include <itpp/itbase.h>
+#include "libBM.h"
+#include "libDC.h"
+
+using namespace itpp;
+
+enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
+
+/*!
+* \brief A Particle Filter prototype
+
+Bayesian Filtering equations hold.
+*/
+class PF : public BM { 
+protected:
+	int n; //number of particles
+	vec w; //particle weights
+	Uniform_RNG URNG; //used for resampling
+	
+public:
+	//! Returns indexes of particles that should be resampled. The ordering MUST guarantee inplace replacement. (Important for MPF.)
+	ivec resample(RESAMPLING_METHOD method = SYSTEMATIC);
+	PF (vec w);
+	//TODO remove or implement bayes()!
+	void bayes(const vec &dt, bool evell){};
+};
+
+/*!
+* \brief Trivial particle filter with proposal density that is not conditioned on the data.
+
+
+*/
+
+class TrivialPF : public PF {
+	Array<vec> ptcls;
+	
+	bool is_proposal;
+	mpdf *prop;
+	mpdf *par;
+	mpdf *obs;
+	
+	public:
+	TrivialPF(mpdf &par, mpdf &obs, mpdf &prop, int n0);
+	TrivialPF(mpdf &par, mpdf &obs, int n0);
+	void bayes(const vec &dt, bool evalll);
+};
+
+class MPF : public TrivialPF {
+	Array<BM> Bms;
+	public:
+	MPF(BM &B, mpdf &prop, mpdf &obs, mpdf &par);
+	void bayes(vec &dt);	
+};
+
+#endif // KF_H
+
+
Index: /bdm/estim/libKF.cpp
===================================================================
--- /bdm/estim/libKF.cpp (revision 14)
+++ /bdm/estim/libKF.cpp (revision 14)
@@ -0,0 +1,51 @@
+#include <itpp/itbase.h>
+#include "libBM.h"
+#include "libKF.h"
+
+using namespace itpp;
+
+using std::endl;
+
+KalmanFull::KalmanFull( mat A0, mat B0, mat C0, mat D0, mat R0, mat Q0, mat P0, vec mu0 ) {
+	dimx = A0.rows();
+	dimu = B0.cols();
+	dimy = C0.rows();
+
+	it_assert_debug( A0.cols()==dimx, "KalmanFull: A is not square" );
+	it_assert_debug( B0.rows()==dimx, "KalmanFull: B is not compatible" );
+	it_assert_debug( C0.cols()==dimx, "KalmanFull: C is not square" );
+	it_assert_debug(( D0.rows()==dimy ) || ( D0.cols()==dimu ),	"KalmanFull: D is not compatible" );
+	it_assert_debug(( R0.cols()==dimy ) || ( R0.rows()==dimy ), "KalmanFull: R is not compatible" );
+	it_assert_debug(( Q0.cols()==dimx ) || ( Q0.rows()==dimx ), "KalmanFull: Q is not compatible" );
+
+	A = A0;
+	B = B0;
+	C = C0;
+	D = D0;
+	R = R0;
+	Q = Q0;
+	mu = mu0;
+	P = P0;
+}
+
+void KalmanFull::bayes( const vec &dt, bool evalll ) {
+	it_assert_debug( dt.length()==( dimy+dimu ),"KalmanFull::bayes wrong size of dt" );
+
+	vec u = dt.get( dimy,dimy+dimu-1 );
+	vec y = dt.get( 0,dimy-1 );
+	//Time update
+	mu = A*mu + B*u;
+	P  = A*P*A.transpose() + Q;
+
+	//Data update
+	_Ry = C*P*C.transpose() + R;
+	_iRy = inv( _Ry );
+	_K = P*C.transpose()*_iRy;
+	P -= _K*C*P; // P = P -KCP;
+	mu += _K*( y-C*mu-D*u );
+};
+
+std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf ) {
+	os << "mu:" << kf.mu << endl << "P:" << kf.P <<endl;
+}
+
Index: /bdm/estim/arx.cpp
===================================================================
--- /bdm/estim/arx.cpp (revision 14)
+++ /bdm/estim/arx.cpp (revision 14)
@@ -0,0 +1,6 @@
+#include <itpp/itbase.h>
+#include "libBM.h"
+
+using namespace itpp;
+
+
Index: /bdm/estim/libKF.h
===================================================================
--- /bdm/estim/libKF.h (revision 14)
+++ /bdm/estim/libKF.h (revision 14)
@@ -0,0 +1,145 @@
+/*!
+  \file
+  \brief Bayesian Filtering for linear Gaussian models (Kalman Filter) and extensions
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#ifndef KF_H
+#define KF_H
+
+#include <itpp/itbase.h>
+#include "libBM.h"
+#include "libDC.h"
+
+
+using namespace itpp;
+
+/*!
+* \brief Basic Kalman filter with full matrices (education purpose only)! Will be deleted soon!
+*/
+class KalmanFull : public BM { 
+	int dimx, dimy, dimu;
+	mat A, B, C, D, R, Q;
+	
+	//cache 
+	mat _Pp, _Ry, _iRy, _K;
+public:
+	//posterior 
+	//! Mean value of the posterior density
+	vec mu;
+	//! Variance of the posterior density
+	mat P;
+
+public:
+	//! Full constructor 
+	KalmanFull ( mat A, mat B, mat C, mat D, mat R, mat Q, mat P0, vec mu0);
+	//! Here dt = [yt;ut] of appropriate dimensions
+	void bayes(const vec &dt, bool evalll=true); 
+
+	friend std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf );
+
+};
+
+
+/*!
+* \brief Kalman filter with covaraince matrices in square root form.
+*/
+template<class sq_T>
+class Kalman : public BM { 
+	int dimx, dimy, dimu;
+	mat A, B, C, D;
+	sq_T R, Q;
+	
+	//cache
+	mat _K;
+	vec _yp;
+	sq_T _Ry,_iRy;
+public:
+	//posterior 
+	//! Mean value of the posterior density
+	vec mu;
+	//! Mean value of the posterior density
+	sq_T P;
+
+public:
+	//! Full constructor 
+	Kalman ( mat A0, mat B0, mat C0, mat D0, sq_T R0, sq_T Q0, sq_T P0, vec mu0 );
+	//! Here dt = [yt;ut] of appropriate dimensions
+	void bayes(const vec &dt, bool evalll=true); 
+
+	friend std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf );
+
+};
+
+//////// INstance
+
+template<class sq_T>
+Kalman<sq_T>::Kalman( mat A0, mat B0, mat C0, mat D0, sq_T R0, sq_T Q0, sq_T P0, vec mu0 ) {
+	dimx = A0.rows();
+	dimu = B0.cols();
+	dimy = C0.rows();
+
+	it_assert_debug( A0.cols()==dimx, "Kalman: A is not square" );
+	it_assert_debug( B0.rows()==dimx, "Kalman: B is not compatible" );
+	it_assert_debug( C0.cols()==dimx, "Kalman: C is not square" );
+	it_assert_debug(( D0.rows()==dimy ) || ( D0.cols()==dimu ),	"Kalman: D is not compatible" );
+	it_assert_debug(( R0.cols()==dimy ) || ( R0.rows()==dimy ), "Kalman: R is not compatible" );
+	it_assert_debug(( Q0.cols()==dimx ) || ( Q0.rows()==dimx ), "Kalman: Q is not compatible" );
+
+	A = A0;
+	B = B0;
+	C = C0;
+	D = D0;
+	R = R0;
+	Q = Q0;
+	mu = mu0;
+	P = P0;
+
+	ll = 0;
+//Fixme should we assign cache??
+	_iRy = eye(dimy); // needed in inv(_iRy)
+}
+
+template<class sq_T>
+void Kalman<sq_T>::bayes( const vec &dt , bool evalll) {
+	it_assert_debug( dt.length()==( dimy+dimu ),"KalmanFull::bayes wrong size of dt" );
+
+	vec u = dt.get( dimy,dimy+dimu-1 );
+	vec y = dt.get( 0,dimy-1 );
+	//Time update
+	mu = A*mu + B*u;
+	//P  = A*P*A.transpose() + Q; in sq_T
+	P.mult_sym( A );
+	P+=Q;
+
+	//Data update
+	//_Ry = C*P*C.transpose() + R; in sq_T
+	_Ry.mult_sym( C, P);
+	_Ry+=R;
+
+	mat Pfull = P.to_mat();
+	
+	_Ry.inv( _iRy ); // result is in _iRy;
+	_K = Pfull*C.transpose()*(_iRy.to_mat());
+	P -= _K*C*Pfull; // P = P -KCP;
+	_yp = y-C*mu-D*u; //y prediction
+	mu += _K*( _yp );
+	
+	if (evalll==true) {
+	ll+= -0.5*(_Ry.cols()*0.79817986835811504957 \
+	+_Ry.logdet() +_iRy.qform(_yp));
+	}
+};
+
+//extern template class Kalman<ldmat>; 
+
+
+#endif // KF_H
+
+
Index: /bdm/estim/libPF.cpp
===================================================================
--- /bdm/estim/libPF.cpp (revision 14)
+++ /bdm/estim/libPF.cpp (revision 14)
@@ -0,0 +1,123 @@
+#include "libPF.h"
+
+using namespace itpp;
+
+using std::endl;
+
+PF::PF(vec w0){w=w0/sum(w0); n=w.length();}
+
+ivec PF::resample( RESAMPLING_METHOD method ) {
+	ivec ind=zeros_i( n );
+	ivec N_babies = zeros_i( n );
+	vec cumDist = cumsum( w );
+	vec u( n );
+	int i,j,parent;
+	double u0;
+
+	switch ( method ) {
+	case MULTINOMIAL:
+		u( n - 1 ) = pow( URNG.sample(), 1.0 / n );
+		for ( i = n - 2;i >= 0;i-- ) {
+			u( i ) = u( i + 1 ) * pow( URNG.sample(), 1.0 / ( i + 1 ) );
+		}
+		break;
+	case STRATIFIED:
+		for ( i = 0;i < n;i++ ) {
+			u( i ) = ( i + URNG.sample() ) / n;
+		}
+		break;
+	case SYSTEMATIC:
+		u0 = URNG.sample();
+		for ( i = 0;i < n;i++ ) {
+			u( i ) = ( i + u0 ) / n;
+		}
+		break;
+	default:
+		it_error( "PF::resample(): Unknown resampling method" );
+	}
+	
+// 	// <<<<<<<<<<<
+// 	using std::cout;
+// 	cout << "resampling:" << endl;
+// 	cout << w <<endl;
+// 	cout << cumDist <<endl;
+// 	cout << u <<endl;
+
+	// U is now full
+	j = 0;
+	for ( i = 0;i < n;i++ ) {
+		while ( u( i ) > cumDist( j ) ) j++;
+		N_babies( j )++;
+	}
+
+// 	// <<<<<<<<<<<
+// 	cout << N_babies <<endl;
+
+	// We have assigned new babies for each Particle
+	// Now, we fill the resulting index such that:
+	// * particles with at least one baby should not move *
+	// This assures that reassignment can be done inplace;
+
+	// find the first parent;
+	parent=0; while(N_babies(parent)==0) parent++;
+	// Build index
+	for ( i = 0;i < n;i++ ) {
+		if ( N_babies( i ) > 0 ){
+			ind( i ) = i;
+			N_babies( i ) --; //this index was now replicated;
+			}
+		else {
+			// test if the parent has been fully replicated 
+			// if yes, find the next one
+			while((N_babies(parent)==0) || (N_babies(parent)==1 && parent>i) ) parent++;
+			// Replicate parent
+			ind (i) = parent;
+			N_babies( parent ) --; //this index was now replicated;
+		}
+/*			// <<<<<<<<<<<
+			cout <<ind << endl;*/
+	}
+	return ind;
+}
+
+TrivialPF::TrivialPF( mpdf &par0, mpdf &obs0, mpdf &prop0, int n0 ) : PF(1/n*ones(n)) {
+	is_proposal = true;
+	prop = &prop0;
+	par = &par0;
+	obs = &obs0;
+}
+TrivialPF::TrivialPF( mpdf &par0, mpdf &obs0,  int n0 ) : PF(1/n*ones(n)) {
+	is_proposal = false;
+	par = &par0;
+	obs = &obs0;
+}
+
+void TrivialPF::bayes( const vec &dt , bool evalll ) {
+	int i;
+	vec oldp;
+	double ll, gl, sum = 0.0;
+	Sort<double> S;
+	ivec ind, iw;
+	/*
+	//generate new samples
+	for ( i=0;i<n;i++ ) {
+		prop->evalcond( ptcls( i ), &prop_cond );
+		ptcls( i ) = prop_cond.sample();
+		gl = prop_cond.eval( ptcls( i ) );
+
+		obs.evalcond( ptcls( i ), &obs_cond );
+		ll = obs_cond.eval( dt );
+		w( i ) *= ll/gl;
+	}
+	//renormalize
+	for ( i=0;i<n;i++ ){sum+=w( i );};
+	w( i ) /=sum; //?
+	//
+	ind = resample();
+	iw = S.sort_index( 0,n-1,w ); // the first one in iw is the strongest
+
+	for ( i=0;i<n;i++ ) {
+		ptcls( i ) = ptcls( i ); //potentionally dangerous!
+	}
+	*/
+}
Index: /bdm/math/libDC.h
===================================================================
--- /bdm/math/libDC.h (revision 14)
+++ /bdm/math/libDC.h (revision 14)
@@ -0,0 +1,175 @@
+/*!
+ * \file
+ * \brief Matrices in decomposed forms (LDL', LU, UDU', etc).
+ * \author Vaclav Smidl.
+ *
+ * -----------------------------------
+ * BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+ *
+ * Using IT++ for numerical operations
+ * -----------------------------------
+ */
+
+#ifndef DC_H
+#define DC_H
+
+#include <itpp/itbase.h>
+
+using namespace itpp;
+
+/*! \brief Virtual class for representation of double symmetric matrices in square-root form.
+
+All operations defined on this class should be optimized for the chosed decomposition.
+*/
+class sqmat {
+public:
+	/*!
+	 * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$.
+	 * @param  v Vector forming the outer product to be added
+	 * @param  w weight of updating; can be negative
+
+	 BLAS-2b operation.
+	 */
+	virtual void opupdt( const vec &v, double w ) =0;
+
+	/*! \brief Conversion to full matrix.
+	*/
+
+	virtual mat to_mat() =0;
+
+	/*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$
+	@param C multiplying matrix,
+	@param trans if true, product $V = C'*V*C$ will be computed instead;
+	*/
+	virtual void mult_sym( const mat &C, bool trans=true ) =0;
+
+	
+	/*!
+	\brief Logarithm of a determinant.
+	
+	*/
+	virtual double logdet() =0;
+
+	/*!
+	\brief Multiplies square root of $V$ by vector $x$.
+	
+	Used e.g. in generating normal samples.
+	*/
+	virtual vec sqrt_mult(vec &v) =0;
+	
+	/*!
+	\brief Evaluates quadratic form $x= v'*V*v$;
+	
+	*/
+	virtual double qform(vec &v) =0;
+
+//	//! easy version of the 
+//	sqmat inv();
+
+	//! Clearing matrix so that it corresponds to zeros.
+	virtual void clear() =0;
+	
+	//! Reimplementing common functions of mat: cols().
+	virtual int cols() =0;
+
+	//! Reimplementing common functions of mat: cols().
+	virtual int rows() =0;
+
+protected:
+	int dim;
+};
+
+
+/*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix.
+
+This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices;
+*/
+class fsqmat: sqmat {
+	void opupdt( const vec &v, double w );
+	mat to_mat();
+	void mult_sym( const mat &C, bool trans=false );
+   void mult_sym( const mat &C, fsqmat &U, bool trans=false );
+	void inv(fsqmat &Inv);
+	void clear();
+	
+	
+	//! Constructor
+	fsqmat(const mat &M);
+	
+	/*! \brief Matrix inversion preserving the chosen form.
+
+	@param Inv a space where the inverse is stored.
+	
+	*/
+	virtual void inv(fsqmat* Inv);
+	
+
+};
+
+class ldmat: sqmat {
+public:
+
+	//! Construct by copy of L and D.
+	ldmat( const mat &L, const vec &D );
+	//! Construct by decomposition of full matrix V.
+	ldmat( mat V );
+	ldmat ();
+
+	// Reimplementation of compulsory operatios
+
+	void opupdt( const vec &v, double w );
+	mat to_mat();
+   void mult_sym( const mat &C, bool trans=false );
+	void add ( const ldmat &ld2, double w=1.0 );
+	double logdet();
+	double qform(vec &v);
+//	sqmat& operator -= ( const sqmat & ld2 );
+	void clear();
+	int cols();
+	int rows();
+	vec sqrt_mult(vec &v);
+
+	/*! \brief Matrix inversion preserving the chosen form.
+
+	@param Inv a space where the inverse is stored.
+	
+	*/
+	virtual void inv(ldmat &Inv);
+	
+	/*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class.
+
+	@param Inv a space where the inverse is stored.
+	
+	*/
+   void mult_sym( const mat &C, ldmat &U, bool trans=false );
+
+	/*! \brief Transforms general $A'D0A$ into pure $L'DL$
+	
+	The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$
+
+	@param A general matrix
+	@param D0 general vector
+	
+	*/
+   void ldform( mat &A, vec &D0 );
+
+	ldmat& operator += (const ldmat &ldA);
+	ldmat& operator -= (const ldmat &ldA);
+	ldmat& operator *= (double x);
+	
+	friend std::ostream &operator<< ( std::ostream &os, ldmat &sq );
+
+protected:
+	vec D;
+	mat L;
+
+};
+
+//////// Operations:
+
+inline ldmat& ldmat::operator += (const ldmat &ldA)  {this->add(ldA);return *this;}
+inline ldmat& ldmat::operator -= (const ldmat &ldA)  {this->add(ldA,-1.0);return *this;}
+inline int ldmat::cols(){return L.cols();}
+inline int ldmat::rows(){return L.rows();}
+
+#endif // DC_H
Index: /bdm/math/libDC.cpp
===================================================================
--- /bdm/math/libDC.cpp (revision 14)
+++ /bdm/math/libDC.cpp (revision 14)
@@ -0,0 +1,361 @@
+#include <itpp/itbase.h>
+#include "libDC.h"
+
+using namespace itpp;
+
+using std::endl;
+
+//! Auxiliary function dydr; dyadic reduction
+void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx );
+
+//! Auxiliary function ltuinv; inversion of a triangular matrix;
+//TODO can be done via: dtrtri.f from lapack
+mat ltuinv( const mat &L );
+
+
+
+ldmat::ldmat( const mat &exL, const vec &exD ) {
+	D = exD;
+	L = exL;
+	dim = exD.length();
+}
+
+ldmat::ldmat() {
+	vec D ;
+	mat L;
+	dim = 0;
+}
+
+ldmat::ldmat( const mat V ) {
+//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!" );
+
+	
+}
+
+void ldmat::opupdt( const vec &v,  double w ) {
+	int dim = D.length();
+	double kr;
+	vec r = v;
+	//beware! it is potentionally dangerous, if ITpp change _behaviour of _data()!
+	double *Lraw = L._data();
+	double *Draw = D._data();
+	double *rraw = r._data();
+
+	it_assert_debug( v.length() == dim, "LD::ldupdt vector v is not compatible with this ld." );
+
+	for ( int i = dim - 1; i >= 0; i-- ) {
+		dydr( rraw, Lraw + i, &w, Draw + i, rraw + i, 0, i, &kr, 1, dim );
+	}
+}
+
+std::ostream &operator<< ( std::ostream &os,  ldmat &ld ) {
+	os << "L:" << ld.L << endl;
+	os << "D:" << ld.D << endl;
+}
+
+mat ldmat::to_mat() {
+	int dim = D.length();
+	mat V( dim, dim );
+	double sum;
+	int r, c, cc;
+
+	for ( r = 0;r < dim;r++ ) { //row cycle
+		for ( c = r;c < dim;c++ ) {
+			//column cycle, using symmetricity => c=r!
+			sum = 0.0;
+			for ( cc = c;cc < dim;cc++ ) { //cycle over the remaining part of the vector
+				sum += L( cc, r ) * D( cc ) * L( cc, c );
+				//here L(cc,r) = L(r,cc)';
+			}
+			V( r, c ) = sum;
+			// symmetricity
+			if ( r != c ) {V( c, r ) = sum;};
+		}
+	}
+	mat V2 = L.transpose()*diag( D )*L;
+	return V2;
+}
+
+
+void ldmat::add( const ldmat &ld2, double w ) {
+	int dim = D.length();
+
+	it_assert_debug( ld2.D.length() == dim, "LD.add() incompatible sizes of LDs;" );
+
+	//Fixme can be done more efficiently either via dydr or ldform
+	for ( int r = 0; r < dim; r++ ) {
+		// Add columns of ld2.L' (i.e. rows of ld2.L) as dyads weighted by ld2.D
+		this->opupdt( ld2.L.get_row( r ), w*ld2.D( r ) );
+	}
+}
+
+void ldmat::clear(){L.clear(); for ( int i=0;i<L.cols();i++ ){L( i,i )=1;}; D.clear();}
+
+void ldmat::inv( ldmat &Inv ) {
+	int dim = D.length();
+	Inv.clear();   //Inv = zero in LD
+	mat U = ltuinv( L );
+
+	//Fixme can be done more efficiently either via dydr or ldform
+	for ( int r = 0; r < dim; r++ ) {
+		// Add columns of U as dyads weighted by 1/D
+		Inv.opupdt( U.get_col( r ), 1.0 / D( r ) );
+	}
+}
+
+void ldmat::mult_sym( const mat &C, bool trans ) {
+
+//TODO better
+
+	it_assert_debug( C.cols()==L.cols(), "ldmat::mult_sym wrong input argument" );
+	mat Ct=C;
+
+	if ( trans==false ) { // return C*this*C'
+		Ct *= this->to_mat();
+		Ct *= C.transpose();
+	} else {	// return C'*this*C
+		Ct = C.transpose();
+		Ct *= this->to_mat();
+		Ct *= C;
+	}
+
+	ldmat Lnew=ldmat( Ct );
+	L = Lnew.L;
+	D = Lnew.D;
+}
+
+void ldmat::mult_sym( const mat &C, ldmat &U, bool trans ) {
+
+//TODO better
+
+//TODO input test
+
+	mat Ct=C;
+
+	if ( trans==false ) { // return C*this*C'
+		Ct *= U.to_mat();
+		Ct *= C.transpose();
+	} else {	// return C'*this*C
+		Ct = C.transpose();
+		Ct *= U.to_mat();
+		Ct *= C;
+	}
+
+	ldmat Lnew=ldmat( Ct );
+	L = Lnew.L;
+	D = Lnew.D;
+}
+
+double ldmat::logdet() {
+	double ldet = 0.0;
+	int i;
+// sum logarithms of diagobal elements
+	for ( i=0; i<D.length(); i++ ){ldet+=log( D( i ) );};
+}
+
+double ldmat::qform( vec &v ) {
+	double x = 0.0, sum;
+	int i,j;
+
+	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 );}
+		x +=D( i )*sum*sum;
+	};
+	return x;
+}
+
+ldmat& ldmat::operator *= ( double x ) {
+	int i;
+	for ( i=0;i<D.length();i++ ){D( i )*=x;};
+}
+
+vec ldmat::sqrt_mult( vec &x ) {
+	int i,j;
+	vec res( dim );
+	double sum;
+	for ( i=0;i<dim;i++ ) {//for each element of result
+		res( i ) = 0.0;
+		for ( j=i;j<dim;j++ ) {//sum D(j)*L(:,i).*x
+			res( i ) += sqrt( D( j ) )*L( j,i )*x( j );
+		}
+	}
+	vec res2 = L.transpose()*diag( sqrt( D ) )*x;
+	return res2;
+}
+
+void ldmat::ldform( mat &A,vec &D0 ) {
+	int m = A.rows();
+	int n = A.cols();
+	int mn = (m<n) ? m :n ;
+
+	it_assert_debug( A.cols()==dim,"ldmat::ldform A is not compatible" );
+	it_assert_debug( D.length()==A.rows(),"ldmat::ldform Vector D must have the length as row count of A" );
+
+	L=concat_vertical( zeros( n,n ), diag( sqrt( D0 ) )*A );
+	D=zeros( n+m );
+	
+	//unnecessary big L and D will be made smaller at the end of file
+	
+	vec w=zeros( n );
+	vec v;
+	double sum, beta, pom;
+
+	int cc=0;
+	int i=n; // +1 in .m
+	int ii,j,jj;
+	while (( i> ( n-mn+1-cc ) )&&( i>1 ) ) {
+		i--;
+		sum = 0.0;
+		v.set_size( m+i-( n-cc ) ); //prepare v
+		for ( ii=n-cc;ii<m+i;i++ ) {
+			sum+= L( ii,i )*L( ii,i );
+			v( ii )=L( ii,i ); //assign v
+		}
+		if ( L( m+i,i )==0 ) {
+			beta = sqrt( sum );
+		} else {
+			beta = L( m+i,i )+sign( L( m+i,i ) )*sqrt( sum ) ;
+		}
+		if ( std::fabs( beta )<eps ) {
+			cc++;
+			L.set_row( n+1-cc, L.get_row( m+i ) );
+			L.set_row( m+i,0 );
+			D( m+i )=0; L( m+i,i )=1;
+			L.set_submatrix( n+1-cc,m+i,i,i,0 );
+			continue;
+		}
+
+		sum-=v( v.length()-1 )*v( v.length()-1 ); //
+		sum/=beta*beta;
+		sum++;
+
+		v/=beta;
+		v( v.length()-1 )=1;
+
+		pom=-2/sum;
+		for ( j=i;i>=0;i-- ) {
+			w( j )=0.0;
+			for ( ii=n-cc;ii<m+i;ii++ ) {
+				w( j )+= v( ii )*L( ii,j );
+			}
+			w( j )*=pom;
+		}
+
+		for ( ii=n-cc;ii<m+i;ii++ ) {
+			for ( jj=0;jj<i-1;jj++ ) {
+				L( ii,jj )+= v( ii )*w( jj );
+			}
+		}
+		for ( ii=n-cc;ii<m+i;ii++ ) {
+			L( ii,i )= 0;
+		}
+		L( m+i,i )+=w( i );
+
+		D( m+i )=L( m+i,i )*L( m+i,i );
+		for ( ii=0;ii<i;ii++ ) {
+			L( m+i,ii )/=L( m+i,i );
+		}
+	}
+	if ( i>0 ) {
+		for ( ii=0;ii<i-1;ii++ ) {
+			jj = D.length()-1-n+ii;
+			L.set_row(jj,0);
+			L(jj,jj)=1;
+		}
+	}
+
+	//cut-out L and D;
+	L.del_rows(0,m-1);
+	D.del(0,m-1);
+}
+
+//////// Auxiliary Functions
+
+mat ltuinv( const mat &L ) {
+	int dim = L.cols();
+	mat Il = eye( dim );
+	int i, j, k, m;
+	double s;
+
+//Fixme blind transcription of ltuinv.m
+	for ( k = 1; k < ( dim );k++ ) {
+		for ( i = 0; i < ( dim - k );i++ ) {
+			j = i + k; //change in .m 1+1=2, here 0+0+1=1
+			s = L( j, i );
+			for ( m = i + 1; m < ( j - 1 ); m++ ) {
+				s += L( m, i ) * Il( j, m );
+			}
+			Il( j, i ) = -s;
+		}
+	}
+
+	return Il;
+}
+
+void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx )
+/********************************************************************
+
+   dydr = dyadic reduction, performs transformation of sum of
+          2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed
+          by R is zeroed. This version allows Dr to be NEGATIVE. Hence the name negdydr or dydr_withneg.
+
+   Parameters :
+     r ... pointer to reduced dyad
+     f ... pointer to reducing dyad
+     Dr .. pointer to the weight of reduced dyad
+     Df .. pointer to the weight of reducing dyad
+     R ... pointer to the element of r, which is to be reduced to
+           zero; the corresponding element of f is assumed to be 1.
+     jl .. lower index of the range within which the dyads are
+           modified
+     ju .. upper index of the range within which the dyads are
+           modified
+     kr .. pointer to the coefficient used in the transformation of r
+           rnew = r + kr*f
+     m  .. number of rows of modified matrix (part of which is r)
+  Remark : Constant mzero means machine zero and should be modified
+           according to the precision of particular machine
+
+                                                 V. Peterka 17-7-89
+
+  Added:
+     mx .. number of rows of modified matrix (part of which is f)  -PN
+
+********************************************************************/
+{
+	int j, jm;
+	double kD, r0;
+	double mzero = 2.2e-16;
+	double threshold = 1e-4;
+
+	if ( fabs( *Dr ) < mzero ) *Dr = 0;
+	r0 = *R;
+	*R = 0.0;
+	kD = *Df;
+	*kr = r0 * *Dr;
+	*Df = kD + r0 * ( *kr );
+	if ( *Df > mzero ) {
+		kD /= *Df;
+		*kr /= *Df;
+	} else {
+		kD = 1.0;
+		*kr = 0.0;
+		if ( *Df < -threshold ) {
+			it_warning( "Problem in dydr: subraction of dyad results in negative definitness. Likely mistake in calling function." );
+		}
+		*Df = 0.0;
+	}
+	*Dr *= kD;
+	jm = mx * jl;
+	for ( j = m * jl; j < m*jh; j += m ) {
+		r[j] -=  r0 * f[jm];
+		f[jm] += *kr * r[j];
+		jm += mx;
+	}
+}
+
+
Index: /bdm/stat/libFN.h
===================================================================
--- /bdm/stat/libFN.h (revision 14)
+++ /bdm/stat/libFN.h (revision 14)
@@ -0,0 +1,27 @@
+//
+// C++ Interface: itpp_ext
+//
+// Description: 
+//
+//
+// Author: smidl <smidl@utia.cas.cz>, (C) 2008
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include <itpp/itbase.h>
+#include "libBM"
+
+using namespace itpp;
+ 
+class constfn : public fnc {
+	RV rv;
+	ivec indexlist; // needed by evalsome
+	vec val;
+
+public:
+	vec eval() {return val;};
+	vec eval(vec &cond) {return val;};
+	
+	vec evalsome(ivec &rvind);
+}
Index: /bdm/stat/libEF.cpp
===================================================================
--- /bdm/stat/libEF.cpp (revision 14)
+++ /bdm/stat/libEF.cpp (revision 14)
@@ -0,0 +1,7 @@
+#include <itpp/itbase.h>
+#include "libEF.h"
+
+using namespace itpp;
+
+using std::cout;
+
Index: /bdm/stat/libBM.cpp
===================================================================
--- /bdm/stat/libBM.cpp (revision 14)
+++ /bdm/stat/libBM.cpp (revision 14)
@@ -0,0 +1,69 @@
+#include <itpp/itbase.h>
+#include "libBM.h"
+#include "itpp_ext.h"
+
+using namespace itpp;
+
+using std::cout;
+
+void RV::init( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs ) {
+	//
+	int i;
+	len = in_ids.length();
+	//PRUDENT_MODE
+	// All vetors should be of same length
+	if ( 	 ( len != in_ids.length() ) || \
+	        ( len != in_names.length() ) || \
+	        ( len != in_sizes.length() ) || \
+	        ( len != in_times.length() ) || \
+	        ( len != in_obs.length() ) ) {
+		it_error( "RV::RV inconsistent length of input vectors." );
+	}
+
+	ids = in_ids;
+	names = in_names;
+	sizes = in_sizes;
+	times = in_times;
+	obs = in_obs;
+	size = 0;
+	for(i=0;i<len;i++){size+=sizes(i);}
+};
+
+RV::RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs ) {
+	init ( in_ids, in_names, in_sizes, in_times, in_obs );
+}
+
+RV::RV () {};
+
+RV::RV ( ivec in_ids ) {
+	
+	len = in_ids.length();
+	Array<std::string> A( len );
+	std::string rvstr = "rv";
+
+	for ( int i = 0; i < len;i++ ) {
+		A( i ) = rvstr + to_str(i);
+	}
+
+	init ( in_ids, A, ones_i( len ), zeros_i( len ), zeros_i( len ) );
+}
+
+RV RV::subselect(ivec ind){
+	return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind), obs(ind));
+}
+
+void RV::t(int delta){ times +=delta;}
+
+RV RV::operator()(ivec ind){
+	return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind), obs(ind));
+}
+
+std::ostream &operator<<( std::ostream &os, const RV &rv ) {
+
+	for ( int i = 0; i < rv.len ;i++ ) {
+		os << rv.ids( i ) << "(" << rv.sizes( i ) << ")" <<  // id(size)=
+		"=" << rv.names( i )  << "_{"  << rv.times( i ) << "}; "; //name_{time}
+	}
+	return os;
+}
+
Index: /bdm/stat/libEF.h
===================================================================
--- /bdm/stat/libEF.h (revision 14)
+++ /bdm/stat/libEF.h (revision 14)
@@ -0,0 +1,168 @@
+/*!
+  \file
+  \brief Probability distributions for Exponential Family models.
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#ifndef EF_H
+#define EF_H
+
+#include <itpp/itbase.h>
+#include "libDC.h"
+#include "libBM.h"
+//#include <std>
+
+using namespace itpp;
+
+/*!
+* \brief General conjugate exponential family posterior density.
+
+* More?...
+*/
+class eEF : public epdf {
+
+public:
+	virtual void tupdate( double phi, mat &vbar, double nubar ) {};
+	virtual void dupdate( mat &v,double nu=1.0 ) {};
+};
+
+class mEF : public mpdf {
+
+public:
+
+};
+
+/*!
+* \brief General exponential family density
+
+* More?...
+*/
+template<class sq_T>
+class enorm : public eEF {
+	int dim;
+	vec mu;
+	sq_T R;
+public:
+	enorm( RV &rv, vec &mu, sq_T &R );
+	enorm();
+	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;
+};
+
+/*!
+ \brief
+*/
+template<class sq_T>
+class mlnorm : public mEF {
+	enorm<sq_T> epdf;
+	mat A;
+public:
+	//! Constructor
+	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 );
+};
+
+////////////////////////
+
+template<class sq_T>
+enorm<sq_T>::enorm( RV &rv, vec &mu0, sq_T &R0 ) {
+	dim = rv.count();
+	mu = mu0;
+	R = R0;
+};
+
+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 ) {
+	//
+};
+
+template<class sq_T>
+vec enorm<sq_T>::sample() {
+	vec x( dim );
+	RNG.sample_vector( dim,x );
+	vec smp = R.sqrt_mult( x );
+
+	smp += mu;
+	return smp;
+};
+
+template<class sq_T>
+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 );
+		pom +=mu;
+		X.set_col( i, pom);
+	}
+	return X;
+};
+
+template<class sq_T>
+double enorm<sq_T>::eval( const vec &val ) {
+	//
+};
+
+
+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 ) {
+	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 smp = epdf.sample();
+	lik = epdf.eval( smp );
+	return smp;
+}
+
+template<class sq_T>
+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 );
+	for ( i=0; i<dim; i++ ) {
+		smp = epdf.sample();
+		lik( i ) = epdf.eval( smp );
+		Smp.set_col( i ,smp );
+	}
+	return Smp;
+}
+
+template<class sq_T>
+void mlnorm<sq_T>::condition( vec &cond ) {
+	epdf.mu = A*cond;
+//R is already assigned;
+}
+
+#endif //EF_H
Index: /bdm/stat/libBM.h
===================================================================
--- /bdm/stat/libBM.h (revision 14)
+++ /bdm/stat/libBM.h (revision 14)
@@ -0,0 +1,116 @@
+/*!
+  \file
+  \brief Bayesian Models (bm) that use Bayes rule to learn from observations
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#ifndef BM_H
+#define BM_H
+
+#include <itpp/itbase.h>
+//#include <std>
+
+using namespace itpp;
+
+/*!
+* \brief Class representing variables, most often random variables
+
+* More?...
+*/
+class RV {
+	int size,len;
+	ivec ids;
+	ivec sizes;
+	ivec times;
+	ivec obs;
+	Array<std::string> names;
+
+private:
+	void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs );
+public:
+	//! Full constructor which is called by the others
+	RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs );
+	//! default constructor
+	RV ( ivec ids );
+	//! Empty constructor will be set later
+	RV ();
+	
+	//! Printing output e.g. for debugging.
+	friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
+
+	//! Return length (number of scalars) of the RV.
+	int count() {return size;}
+	//TODO why not inline and later??
+	
+	//! Find indexes of another rv in self
+	ivec find(RV rv2);
+	//! Add (concat) another variable to the current one
+	RV add(RV rv2);
+	//! Subtract  another variable from the current one
+	RV subt(RV rv2);
+	//! Select only variables at indeces ind
+	RV subselect(ivec ind);
+	//! Select only variables at indeces ind
+	RV operator()(ivec ind);
+	//! Generate new \c RV with \c time shifted by delta.
+	void t(int delta);
+};
+
+
+
+
+//! Class representing function of variables
+class fnc {
+	RV rv;
+};
+
+//! Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
+class BM {
+public:
+	//!Logarithm of marginalized data likelihood.
+	double ll;
+
+	/*! \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.
+class epdf {
+	RV rv;
+public:
+	//! Returns the required moment of the epdf
+//	virtual vec moment ( const int order = 1 );
+	//! Returns a sample from the density, $x \sim epdf(rv)$
+	virtual vec sample (){};
+	virtual double eval(const vec &val){};
+};
+
+//! Conditional probability density, e.g. modeling some dependencies.
+class mpdf {
+	//! modeled random variable
+	RV rv;
+	//! random variable in condition
+	RV rvc;
+public:
+
+	//! Returns the required moment of the epdf
+//	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 void condition (vec &cond){};
+};
+
+
+
+#endif // BM_H
Index: /bdm/stat/libFN.cpp
===================================================================
--- /bdm/stat/libFN.cpp (revision 14)
+++ /bdm/stat/libFN.cpp (revision 14)
@@ -0,0 +1,14 @@
+#include <itpp/itbase.h>
+#include "libFN.h"
+
+using namespace itpp;
+
+using std::endl;
+
+constfn::constfn(RV rv){
+
+}
+
+vec constfn::evalsome(ivec rvind){
+	
+}
Index: /tests/test0.cpp
===================================================================
--- /tests/test0.cpp (revision 14)
+++ /tests/test0.cpp (revision 14)
@@ -0,0 +1,47 @@
+#include <itpp/itbase.h>
+#include "libBM.h"
+#include "libDC.h"
+
+using namespace itpp;
+
+//These lines are needed for use of cout and endl
+using std::cout;
+using std::endl;
+
+int main()
+{
+
+	RV th = RV ( "1 2","{a b }","1 1","0 0","0 0" );
+	RV r = RV ( "3 4" );
+	int i =0;
+
+	cout << th << r << endl;
+
+	ldmat ld = ldmat("1 0;-0.1 1","1.1 1.3");
+	vec v = "1 -0.1";
+	
+	cout << "Test to_mat"<<endl;
+//	cout << ld << endl;
+	cout << ld.to_mat()<<endl;
+
+	cout << "Test opupdt" <<endl;
+	ldmat ldup = ld;
+	ldup.opupdt(v,1);
+	cout << ldup.to_mat()<<endl;
+
+	cout << "Test +="<<endl;
+	ldmat ld2x = ld;
+	ld2x+=ld;
+	cout << ld.to_mat() << ld2x.to_mat() <<endl;
+
+	cout << "Test ldinv()"<<endl;
+	ldmat Il = ld; 
+	ld.inv(Il); // 
+	mat I = Il.to_mat()*ld.to_mat();
+		cout << "ld:"<<Il.to_mat() << "eye:"<< I <<endl;
+
+
+	//Exit program:
+	return 0;
+
+}
Index: /tests/testSmp.cpp
===================================================================
--- /tests/testSmp.cpp (revision 14)
+++ /tests/testSmp.cpp (revision 14)
@@ -0,0 +1,34 @@
+#include <itpp/itbase.h>
+#include "libEF.h"
+
+using namespace itpp;
+
+//These lines are needed for use of cout and endl
+using std::cout;
+using std::endl;
+
+int main() {
+
+	//RNG_randomize();
+
+	RV rv("1","{x }","2","0","0");
+	int N = 10000; //number of samples
+	vec mu0 = "1.5 1.7";
+	ldmat R = ldmat("1.2 0.3; 0.3 5");
+	enorm<ldmat> eN(rv,mu0,R);
+	
+	mat Smp = eN.sample(N);
+	cout << "True:" <<endl;
+	cout << "mu:" << mu0 ;
+	cout << "R:" << R ;
+	cout << "R:" << R.to_mat() ;
+	
+	vec Emu = Smp*ones(N) /N ;
+	cout << "Empirical:" <<endl;
+	cout << "mu:" << Emu;
+	cout << "R:" << (Smp*Smp.transpose())/N - Emu*Emu.transpose() <<endl;
+
+	//Exit program:
+	return 0;
+
+}
Index: /tests/testKF.cpp
===================================================================
--- /tests/testKF.cpp (revision 14)
+++ /tests/testKF.cpp (revision 14)
@@ -0,0 +1,63 @@
+#include <itpp/itbase.h>
+#include "libKF.h"
+#include "libDC.h"
+
+using namespace itpp;
+
+//These lines are needed for use of cout and endl
+using std::cout;
+using std::endl;
+
+int main() {
+
+	// Klaman filter
+	mat A, B,C,D,R,Q,P0;
+	vec mu0;
+	mat Mu0;;
+	// input from Matlab
+	it_file fin( "testKF.it" );
+
+	mat Dt, Xt,Xt2;
+	int Ndat;
+
+	fin >> Name( "d" ) >> Dt;
+	fin.seek( "A" ); 
+	fin >> A;
+	fin.seek( "B" ); 
+	fin >> B;
+	fin.seek( "C" ); 
+	fin >> C;
+	fin.seek( "D" ); 
+	fin >> D;
+	fin.seek( "R" ); 
+	fin >> R;
+	fin.seek( "Q" ); fin >> Q;
+	fin.seek( "P0" ); fin >> P0;
+	fin.seek( "mu0" ); fin >> Mu0; mu0=Mu0;
+	
+	Ndat = Dt.cols();
+	Xt=zeros( 2,Ndat );
+	Xt2=zeros( 2,Ndat );
+
+	Kalman<ldmat> KF( A,B,C,D,ldmat(R),ldmat(Q),ldmat(P0),mu0 );
+//	cout << KF;
+	KalmanFull KF2( A,B,C,D,R,Q,P0,mu0 );
+
+	Xt.set_col( 0,KF.mu );
+	Xt2.set_col( 0,KF.mu );
+	for ( int t=1;t<Ndat;t++ ) {
+		KF.bayes( Dt.get_col( t ));
+		KF2.bayes( Dt.get_col( t ));
+		Xt.set_col(t,KF.mu);
+		Xt2.set_col(t,KF2.mu);
+//		Kmu = KF.mu;
+//		cout <<  "t:" <<t<< "  " << dt<<"  "<<Kmu <<endl;
+	}
+
+	it_file fou( "testKF_res.it" );
+	fou << Name("xth") << Xt;
+	fou << Name("xth2") << Xt2;
+	//Exit program:
+	return 0;
+
+}
Index: /tests/testPF.cpp
===================================================================
--- /tests/testPF.cpp (revision 14)
+++ /tests/testPF.cpp (revision 14)
@@ -0,0 +1,41 @@
+#include <itpp/itbase.h>
+#include "libPF.h"
+#include "libDC.h"
+#include "libEF.h"
+
+using namespace itpp;
+
+//These lines are needed for use of cout and endl
+using std::cout;
+using std::endl;
+
+int main() {
+
+	RV x("1");
+	RV xm=x; xm.t(-1);
+	RV y("2");
+	
+	mat A = "1";
+	ldmat R("1");
+	
+	
+	vec ptc = randn(10);
+	vec w = exp(-0.5*(pow(ptc,2))/0.2);
+	
+	cout << "p:" << ptc << endl;
+	cout << "w:" << w << endl;
+	
+	PF pf(w);
+	ivec ind = pf.resample();
+	
+	cout << ind << endl;
+	/*
+	mlnorm<ldmat> obs(x,xm,A,R);
+	mlnorm<ldmat> par(y,x,A,R);
+	
+	TrivialPF TPF(obs,par,10);
+	*/
+	//Exit program:
+	return 0;
+
+}
Index: /mixpp.kdevelop
===================================================================
--- /mixpp.kdevelop (revision 11)
+++ /mixpp.kdevelop (revision 14)
@@ -18,5 +18,5 @@
     <run>
       <directoryradio>executable</directoryradio>
-      <mainprogram>/home/smidl/work/mixpp/testPF</mainprogram>
+      <mainprogram>/home/smidl/work/mixpp/testSmp</mainprogram>
       <programargs></programargs>
       <globaldebugarguments></globaldebugarguments>
@@ -221,11 +221,10 @@
     </groups>
     <tree>
-      <showvcsfields>false</showvcsfields>
-      <hidenonprojectfiles>true</hidenonprojectfiles>
       <hidepatterns>*.o,*.lo,CVS</hidepatterns>
+      <hidenonprojectfiles>false</hidenonprojectfiles>
     </tree>
   </kdevfileview>
   <ctagspart>
-    <customArguments/>
+    <customArguments></customArguments>
     <customTagfilePath>/home/smidl/work/mixpp/tags</customTagfilePath>
     <activeTagsFiles/>
Index: /stSmp.cpp
===================================================================
--- /testSmp.cpp (revision 13)
+++  (revision )
@@ -1,34 +1,0 @@
-#include <itpp/itbase.h>
-#include "libEF.h"
-
-using namespace itpp;
-
-//These lines are needed for use of cout and endl
-using std::cout;
-using std::endl;
-
-int main() {
-
-	//RNG_randomize();
-
-	RV rv("1","{x }","2","0","0");
-	int N = 10000; //number of samples
-	vec mu0 = "1.5 1.7";
-	ldmat R = ldmat("1.2 0.3; 0.3 5");
-	enorm<ldmat> eN(rv,mu0,R);
-	
-	mat Smp = eN.sample(N);
-	cout << "True:" <<endl;
-	cout << "mu:" << mu0 ;
-	cout << "R:" << R ;
-	cout << "R:" << R.to_mat() ;
-	
-	vec Emu = Smp*ones(N) /N ;
-	cout << "Empirical:" <<endl;
-	cout << "mu:" << Emu;
-	cout << "R:" << (Smp*Smp.transpose())/N - Emu*Emu.transpose() <<endl;
-
-	//Exit program:
-	return 0;
-
-}
Index: /st0.cpp
===================================================================
--- /test0.cpp (revision 7)
+++  (revision )
@@ -1,47 +1,0 @@
-#include <itpp/itbase.h>
-#include "libBM.h"
-#include "libDC.h"
-
-using namespace itpp;
-
-//These lines are needed for use of cout and endl
-using std::cout;
-using std::endl;
-
-int main()
-{
-
-	RV th = RV ( "1 2","{a b }","1 1","0 0","0 0" );
-	RV r = RV ( "3 4" );
-	int i =0;
-
-	cout << th << r << endl;
-
-	ldmat ld = ldmat("1 0;-0.1 1","1.1 1.3");
-	vec v = "1 -0.1";
-	
-	cout << "Test to_mat"<<endl;
-//	cout << ld << endl;
-	cout << ld.to_mat()<<endl;
-
-	cout << "Test opupdt" <<endl;
-	ldmat ldup = ld;
-	ldup.opupdt(v,1);
-	cout << ldup.to_mat()<<endl;
-
-	cout << "Test +="<<endl;
-	ldmat ld2x = ld;
-	ld2x+=ld;
-	cout << ld.to_mat() << ld2x.to_mat() <<endl;
-
-	cout << "Test ldinv()"<<endl;
-	ldmat Il = ld; 
-	ld.inv(Il); // 
-	mat I = Il.to_mat()*ld.to_mat();
-		cout << "ld:"<<Il.to_mat() << "eye:"<< I <<endl;
-
-
-	//Exit program:
-	return 0;
-
-}
Index: /bFN.h
===================================================================
--- /libFN.h (revision 11)
+++  (revision )
@@ -1,27 +1,0 @@
-//
-// C++ Interface: itpp_ext
-//
-// Description: 
-//
-//
-// Author: smidl <smidl@utia.cas.cz>, (C) 2008
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-#include <itpp/itbase.h>
-#include "libBM"
-
-using namespace itpp;
- 
-class constfn : public fnc {
-	RV rv;
-	ivec indexlist; // needed by evalsome
-	vec val;
-
-public:
-	vec eval() {return val;};
-	vec eval(vec &cond) {return val;};
-	
-	vec evalsome(ivec &rvind);
-}
Index: /stKF.cpp
===================================================================
--- /testKF.cpp (revision 11)
+++  (revision )
@@ -1,63 +1,0 @@
-#include <itpp/itbase.h>
-#include "libKF.h"
-#include "libDC.h"
-
-using namespace itpp;
-
-//These lines are needed for use of cout and endl
-using std::cout;
-using std::endl;
-
-int main() {
-
-	// Klaman filter
-	mat A, B,C,D,R,Q,P0;
-	vec mu0;
-	mat Mu0;;
-	// input from Matlab
-	it_file fin( "testKF.it" );
-
-	mat Dt, Xt,Xt2;
-	int Ndat;
-
-	fin >> Name( "d" ) >> Dt;
-	fin.seek( "A" ); 
-	fin >> A;
-	fin.seek( "B" ); 
-	fin >> B;
-	fin.seek( "C" ); 
-	fin >> C;
-	fin.seek( "D" ); 
-	fin >> D;
-	fin.seek( "R" ); 
-	fin >> R;
-	fin.seek( "Q" ); fin >> Q;
-	fin.seek( "P0" ); fin >> P0;
-	fin.seek( "mu0" ); fin >> Mu0; mu0=Mu0;
-	
-	Ndat = Dt.cols();
-	Xt=zeros( 2,Ndat );
-	Xt2=zeros( 2,Ndat );
-
-	Kalman<ldmat> KF( A,B,C,D,ldmat(R),ldmat(Q),ldmat(P0),mu0 );
-//	cout << KF;
-	KalmanFull KF2( A,B,C,D,R,Q,P0,mu0 );
-
-	Xt.set_col( 0,KF.mu );
-	Xt2.set_col( 0,KF.mu );
-	for ( int t=1;t<Ndat;t++ ) {
-		KF.bayes( Dt.get_col( t ));
-		KF2.bayes( Dt.get_col( t ));
-		Xt.set_col(t,KF.mu);
-		Xt2.set_col(t,KF2.mu);
-//		Kmu = KF.mu;
-//		cout <<  "t:" <<t<< "  " << dt<<"  "<<Kmu <<endl;
-	}
-
-	it_file fou( "testKF_res.it" );
-	fou << Name("xth") << Xt;
-	fou << Name("xth2") << Xt2;
-	//Exit program:
-	return 0;
-
-}
Index: /bPF.h
===================================================================
--- /libPF.h (revision 11)
+++  (revision )
@@ -1,72 +1,0 @@
-/*!
-  \file
-  \brief Bayesian Filtering using stochastic sampling (Particle Filters)
-  \author Vaclav Smidl.
-
-  -----------------------------------
-  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
-
-  Using IT++ for numerical operations
-  -----------------------------------
-*/
-
-#ifndef PF_H
-#define PF_H
-
-#include <itpp/itbase.h>
-#include "libBM.h"
-#include "libDC.h"
-
-using namespace itpp;
-
-enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
-
-/*!
-* \brief A Particle Filter prototype
-
-Bayesian Filtering equations hold.
-*/
-class PF : public BM { 
-protected:
-	int n; //number of particles
-	vec w; //particle weights
-	Uniform_RNG URNG; //used for resampling
-	
-public:
-	//! Returns indexes of particles that should be resampled. The ordering MUST guarantee inplace replacement. (Important for MPF.)
-	ivec resample(RESAMPLING_METHOD method = SYSTEMATIC);
-	PF (vec w);
-	//TODO remove or implement bayes()!
-	void bayes(const vec &dt, bool evell){};
-};
-
-/*!
-* \brief Trivial particle filter with proposal density that is not conditioned on the data.
-
-
-*/
-
-class TrivialPF : public PF {
-	Array<vec> ptcls;
-	
-	bool is_proposal;
-	mpdf *prop;
-	mpdf *par;
-	mpdf *obs;
-	
-	public:
-	TrivialPF(mpdf &par, mpdf &obs, mpdf &prop, int n0);
-	TrivialPF(mpdf &par, mpdf &obs, int n0);
-	void bayes(const vec &dt, bool evalll);
-};
-
-class MPF : public TrivialPF {
-	Array<BM> Bms;
-	public:
-	MPF(BM &B, mpdf &prop, mpdf &obs, mpdf &par);
-	void bayes(vec &dt);	
-};
-
-#endif // KF_H
-
-
Index: /stPF.cpp
===================================================================
--- /testPF.cpp (revision 11)
+++  (revision )
@@ -1,41 +1,0 @@
-#include <itpp/itbase.h>
-#include "libPF.h"
-#include "libDC.h"
-#include "libEF.h"
-
-using namespace itpp;
-
-//These lines are needed for use of cout and endl
-using std::cout;
-using std::endl;
-
-int main() {
-
-	RV x("1");
-	RV xm=x; xm.t(-1);
-	RV y("2");
-	
-	mat A = "1";
-	ldmat R("1");
-	
-	
-	vec ptc = randn(10);
-	vec w = exp(-0.5*(pow(ptc,2))/0.2);
-	
-	cout << "p:" << ptc << endl;
-	cout << "w:" << w << endl;
-	
-	PF pf(w);
-	ivec ind = pf.resample();
-	
-	cout << ind << endl;
-	/*
-	mlnorm<ldmat> obs(x,xm,A,R);
-	mlnorm<ldmat> par(y,x,A,R);
-	
-	TrivialPF TPF(obs,par,10);
-	*/
-	//Exit program:
-	return 0;
-
-}
Index: /bDC.cpp
===================================================================
--- /libDC.cpp (revision 12)
+++  (revision )
@@ -1,361 +1,0 @@
-#include <itpp/itbase.h>
-#include "libDC.h"
-
-using namespace itpp;
-
-using std::endl;
-
-//! Auxiliary function dydr; dyadic reduction
-void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx );
-
-//! Auxiliary function ltuinv; inversion of a triangular matrix;
-//TODO can be done via: dtrtri.f from lapack
-mat ltuinv( const mat &L );
-
-
-
-ldmat::ldmat( const mat &exL, const vec &exD ) {
-	D = exD;
-	L = exL;
-	dim = exD.length();
-}
-
-ldmat::ldmat() {
-	vec D ;
-	mat L;
-	dim = 0;
-}
-
-ldmat::ldmat( const mat V ) {
-//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!" );
-
-	
-}
-
-void ldmat::opupdt( const vec &v,  double w ) {
-	int dim = D.length();
-	double kr;
-	vec r = v;
-	//beware! it is potentionally dangerous, if ITpp change _behaviour of _data()!
-	double *Lraw = L._data();
-	double *Draw = D._data();
-	double *rraw = r._data();
-
-	it_assert_debug( v.length() == dim, "LD::ldupdt vector v is not compatible with this ld." );
-
-	for ( int i = dim - 1; i >= 0; i-- ) {
-		dydr( rraw, Lraw + i, &w, Draw + i, rraw + i, 0, i, &kr, 1, dim );
-	}
-}
-
-std::ostream &operator<< ( std::ostream &os,  ldmat &ld ) {
-	os << "L:" << ld.L << endl;
-	os << "D:" << ld.D << endl;
-}
-
-mat ldmat::to_mat() {
-	int dim = D.length();
-	mat V( dim, dim );
-	double sum;
-	int r, c, cc;
-
-	for ( r = 0;r < dim;r++ ) { //row cycle
-		for ( c = r;c < dim;c++ ) {
-			//column cycle, using symmetricity => c=r!
-			sum = 0.0;
-			for ( cc = c;cc < dim;cc++ ) { //cycle over the remaining part of the vector
-				sum += L( cc, r ) * D( cc ) * L( cc, c );
-				//here L(cc,r) = L(r,cc)';
-			}
-			V( r, c ) = sum;
-			// symmetricity
-			if ( r != c ) {V( c, r ) = sum;};
-		}
-	}
-	mat V2 = L.transpose()*diag( D )*L;
-	return V2;
-}
-
-
-void ldmat::add( const ldmat &ld2, double w ) {
-	int dim = D.length();
-
-	it_assert_debug( ld2.D.length() == dim, "LD.add() incompatible sizes of LDs;" );
-
-	//Fixme can be done more efficiently either via dydr or ldform
-	for ( int r = 0; r < dim; r++ ) {
-		// Add columns of ld2.L' (i.e. rows of ld2.L) as dyads weighted by ld2.D
-		this->opupdt( ld2.L.get_row( r ), w*ld2.D( r ) );
-	}
-}
-
-void ldmat::clear(){L.clear(); for ( int i=0;i<L.cols();i++ ){L( i,i )=1;}; D.clear();}
-
-void ldmat::inv( ldmat &Inv ) {
-	int dim = D.length();
-	Inv.clear();   //Inv = zero in LD
-	mat U = ltuinv( L );
-
-	//Fixme can be done more efficiently either via dydr or ldform
-	for ( int r = 0; r < dim; r++ ) {
-		// Add columns of U as dyads weighted by 1/D
-		Inv.opupdt( U.get_col( r ), 1.0 / D( r ) );
-	}
-}
-
-void ldmat::mult_sym( const mat &C, bool trans ) {
-
-//TODO better
-
-	it_assert_debug( C.cols()==L.cols(), "ldmat::mult_sym wrong input argument" );
-	mat Ct=C;
-
-	if ( trans==false ) { // return C*this*C'
-		Ct *= this->to_mat();
-		Ct *= C.transpose();
-	} else {	// return C'*this*C
-		Ct = C.transpose();
-		Ct *= this->to_mat();
-		Ct *= C;
-	}
-
-	ldmat Lnew=ldmat( Ct );
-	L = Lnew.L;
-	D = Lnew.D;
-}
-
-void ldmat::mult_sym( const mat &C, ldmat &U, bool trans ) {
-
-//TODO better
-
-//TODO input test
-
-	mat Ct=C;
-
-	if ( trans==false ) { // return C*this*C'
-		Ct *= U.to_mat();
-		Ct *= C.transpose();
-	} else {	// return C'*this*C
-		Ct = C.transpose();
-		Ct *= U.to_mat();
-		Ct *= C;
-	}
-
-	ldmat Lnew=ldmat( Ct );
-	L = Lnew.L;
-	D = Lnew.D;
-}
-
-double ldmat::logdet() {
-	double ldet = 0.0;
-	int i;
-// sum logarithms of diagobal elements
-	for ( i=0; i<D.length(); i++ ){ldet+=log( D( i ) );};
-}
-
-double ldmat::qform( vec &v ) {
-	double x = 0.0, sum;
-	int i,j;
-
-	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 );}
-		x +=D( i )*sum*sum;
-	};
-	return x;
-}
-
-ldmat& ldmat::operator *= ( double x ) {
-	int i;
-	for ( i=0;i<D.length();i++ ){D( i )*=x;};
-}
-
-vec ldmat::sqrt_mult( vec &x ) {
-	int i,j;
-	vec res( dim );
-	double sum;
-	for ( i=0;i<dim;i++ ) {//for each element of result
-		res( i ) = 0.0;
-		for ( j=i;j<dim;j++ ) {//sum D(j)*L(:,i).*x
-			res( i ) += sqrt( D( j ) )*L( j,i )*x( j );
-		}
-	}
-	vec res2 = L.transpose()*diag( sqrt( D ) )*x;
-	return res2;
-}
-
-void ldmat::ldform( mat &A,vec &D0 ) {
-	int m = A.rows();
-	int n = A.cols();
-	int mn = (m<n) ? m :n ;
-
-	it_assert_debug( A.cols()==dim,"ldmat::ldform A is not compatible" );
-	it_assert_debug( D.length()==A.rows(),"ldmat::ldform Vector D must have the length as row count of A" );
-
-	L=concat_vertical( zeros( n,n ), diag( sqrt( D0 ) )*A );
-	D=zeros( n+m );
-	
-	//unnecessary big L and D will be made smaller at the end of file
-	
-	vec w=zeros( n );
-	vec v;
-	double sum, beta, pom;
-
-	int cc=0;
-	int i=n; // +1 in .m
-	int ii,j,jj;
-	while (( i> ( n-mn+1-cc ) )&&( i>1 ) ) {
-		i--;
-		sum = 0.0;
-		v.set_size( m+i-( n-cc ) ); //prepare v
-		for ( ii=n-cc;ii<m+i;i++ ) {
-			sum+= L( ii,i )*L( ii,i );
-			v( ii )=L( ii,i ); //assign v
-		}
-		if ( L( m+i,i )==0 ) {
-			beta = sqrt( sum );
-		} else {
-			beta = L( m+i,i )+sign( L( m+i,i ) )*sqrt( sum ) ;
-		}
-		if ( std::fabs( beta )<eps ) {
-			cc++;
-			L.set_row( n+1-cc, L.get_row( m+i ) );
-			L.set_row( m+i,0 );
-			D( m+i )=0; L( m+i,i )=1;
-			L.set_submatrix( n+1-cc,m+i,i,i,0 );
-			continue;
-		}
-
-		sum-=v( v.length()-1 )*v( v.length()-1 ); //
-		sum/=beta*beta;
-		sum++;
-
-		v/=beta;
-		v( v.length()-1 )=1;
-
-		pom=-2/sum;
-		for ( j=i;i>=0;i-- ) {
-			w( j )=0.0;
-			for ( ii=n-cc;ii<m+i;ii++ ) {
-				w( j )+= v( ii )*L( ii,j );
-			}
-			w( j )*=pom;
-		}
-
-		for ( ii=n-cc;ii<m+i;ii++ ) {
-			for ( jj=0;jj<i-1;jj++ ) {
-				L( ii,jj )+= v( ii )*w( jj );
-			}
-		}
-		for ( ii=n-cc;ii<m+i;ii++ ) {
-			L( ii,i )= 0;
-		}
-		L( m+i,i )+=w( i );
-
-		D( m+i )=L( m+i,i )*L( m+i,i );
-		for ( ii=0;ii<i;ii++ ) {
-			L( m+i,ii )/=L( m+i,i );
-		}
-	}
-	if ( i>0 ) {
-		for ( ii=0;ii<i-1;ii++ ) {
-			jj = D.length()-1-n+ii;
-			L.set_row(jj,0);
-			L(jj,jj)=1;
-		}
-	}
-
-	//cut-out L and D;
-	L.del_rows(0,m-1);
-	D.del(0,m-1);
-}
-
-//////// Auxiliary Functions
-
-mat ltuinv( const mat &L ) {
-	int dim = L.cols();
-	mat Il = eye( dim );
-	int i, j, k, m;
-	double s;
-
-//Fixme blind transcription of ltuinv.m
-	for ( k = 1; k < ( dim );k++ ) {
-		for ( i = 0; i < ( dim - k );i++ ) {
-			j = i + k; //change in .m 1+1=2, here 0+0+1=1
-			s = L( j, i );
-			for ( m = i + 1; m < ( j - 1 ); m++ ) {
-				s += L( m, i ) * Il( j, m );
-			}
-			Il( j, i ) = -s;
-		}
-	}
-
-	return Il;
-}
-
-void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx )
-/********************************************************************
-
-   dydr = dyadic reduction, performs transformation of sum of
-          2 dyads r*Dr*r'+ f*Df*f' so that the element of r pointed
-          by R is zeroed. This version allows Dr to be NEGATIVE. Hence the name negdydr or dydr_withneg.
-
-   Parameters :
-     r ... pointer to reduced dyad
-     f ... pointer to reducing dyad
-     Dr .. pointer to the weight of reduced dyad
-     Df .. pointer to the weight of reducing dyad
-     R ... pointer to the element of r, which is to be reduced to
-           zero; the corresponding element of f is assumed to be 1.
-     jl .. lower index of the range within which the dyads are
-           modified
-     ju .. upper index of the range within which the dyads are
-           modified
-     kr .. pointer to the coefficient used in the transformation of r
-           rnew = r + kr*f
-     m  .. number of rows of modified matrix (part of which is r)
-  Remark : Constant mzero means machine zero and should be modified
-           according to the precision of particular machine
-
-                                                 V. Peterka 17-7-89
-
-  Added:
-     mx .. number of rows of modified matrix (part of which is f)  -PN
-
-********************************************************************/
-{
-	int j, jm;
-	double kD, r0;
-	double mzero = 2.2e-16;
-	double threshold = 1e-4;
-
-	if ( fabs( *Dr ) < mzero ) *Dr = 0;
-	r0 = *R;
-	*R = 0.0;
-	kD = *Df;
-	*kr = r0 * *Dr;
-	*Df = kD + r0 * ( *kr );
-	if ( *Df > mzero ) {
-		kD /= *Df;
-		*kr /= *Df;
-	} else {
-		kD = 1.0;
-		*kr = 0.0;
-		if ( *Df < -threshold ) {
-			it_warning( "Problem in dydr: subraction of dyad results in negative definitness. Likely mistake in calling function." );
-		}
-		*Df = 0.0;
-	}
-	*Dr *= kD;
-	jm = mx * jl;
-	for ( j = m * jl; j < m*jh; j += m ) {
-		r[j] -=  r0 * f[jm];
-		f[jm] += *kr * r[j];
-		jm += mx;
-	}
-}
-
-
Index: /bEF.cpp
===================================================================
--- /libEF.cpp (revision 13)
+++  (revision )
@@ -1,7 +1,0 @@
-#include <itpp/itbase.h>
-#include "libEF.h"
-
-using namespace itpp;
-
-using std::cout;
-
Index: /bDC.h
===================================================================
--- /libDC.h (revision 12)
+++  (revision )
@@ -1,175 +1,0 @@
-/*!
- * \file
- * \brief Matrices in decomposed forms (LDL', LU, UDU', etc).
- * \author Vaclav Smidl.
- *
- * -----------------------------------
- * BDM++ - C++ library for Bayesian Decision Making under Uncertainty
- *
- * Using IT++ for numerical operations
- * -----------------------------------
- */
-
-#ifndef DC_H
-#define DC_H
-
-#include <itpp/itbase.h>
-
-using namespace itpp;
-
-/*! \brief Virtual class for representation of double symmetric matrices in square-root form.
-
-All operations defined on this class should be optimized for the chosed decomposition.
-*/
-class sqmat {
-public:
-	/*!
-	 * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$.
-	 * @param  v Vector forming the outer product to be added
-	 * @param  w weight of updating; can be negative
-
-	 BLAS-2b operation.
-	 */
-	virtual void opupdt( const vec &v, double w ) =0;
-
-	/*! \brief Conversion to full matrix.
-	*/
-
-	virtual mat to_mat() =0;
-
-	/*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$
-	@param C multiplying matrix,
-	@param trans if true, product $V = C'*V*C$ will be computed instead;
-	*/
-	virtual void mult_sym( const mat &C, bool trans=true ) =0;
-
-	
-	/*!
-	\brief Logarithm of a determinant.
-	
-	*/
-	virtual double logdet() =0;
-
-	/*!
-	\brief Multiplies square root of $V$ by vector $x$.
-	
-	Used e.g. in generating normal samples.
-	*/
-	virtual vec sqrt_mult(vec &v) =0;
-	
-	/*!
-	\brief Evaluates quadratic form $x= v'*V*v$;
-	
-	*/
-	virtual double qform(vec &v) =0;
-
-//	//! easy version of the 
-//	sqmat inv();
-
-	//! Clearing matrix so that it corresponds to zeros.
-	virtual void clear() =0;
-	
-	//! Reimplementing common functions of mat: cols().
-	virtual int cols() =0;
-
-	//! Reimplementing common functions of mat: cols().
-	virtual int rows() =0;
-
-protected:
-	int dim;
-};
-
-
-/*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix.
-
-This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices;
-*/
-class fsqmat: sqmat {
-	void opupdt( const vec &v, double w );
-	mat to_mat();
-	void mult_sym( const mat &C, bool trans=false );
-   void mult_sym( const mat &C, fsqmat &U, bool trans=false );
-	void inv(fsqmat &Inv);
-	void clear();
-	
-	
-	//! Constructor
-	fsqmat(const mat &M);
-	
-	/*! \brief Matrix inversion preserving the chosen form.
-
-	@param Inv a space where the inverse is stored.
-	
-	*/
-	virtual void inv(fsqmat* Inv);
-	
-
-};
-
-class ldmat: sqmat {
-public:
-
-	//! Construct by copy of L and D.
-	ldmat( const mat &L, const vec &D );
-	//! Construct by decomposition of full matrix V.
-	ldmat( mat V );
-	ldmat ();
-
-	// Reimplementation of compulsory operatios
-
-	void opupdt( const vec &v, double w );
-	mat to_mat();
-   void mult_sym( const mat &C, bool trans=false );
-	void add ( const ldmat &ld2, double w=1.0 );
-	double logdet();
-	double qform(vec &v);
-//	sqmat& operator -= ( const sqmat & ld2 );
-	void clear();
-	int cols();
-	int rows();
-	vec sqrt_mult(vec &v);
-
-	/*! \brief Matrix inversion preserving the chosen form.
-
-	@param Inv a space where the inverse is stored.
-	
-	*/
-	virtual void inv(ldmat &Inv);
-	
-	/*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class.
-
-	@param Inv a space where the inverse is stored.
-	
-	*/
-   void mult_sym( const mat &C, ldmat &U, bool trans=false );
-
-	/*! \brief Transforms general $A'D0A$ into pure $L'DL$
-	
-	The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$
-
-	@param A general matrix
-	@param D0 general vector
-	
-	*/
-   void ldform( mat &A, vec &D0 );
-
-	ldmat& operator += (const ldmat &ldA);
-	ldmat& operator -= (const ldmat &ldA);
-	ldmat& operator *= (double x);
-	
-	friend std::ostream &operator<< ( std::ostream &os, ldmat &sq );
-
-protected:
-	vec D;
-	mat L;
-
-};
-
-//////// Operations:
-
-inline ldmat& ldmat::operator += (const ldmat &ldA)  {this->add(ldA);return *this;}
-inline ldmat& ldmat::operator -= (const ldmat &ldA)  {this->add(ldA,-1.0);return *this;}
-inline int ldmat::cols(){return L.cols();}
-inline int ldmat::rows(){return L.rows();}
-
-#endif // DC_H
Index: /bBM.cpp
===================================================================
--- /libBM.cpp (revision 12)
+++  (revision )
@@ -1,69 +1,0 @@
-#include <itpp/itbase.h>
-#include "libBM.h"
-#include "itpp_ext.h"
-
-using namespace itpp;
-
-using std::cout;
-
-void RV::init( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs ) {
-	//
-	int i;
-	len = in_ids.length();
-	//PRUDENT_MODE
-	// All vetors should be of same length
-	if ( 	 ( len != in_ids.length() ) || \
-	        ( len != in_names.length() ) || \
-	        ( len != in_sizes.length() ) || \
-	        ( len != in_times.length() ) || \
-	        ( len != in_obs.length() ) ) {
-		it_error( "RV::RV inconsistent length of input vectors." );
-	}
-
-	ids = in_ids;
-	names = in_names;
-	sizes = in_sizes;
-	times = in_times;
-	obs = in_obs;
-	size = 0;
-	for(i=0;i<len;i++){size+=sizes(i);}
-};
-
-RV::RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs ) {
-	init ( in_ids, in_names, in_sizes, in_times, in_obs );
-}
-
-RV::RV () {};
-
-RV::RV ( ivec in_ids ) {
-	
-	len = in_ids.length();
-	Array<std::string> A( len );
-	std::string rvstr = "rv";
-
-	for ( int i = 0; i < len;i++ ) {
-		A( i ) = rvstr + to_str(i);
-	}
-
-	init ( in_ids, A, ones_i( len ), zeros_i( len ), zeros_i( len ) );
-}
-
-RV RV::subselect(ivec ind){
-	return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind), obs(ind));
-}
-
-void RV::t(int delta){ times +=delta;}
-
-RV RV::operator()(ivec ind){
-	return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind), obs(ind));
-}
-
-std::ostream &operator<<( std::ostream &os, const RV &rv ) {
-
-	for ( int i = 0; i < rv.len ;i++ ) {
-		os << rv.ids( i ) << "(" << rv.sizes( i ) << ")" <<  // id(size)=
-		"=" << rv.names( i )  << "_{"  << rv.times( i ) << "}; "; //name_{time}
-	}
-	return os;
-}
-
Index: /bEF.h
===================================================================
--- /libEF.h (revision 12)
+++  (revision )
@@ -1,168 +1,0 @@
-/*!
-  \file
-  \brief Probability distributions for Exponential Family models.
-  \author Vaclav Smidl.
-
-  -----------------------------------
-  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
-
-  Using IT++ for numerical operations
-  -----------------------------------
-*/
-
-#ifndef EF_H
-#define EF_H
-
-#include <itpp/itbase.h>
-#include "libDC.h"
-#include "libBM.h"
-//#include <std>
-
-using namespace itpp;
-
-/*!
-* \brief General conjugate exponential family posterior density.
-
-* More?...
-*/
-class eEF : public epdf {
-
-public:
-	virtual void tupdate( double phi, mat &vbar, double nubar ) {};
-	virtual void dupdate( mat &v,double nu=1.0 ) {};
-};
-
-class mEF : public mpdf {
-
-public:
-
-};
-
-/*!
-* \brief General exponential family density
-
-* More?...
-*/
-template<class sq_T>
-class enorm : public eEF {
-	int dim;
-	vec mu;
-	sq_T R;
-public:
-	enorm( RV &rv, vec &mu, sq_T &R );
-	enorm();
-	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;
-};
-
-/*!
- \brief
-*/
-template<class sq_T>
-class mlnorm : public mEF {
-	enorm<sq_T> epdf;
-	mat A;
-public:
-	//! Constructor
-	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 );
-};
-
-////////////////////////
-
-template<class sq_T>
-enorm<sq_T>::enorm( RV &rv, vec &mu0, sq_T &R0 ) {
-	dim = rv.count();
-	mu = mu0;
-	R = R0;
-};
-
-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 ) {
-	//
-};
-
-template<class sq_T>
-vec enorm<sq_T>::sample() {
-	vec x( dim );
-	RNG.sample_vector( dim,x );
-	vec smp = R.sqrt_mult( x );
-
-	smp += mu;
-	return smp;
-};
-
-template<class sq_T>
-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 );
-		pom +=mu;
-		X.set_col( i, pom);
-	}
-	return X;
-};
-
-template<class sq_T>
-double enorm<sq_T>::eval( const vec &val ) {
-	//
-};
-
-
-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 ) {
-	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 smp = epdf.sample();
-	lik = epdf.eval( smp );
-	return smp;
-}
-
-template<class sq_T>
-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 );
-	for ( i=0; i<dim; i++ ) {
-		smp = epdf.sample();
-		lik( i ) = epdf.eval( smp );
-		Smp.set_col( i ,smp );
-	}
-	return Smp;
-}
-
-template<class sq_T>
-void mlnorm<sq_T>::condition( vec &cond ) {
-	epdf.mu = A*cond;
-//R is already assigned;
-}
-
-#endif //EF_H
Index: /bKF.cpp
===================================================================
--- /libKF.cpp (revision 7)
+++  (revision )
@@ -1,51 +1,0 @@
-#include <itpp/itbase.h>
-#include "libBM.h"
-#include "libKF.h"
-
-using namespace itpp;
-
-using std::endl;
-
-KalmanFull::KalmanFull( mat A0, mat B0, mat C0, mat D0, mat R0, mat Q0, mat P0, vec mu0 ) {
-	dimx = A0.rows();
-	dimu = B0.cols();
-	dimy = C0.rows();
-
-	it_assert_debug( A0.cols()==dimx, "KalmanFull: A is not square" );
-	it_assert_debug( B0.rows()==dimx, "KalmanFull: B is not compatible" );
-	it_assert_debug( C0.cols()==dimx, "KalmanFull: C is not square" );
-	it_assert_debug(( D0.rows()==dimy ) || ( D0.cols()==dimu ),	"KalmanFull: D is not compatible" );
-	it_assert_debug(( R0.cols()==dimy ) || ( R0.rows()==dimy ), "KalmanFull: R is not compatible" );
-	it_assert_debug(( Q0.cols()==dimx ) || ( Q0.rows()==dimx ), "KalmanFull: Q is not compatible" );
-
-	A = A0;
-	B = B0;
-	C = C0;
-	D = D0;
-	R = R0;
-	Q = Q0;
-	mu = mu0;
-	P = P0;
-}
-
-void KalmanFull::bayes( const vec &dt, bool evalll ) {
-	it_assert_debug( dt.length()==( dimy+dimu ),"KalmanFull::bayes wrong size of dt" );
-
-	vec u = dt.get( dimy,dimy+dimu-1 );
-	vec y = dt.get( 0,dimy-1 );
-	//Time update
-	mu = A*mu + B*u;
-	P  = A*P*A.transpose() + Q;
-
-	//Data update
-	_Ry = C*P*C.transpose() + R;
-	_iRy = inv( _Ry );
-	_K = P*C.transpose()*_iRy;
-	P -= _K*C*P; // P = P -KCP;
-	mu += _K*( y-C*mu-D*u );
-};
-
-std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf ) {
-	os << "mu:" << kf.mu << endl << "P:" << kf.P <<endl;
-}
-
Index: /bFN.cpp
===================================================================
--- /libFN.cpp (revision 11)
+++  (revision )
@@ -1,14 +1,0 @@
-#include <itpp/itbase.h>
-#include "libFN.h"
-
-using namespace itpp;
-
-using std::endl;
-
-constfn::constfn(RV rv){
-
-}
-
-vec constfn::evalsome(ivec rvind){
-	
-}
Index: /bBM.h
===================================================================
--- /libBM.h (revision 12)
+++  (revision )
@@ -1,116 +1,0 @@
-/*!
-  \file
-  \brief Bayesian Models (bm) that use Bayes rule to learn from observations
-  \author Vaclav Smidl.
-
-  -----------------------------------
-  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
-
-  Using IT++ for numerical operations
-  -----------------------------------
-*/
-
-#ifndef BM_H
-#define BM_H
-
-#include <itpp/itbase.h>
-//#include <std>
-
-using namespace itpp;
-
-/*!
-* \brief Class representing variables, most often random variables
-
-* More?...
-*/
-class RV {
-	int size,len;
-	ivec ids;
-	ivec sizes;
-	ivec times;
-	ivec obs;
-	Array<std::string> names;
-
-private:
-	void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs );
-public:
-	//! Full constructor which is called by the others
-	RV ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs );
-	//! default constructor
-	RV ( ivec ids );
-	//! Empty constructor will be set later
-	RV ();
-	
-	//! Printing output e.g. for debugging.
-	friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
-
-	//! Return length (number of scalars) of the RV.
-	int count() {return size;}
-	//TODO why not inline and later??
-	
-	//! Find indexes of another rv in self
-	ivec find(RV rv2);
-	//! Add (concat) another variable to the current one
-	RV add(RV rv2);
-	//! Subtract  another variable from the current one
-	RV subt(RV rv2);
-	//! Select only variables at indeces ind
-	RV subselect(ivec ind);
-	//! Select only variables at indeces ind
-	RV operator()(ivec ind);
-	//! Generate new \c RV with \c time shifted by delta.
-	void t(int delta);
-};
-
-
-
-
-//! Class representing function of variables
-class fnc {
-	RV rv;
-};
-
-//! Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
-class BM {
-public:
-	//!Logarithm of marginalized data likelihood.
-	double ll;
-
-	/*! \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.
-class epdf {
-	RV rv;
-public:
-	//! Returns the required moment of the epdf
-//	virtual vec moment ( const int order = 1 );
-	//! Returns a sample from the density, $x \sim epdf(rv)$
-	virtual vec sample (){};
-	virtual double eval(const vec &val){};
-};
-
-//! Conditional probability density, e.g. modeling some dependencies.
-class mpdf {
-	//! modeled random variable
-	RV rv;
-	//! random variable in condition
-	RV rvc;
-public:
-
-	//! Returns the required moment of the epdf
-//	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 void condition (vec &cond){};
-};
-
-
-
-#endif // BM_H
Index: /x.cpp
===================================================================
--- /arx.cpp (revision 13)
+++  (revision )
@@ -1,6 +1,0 @@
-#include <itpp/itbase.h>
-#include "libBM.h"
-
-using namespace itpp;
-
-
Index: /bKF.h
===================================================================
--- /libKF.h (revision 8)
+++  (revision )
@@ -1,145 +1,0 @@
-/*!
-  \file
-  \brief Bayesian Filtering for linear Gaussian models (Kalman Filter) and extensions
-  \author Vaclav Smidl.
-
-  -----------------------------------
-  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
-
-  Using IT++ for numerical operations
-  -----------------------------------
-*/
-
-#ifndef KF_H
-#define KF_H
-
-#include <itpp/itbase.h>
-#include "libBM.h"
-#include "libDC.h"
-
-
-using namespace itpp;
-
-/*!
-* \brief Basic Kalman filter with full matrices (education purpose only)! Will be deleted soon!
-*/
-class KalmanFull : public BM { 
-	int dimx, dimy, dimu;
-	mat A, B, C, D, R, Q;
-	
-	//cache 
-	mat _Pp, _Ry, _iRy, _K;
-public:
-	//posterior 
-	//! Mean value of the posterior density
-	vec mu;
-	//! Variance of the posterior density
-	mat P;
-
-public:
-	//! Full constructor 
-	KalmanFull ( mat A, mat B, mat C, mat D, mat R, mat Q, mat P0, vec mu0);
-	//! Here dt = [yt;ut] of appropriate dimensions
-	void bayes(const vec &dt, bool evalll=true); 
-
-	friend std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf );
-
-};
-
-
-/*!
-* \brief Kalman filter with covaraince matrices in square root form.
-*/
-template<class sq_T>
-class Kalman : public BM { 
-	int dimx, dimy, dimu;
-	mat A, B, C, D;
-	sq_T R, Q;
-	
-	//cache
-	mat _K;
-	vec _yp;
-	sq_T _Ry,_iRy;
-public:
-	//posterior 
-	//! Mean value of the posterior density
-	vec mu;
-	//! Mean value of the posterior density
-	sq_T P;
-
-public:
-	//! Full constructor 
-	Kalman ( mat A0, mat B0, mat C0, mat D0, sq_T R0, sq_T Q0, sq_T P0, vec mu0 );
-	//! Here dt = [yt;ut] of appropriate dimensions
-	void bayes(const vec &dt, bool evalll=true); 
-
-	friend std::ostream &operator<< ( std::ostream &os, const KalmanFull &kf );
-
-};
-
-//////// INstance
-
-template<class sq_T>
-Kalman<sq_T>::Kalman( mat A0, mat B0, mat C0, mat D0, sq_T R0, sq_T Q0, sq_T P0, vec mu0 ) {
-	dimx = A0.rows();
-	dimu = B0.cols();
-	dimy = C0.rows();
-
-	it_assert_debug( A0.cols()==dimx, "Kalman: A is not square" );
-	it_assert_debug( B0.rows()==dimx, "Kalman: B is not compatible" );
-	it_assert_debug( C0.cols()==dimx, "Kalman: C is not square" );
-	it_assert_debug(( D0.rows()==dimy ) || ( D0.cols()==dimu ),	"Kalman: D is not compatible" );
-	it_assert_debug(( R0.cols()==dimy ) || ( R0.rows()==dimy ), "Kalman: R is not compatible" );
-	it_assert_debug(( Q0.cols()==dimx ) || ( Q0.rows()==dimx ), "Kalman: Q is not compatible" );
-
-	A = A0;
-	B = B0;
-	C = C0;
-	D = D0;
-	R = R0;
-	Q = Q0;
-	mu = mu0;
-	P = P0;
-
-	ll = 0;
-//Fixme should we assign cache??
-	_iRy = eye(dimy); // needed in inv(_iRy)
-}
-
-template<class sq_T>
-void Kalman<sq_T>::bayes( const vec &dt , bool evalll) {
-	it_assert_debug( dt.length()==( dimy+dimu ),"KalmanFull::bayes wrong size of dt" );
-
-	vec u = dt.get( dimy,dimy+dimu-1 );
-	vec y = dt.get( 0,dimy-1 );
-	//Time update
-	mu = A*mu + B*u;
-	//P  = A*P*A.transpose() + Q; in sq_T
-	P.mult_sym( A );
-	P+=Q;
-
-	//Data update
-	//_Ry = C*P*C.transpose() + R; in sq_T
-	_Ry.mult_sym( C, P);
-	_Ry+=R;
-
-	mat Pfull = P.to_mat();
-	
-	_Ry.inv( _iRy ); // result is in _iRy;
-	_K = Pfull*C.transpose()*(_iRy.to_mat());
-	P -= _K*C*Pfull; // P = P -KCP;
-	_yp = y-C*mu-D*u; //y prediction
-	mu += _K*( _yp );
-	
-	if (evalll==true) {
-	ll+= -0.5*(_Ry.cols()*0.79817986835811504957 \
-	+_Ry.logdet() +_iRy.qform(_yp));
-	}
-};
-
-//extern template class Kalman<ldmat>; 
-
-
-#endif // KF_H
-
-
Index: /bPF.cpp
===================================================================
--- /libPF.cpp (revision 11)
+++  (revision )
@@ -1,123 +1,0 @@
-#include "libPF.h"
-
-using namespace itpp;
-
-using std::endl;
-
-PF::PF(vec w0){w=w0/sum(w0); n=w.length();}
-
-ivec PF::resample( RESAMPLING_METHOD method ) {
-	ivec ind=zeros_i( n );
-	ivec N_babies = zeros_i( n );
-	vec cumDist = cumsum( w );
-	vec u( n );
-	int i,j,parent;
-	double u0;
-
-	switch ( method ) {
-	case MULTINOMIAL:
-		u( n - 1 ) = pow( URNG.sample(), 1.0 / n );
-		for ( i = n - 2;i >= 0;i-- ) {
-			u( i ) = u( i + 1 ) * pow( URNG.sample(), 1.0 / ( i + 1 ) );
-		}
-		break;
-	case STRATIFIED:
-		for ( i = 0;i < n;i++ ) {
-			u( i ) = ( i + URNG.sample() ) / n;
-		}
-		break;
-	case SYSTEMATIC:
-		u0 = URNG.sample();
-		for ( i = 0;i < n;i++ ) {
-			u( i ) = ( i + u0 ) / n;
-		}
-		break;
-	default:
-		it_error( "PF::resample(): Unknown resampling method" );
-	}
-	
-// 	// <<<<<<<<<<<
-// 	using std::cout;
-// 	cout << "resampling:" << endl;
-// 	cout << w <<endl;
-// 	cout << cumDist <<endl;
-// 	cout << u <<endl;
-
-	// U is now full
-	j = 0;
-	for ( i = 0;i < n;i++ ) {
-		while ( u( i ) > cumDist( j ) ) j++;
-		N_babies( j )++;
-	}
-
-// 	// <<<<<<<<<<<
-// 	cout << N_babies <<endl;
-
-	// We have assigned new babies for each Particle
-	// Now, we fill the resulting index such that:
-	// * particles with at least one baby should not move *
-	// This assures that reassignment can be done inplace;
-
-	// find the first parent;
-	parent=0; while(N_babies(parent)==0) parent++;
-	// Build index
-	for ( i = 0;i < n;i++ ) {
-		if ( N_babies( i ) > 0 ){
-			ind( i ) = i;
-			N_babies( i ) --; //this index was now replicated;
-			}
-		else {
-			// test if the parent has been fully replicated 
-			// if yes, find the next one
-			while((N_babies(parent)==0) || (N_babies(parent)==1 && parent>i) ) parent++;
-			// Replicate parent
-			ind (i) = parent;
-			N_babies( parent ) --; //this index was now replicated;
-		}
-/*			// <<<<<<<<<<<
-			cout <<ind << endl;*/
-	}
-	return ind;
-}
-
-TrivialPF::TrivialPF( mpdf &par0, mpdf &obs0, mpdf &prop0, int n0 ) : PF(1/n*ones(n)) {
-	is_proposal = true;
-	prop = &prop0;
-	par = &par0;
-	obs = &obs0;
-}
-TrivialPF::TrivialPF( mpdf &par0, mpdf &obs0,  int n0 ) : PF(1/n*ones(n)) {
-	is_proposal = false;
-	par = &par0;
-	obs = &obs0;
-}
-
-void TrivialPF::bayes( const vec &dt , bool evalll ) {
-	int i;
-	vec oldp;
-	double ll, gl, sum = 0.0;
-	Sort<double> S;
-	ivec ind, iw;
-	/*
-	//generate new samples
-	for ( i=0;i<n;i++ ) {
-		prop->evalcond( ptcls( i ), &prop_cond );
-		ptcls( i ) = prop_cond.sample();
-		gl = prop_cond.eval( ptcls( i ) );
-
-		obs.evalcond( ptcls( i ), &obs_cond );
-		ll = obs_cond.eval( dt );
-		w( i ) *= ll/gl;
-	}
-	//renormalize
-	for ( i=0;i<n;i++ ){sum+=w( i );};
-	w( i ) /=sum; //?
-	//
-	ind = resample();
-	iw = S.sort_index( 0,n-1,w ); // the first one in iw is the strongest
-
-	for ( i=0;i<n;i++ ) {
-		ptcls( i ) = ptcls( i ); //potentionally dangerous!
-	}
-	*/
-}
