Index: library/bdm/design/ctrlbase.h
===================================================================
--- library/bdm/design/ctrlbase.h (revision 586)
+++ library/bdm/design/ctrlbase.h (revision 598)
@@ -96,5 +96,5 @@
 			pre_qr.set_submatrix(0,0,s*pr);
 			pre_qr.set_submatrix(dimx+dimu+dimy, dimu+dimx, -Qy*y_req);
-			if (!qr(pre_qr,post_qr)) bdm_warning("QR in LQG unstable");
+			if (!qr(pre_qr,post_qr)){ bdm_warning("QR in LQG unstable");}
 			triu(post_qr);
 	// hn(m+1:2*m+n+r,m+1:2*m+n+r);
Index: library/bdm/base/bdmbase.h
===================================================================
--- library/bdm/base/bdmbase.h (revision 596)
+++ library/bdm/base/bdmbase.h (revision 598)
@@ -99,7 +99,13 @@
 
 	private:
+	  enum enum_dummy {dummy};
 		//! auxiliary function used in constructor
 		void init (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times);
 		int init (const string &name, int size);
+		//! Private constructor from IDs, potentially dangerous since all ids must be valid!
+		//! dummy is there to prevent confusion with RV(" string");
+		explicit RV (const ivec &ids0, enum_dummy dum): dsize(0), len(ids0.length()), ids(ids0), times(zeros_i(ids0.length())) {
+		  dsize = countsize();
+		}
 	public:
 		//! \name Constructors
@@ -171,4 +177,6 @@
 		//! Find indices of self in another rv, \return ivec of the same size as self.
 		ivec findself (const RV &rv2) const;
+		//! Find indices of self in another rv, ignore time, \return ivec of the same size as self.
+		ivec findself_ids (const RV &rv2) const;
 		//! Compare if \c rv2 is identical to this \c RV
 		bool equal (const RV &rv2) const;
@@ -192,4 +200,20 @@
 		//!@}
 
+		//! @{ \name Auxiliary functions
+		//! returns rvs with time set to 0 and removed duplicates
+		RV remove_time() const {  return RV(unique(ids), dummy);	}
+		//! return rvs with expanded delayes and sorted in the order of: \f$ [ rv_{0}, rv_{-1}, rv_{
+		RV expand_delayes() const {
+		  RV rvt=this->remove_time(); //rv at t=0
+		  RV tmp = rvt;
+		  int td = mint();
+		  for ( int i = -1; i >= td; i-- ) {
+			rvt.t ( -1 );
+			tmp.add ( rvt ); //shift u1
+		  }
+		  return tmp;
+		}
+		//!@}
+		
 		//!\name Relation to vectors
 		//!@{
@@ -200,4 +224,6 @@
 		//! Then, data can be copied via: data_of_this = cdata(ind);
 		ivec dataind (const RV &crv) const;
+		//! same as dataind but this time crv should not be complete supperset of rv.
+		ivec dataind_part (const RV &crv) const;
 		//! generate mutual indices when copying data between self and crv.
 		//! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
@@ -461,8 +487,8 @@
 		//! @{
 
-		RV _rv() const {
+		const RV& _rv() const {
 			return ep->_rv();
 		}
-		RV _rvc() {
+		const RV& _rvc() {
 			return rvc;
 		}
@@ -570,5 +596,5 @@
 		//! val-to-val link, indices of the upper val
 		ivec v2v_up;
-
+		
 	public:
 		//! Constructor
@@ -588,6 +614,12 @@
 		//! Get val for myself from val of "Up"
 		vec pushdown (const vec &val_up) {
+		  vec tmp(downsize);
+		  filldown(val_up, tmp);
+			return tmp;
+		}
+		//! Get val for vector val_down from val of "Up"
+		void filldown (const vec &val_up, vec &val_down) {
 			bdm_assert_debug (upsize == val_up.length(), "Wrong val_up");
-			return get_vec (val_up, v2v_up);
+			val_down=val_up(v2v_up);
 		}
 
@@ -598,4 +630,98 @@
 			set_subvector (val_up, v2v_up, val);
 		}
+		//! access functions
+		int _upsize(){return upsize;}
+		//! access functions
+		int _downsize(){return downsize;}
+};
+
+/*! Extension of datalink to fill only part of Down
+*/
+class datalink_part : public datalink{
+  protected:
+	//! indeces of values in vector downsize
+	ivec v2v_down;
+  public:
+	void set_connection(const RV &rv, const RV &rv_up);
+		//! Get val for vector val_down from val of "Up"
+		void filldown (const vec &val_up, vec &val_down) {
+			set_subvector(val_down, v2v_down, val_up(v2v_up));
+		}
+};
+
+/*! \brief Datalink that buffers delayed values - do not forget to call step()
+
+Up is current data, Down is their subset with possibly delayed values
+*/
+class datalink_buffered: public datalink_part{
+  protected:
+	//! History, ordered as \f$[Up_{t-1},Up_{t-2}, \ldots]\f$
+	vec history;
+	//! rv of the history
+	RV Hrv;
+	//! h2v : indeces in down 
+	ivec h2v_down;
+	//! h2v : indeces in history
+	ivec h2v_hist;
+  public:
+	
+    datalink_buffered():datalink_part(),history(0), h2v_down(0), h2v_hist(0){};
+	//! push current data to history
+	void step(const vec &val_up){
+		if (Hrv._dsize()>0){
+		  history.shift_right ( 0, Hrv._dsize() );
+		  history.set_subvector ( 0, val_up ); // write U after Drv
+		}
+	}
+		//! Get val for myself from val of "Up"
+		vec pushdown (const vec &val_up) {
+		  vec tmp(downsize);
+		  filldown(val_up, tmp);
+			return tmp;
+		}
+	
+	void filldown(const vec &val_up, vec &val_down){
+	  bdm_assert_debug(val_down.length()>=downsize, "short val_down");
+	
+	  set_subvector(val_down, v2v_down, val_up(v2v_up)); // copy direct values
+	  set_subvector(val_down, h2v_down, history(h2v_hist)); // copy delayed values
+	}
+	
+	void set_connection(const RV &rv, const RV &rv_up){
+	  // create link between up and down
+	  datalink_part::set_connection(rv, rv_up);
+	  
+	  // create rvs of history 
+	  // we can store only what we get in rv_up - everything else is removed
+	  ivec valid_ids=rv.findself_ids(rv_up);
+	  RV rv_hist = rv.subselect( find(valid_ids>=0) ); // select only rvs that are in rv_up
+	  rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0
+	  Hrv=rv_hist.subt(rv.remove_time());   // remove time 0
+	  history = zeros(Hrv._dsize());
+	  
+	  Hrv.dataind(rv,h2v_hist,h2v_down);
+	  
+	  downsize = v2v_down.length()+h2v_down.length();
+	  upsize = v2v_up.length();
+	}
+};
+
+//! buffered datalink from 2 vectors to 1
+class datalink_2to1_buffered{
+  protected:
+	datalink_buffered dl1;
+	datalink_buffered dl2;
+  public:
+	//! set connection between RVs
+	void set_connection(const RV &rv, const RV &rv_up1, const RV &rv_up2){
+	  dl1.set_connection(rv,rv_up1);
+	  dl2.set_connection(rv,rv_up2);
+	}
+	void filldown(const vec &val1, const vec &val2, vec &val_down){
+	  bdm_assert_debug(val_down.length()>=dl1._downsize()+dl2._downsize(), "short val_down");
+	  dl1.filldown(val1, val_down);
+	  dl2.filldown(val2, val_down);
+	}
+	void step(const vec &dt, const vec &ut) {dl1.step(dt); dl2.step(ut);}
 };
 
@@ -805,5 +931,5 @@
 		}
 		//!access function
-		virtual RV _drv() const {
+		virtual const RV& _drv() const {
 	//		return concat (Drv, Urv);// why!!!
 			return Drv;// why!!!
Index: library/bdm/base/datasources.h
===================================================================
--- library/bdm/base/datasources.h (revision 565)
+++ library/bdm/base/datasources.h (revision 598)
@@ -28,32 +28,128 @@
 */
 class MemDS : public DS {
-protected:
-	//! internal matrix of data
-	mat Data;
-	//! active column in the Data matrix
-	int time;
-	//!  vector of rows that are presented in Dt
-	ivec rowid;
-	//! vector of delays that are presented in Dt
-	ivec delays;
-
-public:
-	void getdata ( vec &dt );
-	void getdata ( vec &dt, const ivec &indeces );
-	void set_rvs ( RV &drv, RV &urv );
-
-	void write ( vec &ut ) {
-		bdm_error ( "MemDS::write is not supported" );
-	}
-
-	void write ( vec &ut, ivec &indices ) {
-		bdm_error ( "MemDS::write is not supported" );
-	}
-
-	void step();
-	//!Default constructor
-	MemDS () {};
-	MemDS ( mat &Dat, ivec &rowid0, ivec &delays0 );
-};
+	protected:
+		//! internal matrix of data
+		mat Data;
+		//! active column in the Data matrix
+		int time;
+		//!  vector of rows that are presented in Dt
+		ivec rowid;
+		//! vector of delays that are presented in Dt
+		ivec delays;
+
+	public:
+		void getdata ( vec &dt );
+		void getdata ( vec &dt, const ivec &indeces );
+		void set_rvs ( RV &drv, RV &urv );
+
+		void write ( vec &ut ) {
+			bdm_error ( "MemDS::write is not supported" );
+		}
+
+		void write ( vec &ut, ivec &indices ) {
+			bdm_error ( "MemDS::write is not supported" );
+		}
+
+		void step();
+		//!Default constructor
+		MemDS () {};
+		MemDS ( mat &Dat, ivec &rowid0, ivec &delays0 );
+};
+
+/*!  \brief Simulate data from a static pdf
+Trivial example of a data source, could be used for tests of some estimation algorithms. For example, simulating data from a mixture model and feeding them to mixture model estimators.
+*/
+
+class EpdfDS: public DS {
+	protected:
+		//! internal pointer to epdf from which we samplecond
+		shared_ptr<epdf> iepdf;
+		//! internal storage of data sample
+		vec dt;
+	public:
+		void step() {
+			dt=iepdf->sample();
+		}
+		void getdata ( vec &dt_out ) {
+			dt_out = dt;
+		}
+		void getdata ( vec &dt_out, const ivec &ids ) {
+			dt_out = dt ( ids );
+		}
+		const RV& _drv() {
+			return iepdf->_rv();
+		}
+
+		/*!
+		\code
+		class = "PdfDS";
+		epdf = {class="epdf_offspring", ...}// list of points
+		\endcode
+
+		*/
+		void from_setting ( const Setting &set ) {
+			iepdf=UI::build<epdf> ( set,"epdf",UI::compulsory );
+			dt = zeros ( iepdf->dimension() );
+		}
+};
+UIREGISTER ( EpdfDS );
+
+/*!  \brief Simulate data from conditional density
+Still having only one density but allowing conditioning on either input or delayed values.
+*/
+class MpdfDS :public DS {
+	protected:
+		//! internal pointer to epdf from which we samplecond
+		shared_ptr<mpdf> impdf;
+		//! internal storage of data sample
+		vec dt;
+		//! input vector
+		vec ut;
+		//! datalink between dt and ut and regressor
+		datalink_2to1_buffered rgrlink;
+		//! numeric values of regressor
+		vec rgr;
+		
+	public:
+		void step() {
+			rgrlink.filldown ( dt,ut,rgr );
+			rgrlink.step(dt,ut);//whist history
+			dt=impdf->samplecond ( rgr );
+		}
+		void getdata ( vec &dt_out ) {
+			dt_out = dt;
+		}
+		void getdata ( vec &dt_out, const ivec &ids ) {
+			dt_out = dt ( ids );
+		}
+		const RV& _drv() const {
+			return impdf->_rv();
+		}
+		void write(const vec &ut0){ut=ut0;}
+
+		/*!
+		\code
+		class = "MpdfDS";
+		mpdf = {class="epdf_offspring", ...}// list of points
+		\endcode
+
+		*/
+		void from_setting ( const Setting &set ) {
+			impdf=UI::build<mpdf> ( set,"mpdf",UI::compulsory );
+			
+			// get unique rvs form rvc
+			RV rgrv0=impdf->_rvc().remove_time();
+			// input is what in not in _rv()
+			RV urv=rgrv0.subt(impdf->_rv()); 
+			set_drv(impdf->_rv(), urv);
+			// connect input and output to rvc
+			rgrlink.set_connection(impdf->_rvc(), Drv,Urv); 
+
+			dt = zeros ( impdf->dimension() );
+			rgr = zeros ( impdf->dimensionc() );
+			ut = zeros(urv._dsize());
+		}
+};
+UIREGISTER ( MpdfDS );
 
 /*! Pseudovirtual class for reading data from files
@@ -62,22 +158,22 @@
 class FileDS: public MemDS {
 
-public:
-	void getdata ( vec &dt ) {
-		dt = Data.get_col ( time );
-	}
-
-	void getdata ( vec &dt, const ivec &indices ) {
-		vec tmp = Data.get_col ( time );
-		dt = tmp ( indices );
-	}
-
-	//! returns number of data in the file;
-	int ndat() {
-		return Data.cols();
-	}
-	//! no sense to log this type
-	void log_add ( logger &L ) {};
-	//! no sense to log this type
-	void logit ( logger &L ) {};
+	public:
+		void getdata ( vec &dt ) {
+			dt = Data.get_col ( time );
+		}
+
+		void getdata ( vec &dt, const ivec &indices ) {
+			vec tmp = Data.get_col ( time );
+			dt = tmp ( indices );
+		}
+
+		//! returns number of data in the file;
+		int ndat() {
+			return Data.cols();
+		}
+		//! no sense to log this type
+		void log_add ( logger &L ) {};
+		//! no sense to log this type
+		void logit ( logger &L ) {};
 };
 
@@ -90,19 +186,19 @@
 class ITppFileDS: public FileDS {
 
-public:
-	ITppFileDS ( const string &fname, const string &varname ) : FileDS() {
-		it_file it ( fname );
-		it << Name ( varname );
-		it >> Data;
-		time = 0;
-		//rowid and delays are ignored
-	};
-
-	ITppFileDS () : FileDS() {
-	};
-
-	void from_setting ( const Setting &set );
-
-	// TODO dodelat void to_setting( Setting &set ) const;
+	public:
+		ITppFileDS ( const string &fname, const string &varname ) : FileDS() {
+			it_file it ( fname );
+			it << Name ( varname );
+			it >> Data;
+			time = 0;
+			//rowid and delays are ignored
+		};
+
+		ITppFileDS () : FileDS() {
+		};
+
+		void from_setting ( const Setting &set );
+
+		// TODO dodelat void to_setting( Setting &set ) const;
 
 };
@@ -120,7 +216,7 @@
 class CsvFileDS: public FileDS {
 
-public:
-	//! Constructor - create DS from a CSV file.
-	CsvFileDS ( const string& fname, const string& orientation = "BY_COL" );
+	public:
+		//! Constructor - create DS from a CSV file.
+		CsvFileDS ( const string& fname, const string& orientation = "BY_COL" );
 };
 
@@ -132,137 +228,137 @@
 */
 class ArxDS : public DS {
-protected:
-	//! Rv of the regressor
-	RV Rrv;
-	//! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$
-	vec H;
-	//! (future) input
-	vec U;
-	//! temporary variable for regressor
-	vec rgr;
-	//! data link: H -> rgr
-	datalink rgrlnk;
-	//! model of Y - linear Gaussian
-	mlnorm<chmat> model;
-	//! options
-	bool opt_L_theta;
-	//! loggers
-	int L_theta;
-	int L_R;
-	int dt_size;
-public:
-	void getdata ( vec &dt ) {
-		dt = H;
-	}
-
-	void getdata ( vec &dt, const ivec &indices ) {
-		dt = H ( indices );
-	}
-
-	void write ( vec &ut ) {
-		U = ut;
-	}
-
-	void write ( vec &ut, const ivec &indices ) {
-		bdm_assert_debug ( ut.length() == indices.length(), "ArxDS" );
-		set_subvector ( U, indices, ut );
-	}
-
-	void step();
-
-	//!Default constructor
-	ArxDS ( ) {};
-	//! Set parameters of the internal model, H is maximum time delay
-	void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 ) {
-		model.set_parameters ( Th0, mu0, sqR0 );
-	};
-	//! Set
-	void set_drv ( const RV &yrv, const RV &urv, const RV &rrv ) {
-		Rrv = rrv;
-		Urv = urv;
-		dt_size = yrv._dsize() + urv._dsize();
-
-		RV drv = concat ( yrv, urv );
-		Drv = drv;
-		int td = rrv.mint();
-		H.set_size ( drv._dsize() * ( -td + 1 ) );
-		U.set_size ( Urv._dsize() );
-		for ( int i = -1; i >= td; i-- ) {
-			drv.t ( -1 );
-			Drv.add ( drv ); //shift u1
-		}
-		rgrlnk.set_connection ( rrv, Drv );
-
-		dtsize = Drv._dsize();
-		utsize = Urv._dsize();
-	}
-	//! set options from a string
-	void set_options ( const string &s ) {
-		opt_L_theta = ( s.find ( "L_theta" ) != string::npos );
-	};
-	virtual void log_add ( logger &L ) {
-		//DS::log_add ( L ); too long!!
-		L_dt = L.add ( Drv ( 0, dt_size ), "" );
-		L_ut = L.add ( Urv, "" );
-
-		mat &A = model._A();
-		mat R = model._R();
-		if ( opt_L_theta ) {
-			L_theta = L.add ( RV ( "{th }", vec_1 ( A.rows() * A.cols() ) ), "t" );
-		}
-		if ( opt_L_theta ) {
-			L_R = L.add ( RV ( "{R }", vec_1 ( R.rows() * R.cols() ) ), "r" );
-		}
-	}
-	virtual void logit ( logger &L ) {
-		//DS::logit ( L );
-		L.logit ( L_dt, H.left ( dt_size ) );
-		L.logit ( L_ut, U );
-
-		mat &A = model._A();
-		mat R = model._R();
-		if ( opt_L_theta ) {
-			L.logit ( L_theta, vec ( A._data(), A.rows() *A.cols() ) );
-		};
-		if ( opt_L_theta ) {
-			L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) );
-		};
-	}
-
-	// TODO dokumentace - aktualizovat
-	/*! UI for ArxDS using factorized description!
-
-	The ArxDS is constructed from a structure with fields:
-	\code
-	system = {
-		type = "ArxDS";
-		// description of y variables
-		y = {type="rv"; names=["y", "u"];};
-		// description of u variable
-		u = {type="rv"; names=[];}
-		// description of regressor
-		rgr = {type="rv";
-			names = ["y","y","y","u"];
-			times = [-1, -2, -3, -1];
-		}
-
-		// theta
-		theta = [0.8, -0.3, 0.4, 1.0,
-				 0.0, 0.0, 0.0, 0.0];
-		// offset (optional)
-		offset = [0.0, 0.0];
-		//variance
-		r = [0.1, 0.0,
-			 0.0, 1.0];
-		//options: L_theta = log value of theta,
-		opt = "L_theta";
-	};
-	\endcode
-
-	Result is ARX data source offering with full history as Drv.
-	*/
-	void from_setting ( const Setting &set );
-
-	// TODO dodelat void to_setting( Setting &set ) const;
+	protected:
+		//! Rv of the regressor
+		RV Rrv;
+		//! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$
+		vec H;
+		//! (future) input
+		vec U;
+		//! temporary variable for regressor
+		vec rgr;
+		//! data link: H -> rgr
+		datalink rgrlnk;
+		//! model of Y - linear Gaussian
+		mlnorm<chmat> model;
+		//! options
+		bool opt_L_theta;
+		//! loggers
+		int L_theta;
+		int L_R;
+		int dt_size;
+	public:
+		void getdata ( vec &dt ) {
+			dt = H;
+		}
+
+		void getdata ( vec &dt, const ivec &indices ) {
+			dt = H ( indices );
+		}
+
+		void write ( vec &ut ) {
+			U = ut;
+		}
+
+		void write ( vec &ut, const ivec &indices ) {
+			bdm_assert_debug ( ut.length() == indices.length(), "ArxDS" );
+			set_subvector ( U, indices, ut );
+		}
+
+		void step();
+
+		//!Default constructor
+		ArxDS ( ) {};
+		//! Set parameters of the internal model, H is maximum time delay
+		void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 ) {
+			model.set_parameters ( Th0, mu0, sqR0 );
+		};
+		//! Set
+		void set_drv ( const RV &yrv, const RV &urv, const RV &rrv ) {
+			Rrv = rrv;
+			Urv = urv;
+			dt_size = yrv._dsize() + urv._dsize();
+
+			RV drv = concat ( yrv, urv );
+			Drv = drv;
+			int td = rrv.mint();
+			H.set_size ( drv._dsize() * ( -td + 1 ) );
+			U.set_size ( Urv._dsize() );
+			for ( int i = -1; i >= td; i-- ) {
+				drv.t ( -1 );
+				Drv.add ( drv ); //shift u1
+			}
+			rgrlnk.set_connection ( rrv, Drv );
+
+			dtsize = Drv._dsize();
+			utsize = Urv._dsize();
+		}
+		//! set options from a string
+		void set_options ( const string &s ) {
+			opt_L_theta = ( s.find ( "L_theta" ) != string::npos );
+		};
+		virtual void log_add ( logger &L ) {
+			//DS::log_add ( L ); too long!!
+			L_dt = L.add ( Drv ( 0, dt_size ), "" );
+			L_ut = L.add ( Urv, "" );
+
+			mat &A = model._A();
+			mat R = model._R();
+			if ( opt_L_theta ) {
+				L_theta = L.add ( RV ( "{th }", vec_1 ( A.rows() * A.cols() ) ), "t" );
+			}
+			if ( opt_L_theta ) {
+				L_R = L.add ( RV ( "{R }", vec_1 ( R.rows() * R.cols() ) ), "r" );
+			}
+		}
+		virtual void logit ( logger &L ) {
+			//DS::logit ( L );
+			L.logit ( L_dt, H.left ( dt_size ) );
+			L.logit ( L_ut, U );
+
+			mat &A = model._A();
+			mat R = model._R();
+			if ( opt_L_theta ) {
+				L.logit ( L_theta, vec ( A._data(), A.rows() *A.cols() ) );
+			};
+			if ( opt_L_theta ) {
+				L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) );
+			};
+		}
+
+		// TODO dokumentace - aktualizovat
+		/*! UI for ArxDS using factorized description!
+
+		The ArxDS is constructed from a structure with fields:
+		\code
+		system = {
+			type = "ArxDS";
+			// description of y variables
+			y = {type="rv"; names=["y", "u"];};
+			// description of u variable
+			u = {type="rv"; names=[];}
+			// description of regressor
+			rgr = {type="rv";
+				names = ["y","y","y","u"];
+				times = [-1, -2, -3, -1];
+			}
+
+			// theta
+			theta = [0.8, -0.3, 0.4, 1.0,
+					 0.0, 0.0, 0.0, 0.0];
+			// offset (optional)
+			offset = [0.0, 0.0];
+			//variance
+			r = [0.1, 0.0,
+				 0.0, 1.0];
+			//options: L_theta = log value of theta,
+			opt = "L_theta";
+		};
+		\endcode
+
+		Result is ARX data source offering with full history as Drv.
+		*/
+		void from_setting ( const Setting &set );
+
+		// TODO dodelat void to_setting( Setting &set ) const;
 };
 
@@ -271,82 +367,82 @@
 
 class stateDS : public DS {
-private:
-	//!conditional pdf of the state evolution \f$ f(x_t|x_{t-1}) \f$
-	shared_ptr<mpdf> IM;
-
-	//!conditional pdf of the observations \f$ f(d_t|x_t) \f$
-	shared_ptr<mpdf> OM;
-
-protected:
-	//! result storage
-	vec dt;
-	//! state storage
-	vec xt;
-	//! input storage
-	vec ut;
-	//! Logger
-	int L_xt;
-
-public:
-	void getdata ( vec &dt0 ) {
-		dt0 = dt;
-	}
-
-	void getdata ( vec &dt0, const ivec &indices ) {
-		dt0 = dt ( indices );
-	}
-
-	stateDS ( const shared_ptr<mpdf> &IM0, const shared_ptr<mpdf> &OM0, int usize ) : IM ( IM0 ), OM ( OM0 ),
-		dt ( OM0->dimension() ), xt ( IM0->dimension() ),
-		ut ( usize ), L_xt(0) { }
-
-	stateDS() : L_xt(0) { }
-
-	virtual void step() {
-		xt = IM->samplecond ( concat ( xt, ut ) );
-		dt = OM->samplecond ( concat ( xt, ut ) );
-	}
-
-	virtual void log_add ( logger &L ) {
-		DS::log_add ( L );
-		L_xt = L.add ( IM->_rv(), "true" );
-	}
-	virtual void logit ( logger &L ) {
-		DS::logit ( L );
-		L.logit ( L_xt, xt );
-	}
-
-	/*! UI for stateDS
-
-	The DS is constructed from a structure with fields:
-	\code
-	system = {
-		type = "stateDS";
-		//Internal model
-		IM = { type = "mpdf"; //<-- valid offspring! e.g. "mlnorm"
-			rv = { //description of x_t
-				names=["name1",...];
-				sizes=[2,1]; // optional default=[1,1...];
-				times=[0,0]; // optional default=[0,0...];
-				}
-			rvu= { //description of  u_t
-				//optional default=empty
-				}
-
-			// remaining fields depending on the chosen type
-			};
-		//Observation model
-		OM = { type = "mpdf-offspring";
-			rv = {}; //description of d_t
-			rvu = {type="internal", path="system.IM.rvu"}; //description of u_t
-
-			//remaining fields
-		}
-	};
-	\endcode
-	*/
-	void from_setting ( const Setting &set );
-
-	// TODO dodelat void to_setting( Setting &set ) const;
+	private:
+		//!conditional pdf of the state evolution \f$ f(x_t|x_{t-1}) \f$
+		shared_ptr<mpdf> IM;
+
+		//!conditional pdf of the observations \f$ f(d_t|x_t) \f$
+		shared_ptr<mpdf> OM;
+
+	protected:
+		//! result storage
+		vec dt;
+		//! state storage
+		vec xt;
+		//! input storage
+		vec ut;
+		//! Logger
+		int L_xt;
+
+	public:
+		void getdata ( vec &dt0 ) {
+			dt0 = dt;
+		}
+
+		void getdata ( vec &dt0, const ivec &indices ) {
+			dt0 = dt ( indices );
+		}
+
+		stateDS ( const shared_ptr<mpdf> &IM0, const shared_ptr<mpdf> &OM0, int usize ) : IM ( IM0 ), OM ( OM0 ),
+				dt ( OM0->dimension() ), xt ( IM0->dimension() ),
+				ut ( usize ), L_xt ( 0 ) { }
+
+		stateDS() : L_xt ( 0 ) { }
+
+		virtual void step() {
+			xt = IM->samplecond ( concat ( xt, ut ) );
+			dt = OM->samplecond ( concat ( xt, ut ) );
+		}
+
+		virtual void log_add ( logger &L ) {
+			DS::log_add ( L );
+			L_xt = L.add ( IM->_rv(), "true" );
+		}
+		virtual void logit ( logger &L ) {
+			DS::logit ( L );
+			L.logit ( L_xt, xt );
+		}
+
+		/*! UI for stateDS
+
+		The DS is constructed from a structure with fields:
+		\code
+		system = {
+			type = "stateDS";
+			//Internal model
+			IM = { type = "mpdf"; //<-- valid offspring! e.g. "mlnorm"
+				rv = { //description of x_t
+					names=["name1",...];
+					sizes=[2,1]; // optional default=[1,1...];
+					times=[0,0]; // optional default=[0,0...];
+					}
+				rvu= { //description of  u_t
+					//optional default=empty
+					}
+
+				// remaining fields depending on the chosen type
+				};
+			//Observation model
+			OM = { type = "mpdf-offspring";
+				rv = {}; //description of d_t
+				rvu = {type="internal", path="system.IM.rvu"}; //description of u_t
+
+				//remaining fields
+			}
+		};
+		\endcode
+		*/
+		void from_setting ( const Setting &set );
+
+		// TODO dodelat void to_setting( Setting &set ) const;
 
 };
Index: library/bdm/base/bdmbase.cpp
===================================================================
--- library/bdm/base/bdmbase.cpp (revision 565)
+++ library/bdm/base/bdmbase.cpp (revision 598)
@@ -191,8 +191,7 @@
 void datalink::set_connection (const RV &rv, const RV &rv_up) {
 	downsize = rv._dsize();
-	upsize = rv_up._dsize();
-	v2v_up = rv.dataind (rv_up);
-
-	bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
+	  upsize = rv_up._dsize();
+	  v2v_up = rv.dataind (rv_up);
+	  bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
 }
 
@@ -201,6 +200,11 @@
 	upsize = us;
 	v2v_up = upind;
-
-	bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
+    bdm_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
+}
+
+void datalink_part::set_connection (const RV &rv, const RV &rv_up) {
+	 rv.dataind (rv_up, v2v_down, v2v_up);
+	 downsize = v2v_down.length();
+	 upsize = v2v_up.length();
 }
 
@@ -321,4 +325,18 @@
 }
 
+ivec RV::findself_ids ( const RV &rv2 ) const {
+	int i, j;
+	ivec tmp = -ones_i ( len );
+	for ( i = 0; i < len; i++ ) {
+		for ( j = 0; j < rv2.length(); j++ ) {
+			if ( ( ids ( i ) == rv2.ids ( j ) ) ) {
+				tmp ( i ) = j;
+				break;
+			}
+		}
+	}
+	return tmp;
+}
+
 void RV::from_setting ( const Setting &set ) {
 	Array<string> A;
Index: library/tests/datasource_test.cpp
===================================================================
--- library/tests/datasource_test.cpp (revision 598)
+++ library/tests/datasource_test.cpp (revision 598)
@@ -0,0 +1,43 @@
+#include "../bdm/base/datasources.h"
+#include "mat_checks.h"
+#include "UnitTest++.h"
+
+using namespace bdm;
+
+TEST ( test_EpdfDS ) {
+	UIFile uif("epdfds.cfg");
+
+	shared_ptr<EpdfDS> ds=UI::build<EpdfDS>(uif, "ds");
+	
+	vec mean = zeros(2);
+	vec dt = zeros(2);
+	for (int i=0; i<100; i++){
+		ds->step();
+		ds->getdata(dt);
+		mean += dt;
+	}
+	CHECK_CLOSE_EX (  vec_2(1.0,2.0) ,mean/100, 1e-2);
+}
+
+TEST ( test_MpdfDS ) {
+	UIFile uif("mpdfds.cfg");
+
+	shared_ptr<MpdfDS> ds=UI::build<MpdfDS>(uif, "ds", UI::compulsory);
+	
+	vec dt = zeros(1);
+	vec ut="1.0";
+	ds->write(ut);
+	for (int i=0; i<100; i++){
+		ds->step();
+	}
+	ds->getdata(dt);
+	CHECK_CLOSE ( -0.2 , dt(0), 1e-5);
+
+	ut="2.0";
+	ds->write(ut);
+	for (int i=0; i<100; i++){
+		ds->step();
+	}
+	ds->getdata(dt);
+	CHECK_CLOSE ( -0.4 , dt(0), 1e-5);
+}
Index: library/tests/epdfds.cfg
===================================================================
--- library/tests/epdfds.cfg (revision 598)
+++ library/tests/epdfds.cfg (revision 598)
@@ -0,0 +1,12 @@
+ds={
+class = "EpdfDS";
+epdf = {
+  class = "enorm<ldmat>";
+  mu = ( "matrix", 2, 1, [ 1.0, 2.0 ] );
+  R =  ( "matrix", 2, 2, [ 1e-4, 0.0, 0.0,  2e-4 ] );
+  rv = {
+    class = "RV";
+    names = ( "a", "b" );
+  };
+};
+};
Index: library/tests/datalink_test.cpp
===================================================================
--- library/tests/datalink_test.cpp (revision 545)
+++ library/tests/datalink_test.cpp (revision 598)
@@ -87,2 +87,25 @@
 	}
 }
+
+TEST ( test_datalink_buffered ) {
+	RV a = RV ( "{dl_a }" );
+	RV b = RV ( "{dlb }" );
+	RV ab=concat(a,b);
+	RV aa = a;
+	a.t(-1);
+	aa.add(a); 
+
+	
+	datalink_buffered dl;
+	dl.set_connection ( concat(aa,b), ab);
+
+	vec val_up ( "1 2" );
+	dl.step(val_up);
+	val_up="3 4";
+	
+	
+	vec p = dl.pushdown ( val_up );
+	vec exp_p = " 3, 1, 4 ";
+	CHECK_EQUAL ( exp_p, p );
+
+}
Index: library/tests/mpdfds.cfg
===================================================================
--- library/tests/mpdfds.cfg (revision 598)
+++ library/tests/mpdfds.cfg (revision 598)
@@ -0,0 +1,18 @@
+ds={
+class = "MpdfDS";
+mpdf = {
+  class = "mlnorm<ldmat>";
+  A = ( "matrix", 1, 2, [ 0.5, -0.1 ] );
+  R = ( "matrix", 1, 1, [ 1e-10 ] );
+  const =  ( "matrix", 1, 1, [ 0.0 ] );
+  rv = {
+    class = "RV";
+    names = ( "a" );
+  };
+  rvc = {
+    class = "RV";
+    names = ( "a", "b" );
+	 times = [-1 , 0];
+  };
+};
+};
Index: library/tests/CMakeLists.txt
===================================================================
--- library/tests/CMakeLists.txt (revision 586)
+++ library/tests/CMakeLists.txt (revision 598)
@@ -39,5 +39,5 @@
 
 # using UnitTest++
-add_executable(testsuite datalink_test.cpp egiw_test.cpp emix_test.cpp epdf_test.cpp logger_test.cpp LQG_test.cpp merger_test.cpp 
+add_executable(testsuite datalink_test.cpp  datasource_test.cpp egiw_test.cpp emix_test.cpp epdf_test.cpp logger_test.cpp LQG_test.cpp merger_test.cpp 
 	mpdf_test.cpp randun_test.cpp rectangular_support_test.cpp rv_test.cpp shared_ptr_test.cpp square_mat_test.cpp testsuite.cpp user_info_test.cpp)
 target_link_libraries(testsuite testutil unittest)
Index: library/tests/epdf_harness.cpp
===================================================================
--- library/tests/epdf_harness.cpp (revision 536)
+++ library/tests/epdf_harness.cpp (revision 598)
@@ -81,11 +81,9 @@
 		// test of pdflog at zero
 		vec zero ( 0 );
-		vec zeron ( hepdf->dimension() );
-		for ( int i = 0; i < zeron.size(); ++i ) {
-			zeron ( i ) = 0;
-		}
+		vec zeron=zeros ( hepdf->dimension() );
 
 		double lpz = hepdf->evallog ( zeron );
-		CHECK_CLOSE_EX ( lpz, mEp.evallogcond ( zeron, zero ), tolerance );
+		double lpzc=mEp.evallogcond ( zeron, zero );
+		CHECK_CLOSE_EX ( lpz, lpzc, tolerance );
 
 		vec lpzv(1);
