Index: library/bdm/estim/arx.h
===================================================================
--- library/bdm/estim/arx.h (revision 700)
+++ library/bdm/estim/arx.h (revision 723)
@@ -98,5 +98,9 @@
 	}
 	//! conditional version of the predictor
-	mlnorm<ldmat>* predictor() const;
+	template<class sq_T>
+	shared_ptr<mlnorm<sq_T> > ml_predictor() const;
+	//! fast version of predicto 
+	template<class sq_T>
+	void ml_predictor_update(mlnorm<sq_T> &pred) const;
 	mlstudent* predictor_student() const;
 	//! Brute force structure estimation.\return indeces of accepted regressors.
@@ -188,4 +192,35 @@
 };
 UIREGISTER(ARXfrg);
+
+
+////////////////////
+template<class sq_T>
+shared_ptr< mlnorm<sq_T> > ARX::ml_predictor ( ) const {
+	shared_ptr< mlnorm<sq_T> > tmp = new mlnorm<sq_T> ( );
+	tmp->set_rv(yrv);
+	tmp->set_rvc(_rvc());
+	
+	ml_predictor_update(*tmp);
+	tmp->validate();
+	return tmp;
+}
+
+template<class sq_T>
+void ARX::ml_predictor_update(mlnorm<sq_T> &pred) const {
+	mat mu ( dimy, posterior()._V().rows() - dimy );
+	mat R ( dimy, dimy );
+	
+	est.mean_mat ( mu, R ); //mu =
+	mu = mu.T();
+	//correction for student-t  -- TODO check if correct!!
+	
+	if ( have_constant) { // constant term
+		//Assume the constant term is the last one:
+		pred.set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), sq_T ( R ) );
+	} else {
+		pred.set_parameters ( mu, zeros ( dimy ), sq_T( R ) );
+	}
+}
+
 };
 #endif // AR_H
Index: library/bdm/estim/kalman.h
===================================================================
--- library/bdm/estim/kalman.h (revision 703)
+++ library/bdm/estim/kalman.h (revision 723)
@@ -27,7 +27,7 @@
  * \brief Basic elements of linear state-space model
 
-Parameter evolution model:\f[ x_t = A x_{t-1} + B u_t + Q^{1/2} e_t \f]
-Observation model: \f[ y_t = C x_{t-1} + C u_t + Q^{1/2} w_t. \f]
-Where $e_t$ and $w_t$ are independent vectors Normal(0,1)-distributed disturbances.
+Parameter evolution model:\f[ x_{t+1} = A x_{t} + B u_t + Q^{1/2} e_t \f]
+Observation model: \f[ y_t = C x_{t} + C u_t + R^{1/2} w_t. \f]
+Where $e_t$ and $w_t$ are mutually independent vectors of Normal(0,1)-distributed disturbances.
  */
 template<class sq_T>
@@ -61,5 +61,5 @@
 			UI::get (C, set, "C", UI::compulsory);
 			UI::get (D, set, "D", UI::compulsory);
-			mat Qtm, Rtm;
+			mat Qtm, Rtm; // full matrices
 			if(!UI::get(Qtm, set, "Q", UI::optional)){
 				vec dq;
@@ -72,5 +72,5 @@
 				Rtm=diag(dr);
 			}
-			R=Rtm; // automatic conversion
+			R=Rtm; // automatic conversion to square-root form
 			Q=Qtm; 
 			
@@ -519,17 +519,23 @@
 
 */
-class StateFromARX: public StateSpace<fsqmat>{
+class StateFromARX: public StateSpace<chmat>{
 	protected:
 		//! remember connection from theta ->A
 		datalink_part th2A;
-		//! xrv
-		RV xrv;
+		//! remember connection from theta ->B
+		datalink_part th2B;
 		//!function adds n diagonal elements from given starting point r,c
 		void diagonal_part(mat &A, int r, int c, int n){
 			for(int i=0; i<n; i++){A(r,c)=1.0; r++;c++;}
 		};
+		//! similar to ARX.have_constant
+		bool have_constant;
 	public:
 		//! set up this object to match given mlnorm
-		void connect_mlnorm(const mlnorm<fsqmat> &ml){
+		//! Note that state-space and common mpdf use different meaning of \f$ _t \f$ in \f$ u_t \f$. 
+		//!While mlnorm typically assumes that \f$ u_t \rightarrow y_t \f$ in state space it is \f$ u_{t-1} \rightarrow y_t \f$
+		//! For consequences in notation of internal variable xt see arx2statespace_notes.lyx.
+		void connect_mlnorm(const mlnorm<chmat> &ml, RV &xrv, RV &urv){
+
 			//get ids of yrv				
 			const RV &yrv = ml._rv();
@@ -537,14 +543,17 @@
 			const RV &rgr = ml._rvc();
 			RV rgr0 = rgr.remove_time();
-			RV urv = rgr0.subt(yrv); 
+			urv = rgr0.subt(yrv); 
 						
 			// create names for state variables
 			xrv=yrv; 
 			
-			int y_multiplicity = rgr.mint(yrv);
-			int y_block_size = yrv.length()*(y_multiplicity+1); // current yt + all delayed yts
-			for (int m=0;m<y_multiplicity;m++){ 
+			int y_multiplicity = -rgr.mint(yrv);
+			int y_block_size = yrv.length()*(y_multiplicity); // current yt + all delayed yts
+			for (int m=0;m<y_multiplicity - 1;m++){ // ========= -1 is important see arx2statespace_notes
 				xrv.add(yrv.copy_t(-m-1)); //add delayed yt
 			}
+			//! temporary RV for connection to ml.rvc, since notation of xrv and ml.rvc does not match
+			RV xrv_ml = xrv.copy_t(-1);
+			
 			// add regressors
 			ivec u_block_sizes(urv.length()); // no_blocks = yt + unique rgr
@@ -552,13 +561,21 @@
 				RV R=urv.subselect(vec_1(r)); //non-delayed element of rgr
 				int r_size=urv.size(r);
-				int r_multiplicity=rgr.mint(R);
+				int r_multiplicity=-rgr.mint(R);
 				u_block_sizes(r)= r_size*r_multiplicity ;
 				for(int m=0;m<r_multiplicity;m++){ 
 					xrv.add(R.copy_t(-m-1)); //add delayed yt
+					xrv_ml.add(R.copy_t(-m-1)); //add delayed yt
 				}
 			}
-			
+			// add constant
+			if (any(ml._mu_const()!=0.0)){
+				have_constant=true;
+				xrv.add(RV("bdm_reserved_constant_one",1));
+			} else {have_constant=false;}
+				
+		
 			// get mapp
-			th2A.set_connection(xrv, ml._rvc());
+			th2A.set_connection(xrv_ml, ml._rvc());
+			th2B.set_connection(urv, ml._rvc());
 			
 			//set matrix sizes
@@ -577,19 +594,37 @@
 			this->C=zeros(yrv._dsize(), xrv._dsize());
 			this->C.set_submatrix(0,0,eye(yrv._dsize()));
-				this->D=zeros(yrv._dsize(),urv._dsize());
-				this->R = zeros(yrv._dsize(),yrv._dsize());
-				this->Q = zeros(xrv._dsize(),xrv._dsize());
-				// Q is set by update
-								
-				update_from(ml);
-				validate();
+			this->D=zeros(yrv._dsize(),urv._dsize());
+			this->R.setCh(zeros(yrv._dsize(),yrv._dsize()));
+			this->Q.setCh(zeros(xrv._dsize(),xrv._dsize()));
+			// Q is set by update
+							
+			update_from(ml);
+			validate();
 		};
 		//! fast function to update parameters from ml - not checked for compatibility!!
-		void update_from(const mlnorm<fsqmat> &ml){
-			
-			vec theta = ml._A().get_row(0); // this 
-			this->Q._M().set_submatrix(0,0,ml._R());
-			
-		}
+		void update_from(const mlnorm<chmat> &ml){
+
+			vec Arow=zeros(A.cols());
+			vec Brow=zeros(B.cols());
+			//  ROW- WISE EVALUATION =====
+			for(int i=0;i<ml._rv()._dsize(); i++){ 
+				
+				vec theta = ml._A().get_row(i); 
+				
+				th2A.filldown(theta,Arow);
+				if (have_constant){
+					// constant is always at the end no need for datalink
+					Arow(A.cols()-1) = ml._mu_const()(i);
+				}
+				this->A.set_row(i,Arow);
+				
+				th2B.filldown(theta,Brow);
+				this->B.set_row(i,Brow);
+			}
+			this->Q._Ch().set_submatrix(0,0,ml.__R()._Ch());
+			
+		};
+		//! access function
+		bool _have_constant() const {return have_constant;}
 };
 
Index: library/bdm/estim/arx.cpp
===================================================================
--- library/bdm/estim/arx.cpp (revision 700)
+++ library/bdm/estim/arx.cpp (revision 723)
@@ -98,22 +98,4 @@
 }
 
-mlnorm<ldmat>* ARX::predictor ( ) const {
-	mat mu ( dimy, posterior()._V().rows() - dimy );
-	mat R ( dimy, dimy );
-	mlnorm<ldmat>* tmp;
-	tmp = new mlnorm<ldmat> ( );
-
-	est.mean_mat ( mu, R ); //mu =
-	mu = mu.T();
-	//correction for student-t  -- TODO check if correct!!
-
-	if ( have_constant) { // constant term
-		//Assume the constant term is the last one:
-		tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) );
-	} else {
-		tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ) );
-	}
-	return tmp;
-}
 
 mlstudent* ARX::predictor_student ( ) const {
Index: library/bdm/itpp_ext.h
===================================================================
--- library/bdm/itpp_ext.h (revision 712)
+++ library/bdm/itpp_ext.h (revision 723)
@@ -45,5 +45,5 @@
 }
 
-const double inf = std::numeric_limits<double>::infinity();
+static const double inf = std::numeric_limits<double>::infinity();
 
 //#if 0
Index: library/bdm/math/chmat.h
===================================================================
--- library/bdm/math/chmat.h (revision 660)
+++ library/bdm/math/chmat.h (revision 723)
@@ -106,5 +106,16 @@
 		copy_vector ( dim*dim, chQ._data(), Ch._data() );
 	}
-
+	//! access function
+	void setCh ( const chmat &Ch0 ) {
+		Ch = Ch0._Ch();
+		dim=Ch0.rows();
+	}
+	//! access function
+	void setCh ( const mat &Ch0 ) {
+		//TODO check if Ch0 is OK!!!
+		Ch = Ch0;
+		dim=Ch0.rows();
+	}
+	
 	//! Access functions
 	void setD ( const vec &nD, int i ) {
Index: library/bdm/design/ctrlbase.cpp
===================================================================
--- library/bdm/design/ctrlbase.cpp (revision 706)
+++ library/bdm/design/ctrlbase.cpp (revision 723)
@@ -3,11 +3,6 @@
 namespace bdm{
 
-void LQG::set_system(shared_ptr<StateSpace<fsqmat> > S0){
+void LQG::set_system(shared_ptr<StateSpace<chmat> > S0){
 	S = S0;
-	dimx= S->_A().rows();
-	dimy= S->_C().rows();
-	dimu= S->_B().cols();
-	pr=zeros(dimx+dimu+dimy,  dimu+dimx+dimu+dimy);
-	pr.set_submatrix(dimx, dimu+dimx, eye(dimu+dimy));
 }
 
@@ -17,37 +12,30 @@
 	y_req=y_req0;
 	horizon = horizon0;
-	validate();
-	initialize();
 }
 
 void LQG::validate() {
+	// set properties from S
+	dimx= S->_A().rows();
+	dimy= S->_C().rows();
+	dimu= S->_B().cols();
+	pr=zeros(dimx+dimu+dimy,  dimu+dimx+dimu+dimy);
+	pr.set_submatrix(dimx, dimu+dimx, eye(dimu+dimy));
+	
+	//penalization
 	bdm_assert ( Qy.cols() == dimy, "LQG: wrong dimensions of Qy " );
 	bdm_assert ( Qu.cols() == dimu, "LQG: wrong dimensions of Qu " );
 	bdm_assert ( y_req.length() == dimy, "LQG: wrong dimensions of y_req " );
+	
+	qux=zeros(dimu,dimx+2*dimu+dimy);     
+	qyx=zeros(dimy, dimx+dimy+dimu);
+	
+	//
+	initial_belmann();
+	// parts of QR
+	post_qr = zeros(pre_qr.rows(), pre_qr.cols());
+	
+	update_system();
 }
 
-void LQG::initialize(){
-	// set parameter matrix
-	pr.set_submatrix(0,0,S->_B());
-	pr.set_submatrix(0,dimu, S->_A());
-	
-	//penalization
-	qux=zeros(dimu,dimx+2*dimu+dimy);     
-	qux.set_submatrix(0,0,Qu);
-	qux.set_submatrix(0,dimx+dimu+dimy,Qu);
-
-	qyx=zeros(dimy, dimx+dimy+dimu);
-	qyx.set_submatrix(0,0,S->_C());
-	qyx.set_submatrix(0,dimx,-eye(dimy));
-	
-	//
-	s=1e-5*eye(dimx+dimu+dimy);
-	// parts of QR
-	hqy=Qy*qyx*pr;
-	
-	// pre_qr
-	pre_qr = concat_vertical(s*pr, concat_vertical(hqy, qux));
-	post_qr = zeros(pre_qr.rows(), pre_qr.cols());
-}
 
 }
Index: library/bdm/design/arx_ctrl.h
===================================================================
--- library/bdm/design/arx_ctrl.h (revision 723)
+++ library/bdm/design/arx_ctrl.h (revision 723)
@@ -0,0 +1,140 @@
+/*!
+  \file
+  \brief Controllers for linear Gaussian systems
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#include "../base/bdmbase.h"
+#include "ctrlbase.h"
+#include "../estim/arx.h"
+#include "../estim/kalman.h"
+//#include <../applications/pmsm/simulator_zdenek/ekf_example/pmsm_mod.h>
+
+namespace bdm{
+
+	//! Controller using ARX model for estimation and LQG designer for control
+	class LQG_ARX : public Controller{
+		protected:
+		//! Internal ARX estimator
+			shared_ptr<ARX> ar;
+		//! Internal LQG designer
+			LQG lq;
+			//! Intermediate StateSpace model
+			shared_ptr<StateFromARX> Stsp;
+			//! AR predictor
+			shared_ptr<mlnorm<chmat> > pred;
+			//! datalink from rvc to ar.bayes
+			datalink rvc2ar_y;
+			//! datalink from rvc to ar.bayes
+			datalink_buffered rvc2ar_cond;
+			
+			//! flag to use iterations spread in time (ist)
+			bool ist;
+			//! flag to use windsurfer approach
+			bool windsurfer;
+		public:
+			//! Default constructor
+			LQG_ARX(): ar(), lq() { }
+
+			//! adaptation is to store arx estimates in stsp
+			void adapt(const vec &data){
+				ar->bayes(rvc2ar_y.pushdown(data), rvc2ar_cond.pushdown(data));
+				rvc2ar_cond.store_data(data);
+			}
+			void redesign(){
+				ar->ml_predictor_update(*pred);
+				Stsp->update_from(*pred);
+				if (windsurfer) {
+					mat Ry=pred->_R();
+					lq.set_control_Qy(inv(Ry));
+				}
+				if (!ist) {
+					lq.initial_belmann();
+				}
+				lq.redesign();
+			}
+			vec ctrlaction(const vec &cond) const{
+				//cond is xt + ut
+				vec state=cond.left(Stsp->_A().rows());
+				if (Stsp->_have_constant()){
+					state(state.length()-1) = 1;
+				}
+				return lq.ctrlaction(state,cond.right(Stsp->_B().cols()));
+			}
+			//! 
+	//! LQG is defined by quadratic loss function
+	/*! \f[ L(y,u) = (y-y_{req})'Q_y (y-y_{req}) + (u-u_{req})' Q_u (u-u_{req}) \f]
+			expected input
+			\code
+		{ class="LQG";
+			arx = {class="ARX", ...} // internal arx model, see 
+			Qy = ("matrix", ...); 
+			Qu = ("matrix", ...); 
+			yreq = [];               // requested output
+			ureq = [];               // requested input value
+		}
+			\endcode
+	 */
+			void from_setting(const Setting &set){
+				ar=UI::build<ARX> (set, "ARX",UI::compulsory);
+
+				mat Qu;
+				mat Qy;
+				
+				UI::get(Qu, set, "Qu", UI::compulsory);
+				UI::get(Qy, set, "Qy", UI::compulsory);
+		
+				vec y_req;
+				if (!UI::get(y_req, set, "yreq", UI::optional))
+					y_req = zeros(ar->_yrv()._dsize());
+				
+				int horizon;
+				UI::get(horizon, set, "horizon", UI::compulsory);
+				lq.set_control_parameters(Qy,Qu,y_req,horizon);
+				
+				int wind;
+				if (UI::get(wind, set, "windsurfer", UI::optional)){
+					windsurfer=wind>0;
+				} else {
+					windsurfer=false;
+				};
+				int ist_;
+				if (UI::get(ist_, set, "ist", UI::optional)){
+					ist=ist_>0;
+				} else {
+					ist=false;
+				};
+				validate();
+			}
+
+		void validate() {
+			// ar is valid
+			pred = ar->ml_predictor<chmat>();
+			Stsp = new StateFromARX;
+			RV xrv;
+			RV urvm; //old ut
+			Stsp->connect_mlnorm(*pred, xrv, urvm);
+			lq.set_system(Stsp);
+			lq.validate();
+			
+			rvc =concat(xrv, urvm);
+			rvc2ar_y.set_connection(ar->_yrv(), rvc);
+			rvc2ar_cond.set_connection(ar->_rvc(), rvc);
+			//datalink from ARX to rvc
+		}
+		void log_register (logger &L, const string &prefix ) {
+			ar->log_register(L,prefix);
+		}
+		void log_write ( ) const { 
+			ar->log_write();
+		}
+		
+	};
+	UIREGISTER(LQG_ARX);
+} // namespace
Index: library/bdm/design/ctrlbase.h
===================================================================
--- library/bdm/design/ctrlbase.h (revision 706)
+++ library/bdm/design/ctrlbase.h (revision 723)
@@ -14,122 +14,142 @@
 #include "../estim/kalman.h"
 
-namespace bdm{
+namespace bdm {
 
-	//! Base class for adaptive controllers 
-	//! The base class is, however, non-adaptive, method \c adapt() is empty.
-	//! \note advanced Controllers will probably include estimator as their internal attribute (e.g. dual controllers)
-	class Controller : public root {
-		protected:
-			//! identifier of the designed action;
-			RV rv;
-			//! identifier of the conditioning variables;
-			RV rvc;
-		public:
-			//! function processing new observations and adapting control strategy accordingly
-			virtual void adapt(const vec &data){};
-			//! function redesigning the control strategy 
-			virtual void redesign(){};
-			//! returns designed control action
-			virtual vec ctrlaction(const vec &cond){return vec(0);}
-			
-			void from_setting(const Setting &set){
-				UI::get(rv,set,"rv",UI::optional);
-				UI::get(rvc,set,"rvc",UI::optional);
-			}
-			//! access function
-			const RV& _rv() {return rv;}
-			//! access function
-			const RV& _rvc() {return rvc;}
-			//! register this controller with given datasource under name "name"
-			virtual void log_register (logger &L, const string &prefix ) { }
-			//! write requested values into the logger
-			virtual void log_write ( ) const { }
-			
-	};
-	
+//! Base class for adaptive controllers
+//! The base class is, however, non-adaptive, method \c adapt() is empty.
+//! \note advanced Controllers will probably include estimator as their internal attribute (e.g. dual controllers)
+class Controller : public root {
+protected:
+	//! identifier of the designed action;
+	RV rv;
+	//! identifier of the conditioning variables - data needed ;
+	RV rvc;
+public:
+	//! function processing new observations and adapting control strategy accordingly
+	virtual void adapt ( const vec &data ) {};
+	//! function redesigning the control strategy
+	virtual void redesign() {};
+	//! returns designed control action
+	virtual vec ctrlaction ( const vec &cond ) const {return vec ( 0 );}
+
+	void from_setting ( const Setting &set ) {
+		UI::get ( rv,set,"rv",UI::optional );
+		UI::get ( rvc,set,"rvc",UI::optional );
+	}
+	//! access function
+	const RV& _rv() {return rv;}
+	//! access function
+	const RV& _rvc() {return rvc;}
+	//! register this controller with given datasource under name "name"
+	virtual void log_register ( logger &L, const string &prefix ) { }
+	//! write requested values into the logger
+	virtual void log_write ( ) const { }
+
+};
+
 //! Linear Quadratic Gaussian designer for constant penalizations and constant target
 //! Its internals are very close to Kalman estimator
 class LQG : public Controller {
-	protected:
-		//! StateSpace model from which we read data
-		shared_ptr<StateSpace<fsqmat> > S;
-		//! required value of the output y at time t (assumed constant)
-		vec y_req;
-		//! required value of the output y at time t (assumed constant)
-		vec u_req;
-		
-		//! Control horizon, set to maxint for infinite horizons
-		int horizon;
-		//! penalization matrix Qy
-		mat Qy;
-		//! penalization matrix Qu
-		mat Qu;
-		//! time of the design step - from horizon->0
-		int td;
-		//! controller parameters
-		mat L;
-		
-		//!@{ \name temporary storage for ricatti - use initialize
-		//! convenience parameters
-		int dimx;
-		//! convenience parameters
-		int dimy;
-		//! convenience parameters
-		int dimu;
-		
-		//!  parameters
-		mat pr;
-		//! penalization 
-		mat qux;
-		//! penalization
-		mat qyx;
-		//! internal quadratic form 
-		mat s;
-		//! penalization
-		mat qy;
-		//! pre_qr part
-		mat hqy;
-		//! pre qr matrix
-		mat pre_qr;
-		//! post qr matrix
-		mat post_qr;
-		//!@}
-		
-	public:
-		//! set system parameters from given matrices
-		void set_system(shared_ptr<StateSpace<fsqmat> > S0);
-		//! set penalization matrices and control horizon
-		void set_control_parameters(const mat &Qy0, const mat &Qu0, const vec &y_req0, int horizon0);
-		//! set system parameters from Kalman filter
-//		void set_system_parameters(const Kalman &K);
-		//! refresh temporary storage - inefficient can be improved
-		void initialize();
-		//! validation procedure
-		void validate();
-		//! function for future use which is called at each time td; Should call initialize()!
-		virtual void update_state(){};
-		//! redesign one step of the 
-		void ricatti_step(){
-			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");}
-			triu(post_qr);
-	// hn(m+1:2*m+n+r,m+1:2*m+n+r);
-			s=post_qr.get(dimu, 2*dimu+dimx+dimy-1, dimu, 2*dimu+dimx+dimy-1);
-		};
-		void redesign(){
-			for(td=horizon; td>0; td--){
-				update_state();
-				ricatti_step();
-			}
-/*			ws=hn(1:m,m+1:2*m+n+r);
-			wsd=hn(1:m,1:m);
-			Lklq=-inv(wsd)*ws;*/
-			L = -inv(post_qr.get(0,dimu-1, 0,dimu-1)) * post_qr.get(0,dimu-1, dimu, 2*dimu+dimx+dimy-1);
+protected:
+	//! StateSpace model from which we read data
+	shared_ptr<StateSpace<chmat> > S;
+	//! required value of the output y at time t (assumed constant)
+	vec y_req;
+	//! required value of the output y at time t (assumed constant)
+	vec u_req;
+
+	//! Control horizon, set to maxint for infinite horizons
+	int horizon;
+	//! penalization matrix Qy
+	chmat Qy;
+	//! penalization matrix Qu
+	chmat Qu;
+	//! time of the design step - from horizon->0
+	int td;
+	//! controller parameters
+	mat L;
+
+	//!@{ \name temporary storage for ricatti - use initialize
+	//! convenience parameters
+	int dimx;
+	//! convenience parameters
+	int dimy;
+	//! convenience parameters
+	int dimu;
+
+	//!  parameters
+	mat pr;
+	//! penalization
+	mat qux;
+	//! penalization
+	mat qyx;
+	//! internal quadratic form
+	mat s;
+	//! penalization
+	mat qy;
+	//! pre_qr part
+	mat hqy;
+	//! pre qr matrix
+	mat pre_qr;
+	//! post qr matrix
+	mat post_qr;
+	//!@}
+
+public:
+	//! set system parameters from given matrices
+	void set_system ( shared_ptr<StateSpace<chmat> > S0 );
+	//! update internal whan system has changed
+	void update_system() {
+		pr.set_submatrix ( 0,0,S->_B() );
+		pr.set_submatrix ( 0,dimu, S->_A() );
+
+		//penalization
+		qux.set_submatrix ( 0,0, Qu._Ch() );
+		qux.set_submatrix ( 0,dimx+dimu+dimy,Qu._Ch() );
+
+		qyx.set_submatrix ( 0,0,S->_C() );
+		qyx.set_submatrix ( 0,dimx,-eye ( dimy ) );
+
+		// parts of QR
+		hqy=Qy.to_mat()*qyx*pr;
+
+		// pre_qr
+		pre_qr = concat_vertical ( s*pr, concat_vertical ( hqy, qux ) );
+	}
+	//! set penalization matrices and control horizon
+	void set_control_parameters ( const mat &Qy0, const mat &Qu0, const vec &y_req0, int horizon0 );
+	//! set penalization matrices and control horizon
+	void set_control_Qy ( const mat &Qy0 ) {Qy=Qy0;}
+	//! refresh temporary storage - inefficient can be improved
+	void initial_belmann() {s=1e-5*eye ( dimx+dimu+dimy );};
+	//! validation procedure
+	void validate();
+	//! function for future use which is called at each time td; Should call initialize()!
+	//! redesign one step of the
+	void ricatti_step() {
+		pre_qr.set_submatrix ( 0,0,s*pr );
+		pre_qr.set_submatrix ( dimx+dimu+dimy, dimu+dimx, -Qy.to_mat()*y_req );
+		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);
+		s=post_qr.get ( dimu, 2*dimu+dimx+dimy-1, dimu, 2*dimu+dimx+dimy-1 );
+	};
+	void redesign() {
+		for ( td=horizon; td>0; td-- ) {
+			update_system();
+			ricatti_step();
 		}
-		//! compute control action
-		vec ctrlaction(const vec &state, const vec &ukm){vec pom=concat(state, ones(dimy), ukm); return L*pom;}
-	} ;
+		/*			ws=hn(1:m,m+1:2*m+n+r);
+					wsd=hn(1:m,1:m);
+					Lklq=-inv(wsd)*ws;*/
+		L = -inv ( post_qr.get ( 0,dimu-1, 0,dimu-1 ) ) * post_qr.get ( 0,dimu-1, dimu, 2*dimu+dimx+dimy-1 );
+	}
+	//! compute control action
+	vec ctrlaction ( const vec &state, const vec &ukm ) const {
+		vec pom=concat ( state, ones ( dimy ), ukm );
+		return L*pom;
+	}
+} ;
 
-			
+
 } // namespace
Index: library/bdm/design/arx2statespace_notes.lyx
===================================================================
--- library/bdm/design/arx2statespace_notes.lyx (revision 723)
+++ library/bdm/design/arx2statespace_notes.lyx (revision 723)
@@ -0,0 +1,151 @@
+#LyX 1.6.4 created this file. For more info see http://www.lyx.org/
+\lyxformat 345
+\begin_document
+\begin_header
+\textclass article
+\use_default_options true
+\language english
+\inputencoding auto
+\font_roman default
+\font_sans default
+\font_typewriter default
+\font_default_family default
+\font_sc false
+\font_osf false
+\font_sf_scale 100
+\font_tt_scale 100
+
+\graphics default
+\paperfontsize default
+\use_hyperref false
+\papersize default
+\use_geometry false
+\use_amsmath 1
+\use_esint 1
+\cite_engine basic
+\use_bibtopic false
+\paperorientation portrait
+\secnumdepth 3
+\tocdepth 3
+\paragraph_separation indent
+\defskip medskip
+\quotes_language english
+\papercolumns 1
+\papersides 1
+\paperpagestyle default
+\tracking_changes false
+\output_changes false
+\author "" 
+\end_header
+
+\begin_body
+
+\begin_layout Standard
+System
+\begin_inset Formula \[
+y_{t}=y_{t-1}+bu_{t}+e_{t}\]
+
+\end_inset
+
+ is assumed to fullfill causality, i.e.
+ 
+\begin_inset Formula $y_{t-1}\rightarrow u_{t}\rightarrow y_{t}$
+\end_inset
+
+ 
+\end_layout
+
+\begin_layout Standard
+with change in notation, i.e.
+ 
+\begin_inset Formula $y_{t-1}\rightarrow u_{t-1}\rightarrow y_{t}$
+\end_inset
+
+ the same system is: 
+\begin_inset Formula \[
+y_{t}=y_{t-1}+bu_{t-1}+e_{t}\]
+
+\end_inset
+
+and can be transformed as
+\begin_inset Formula \begin{eqnarray*}
+x_{t+1} & = & x_{t}+bu_{t}+e_{t}\\
+y_{t} & = & x_{t}\end{eqnarray*}
+
+\end_inset
+
+Morale:
+\end_layout
+
+\begin_layout Itemize
+for transcription from ARX to State-Space, the state is one time-step shorter
+ than the regressor in ARX!
+\end_layout
+
+\begin_layout Subsection*
+Delays in 
+\begin_inset Formula $u$
+\end_inset
+
+:
+\end_layout
+
+\begin_layout Standard
+System
+\begin_inset Formula \[
+y_{t}=y_{t-1}+bu_{t}+cu_{t-1}+e_{t}\]
+
+\end_inset
+
+is
+\begin_inset Formula \begin{eqnarray*}
+\left[\begin{array}{c}
+x_{t+1}\\
+u_{t}\end{array}\right] & =\left[\begin{array}{cc}
+1 & c\\
+0 & 0\end{array}\right] & \left[\begin{array}{c}
+x_{t}\\
+u_{t-1}\end{array}\right]+\left[\begin{array}{c}
+b\\
+1\end{array}\right]u_{t}+e_{t}\\
+y_{t} & = & x_{t}\end{eqnarray*}
+
+\end_inset
+
+Notation in BDM:
+\end_layout
+
+\begin_layout Standard
+xrv = [y{0}, u{-1}]
+\end_layout
+
+\begin_layout Subsection*
+Constant term
+\end_layout
+
+\begin_layout Standard
+System
+\begin_inset Formula \[
+y_{t}=y_{t-1}+bu_{t-1}+g+e_{t}\]
+
+\end_inset
+
+must be modelled as 
+\begin_inset Formula \[
+\left[\begin{array}{c}
+x_{t+1}\\
+one\end{array}\right]=\left[\begin{array}{cc}
+1 & c\\
+0 & 0\end{array}\right]\left[\begin{array}{c}
+x_{t}\\
+one\end{array}\right]+\left[\begin{array}{c}
+b\\
+1\end{array}\right]u_{t}+e_{t}\]
+
+\end_inset
+
+where one is an extra state variable.
+\end_layout
+
+\end_body
+\end_document
Index: library/bdm/base/bdmbase.h
===================================================================
--- library/bdm/base/bdmbase.h (revision 713)
+++ library/bdm/base/bdmbase.h (revision 723)
@@ -796,5 +796,5 @@
 	datalink_buffered() : datalink_part(), history ( 0 ), h2v_down ( 0 ), h2v_hist ( 0 ) {};
 	//! push current data to history
-	void step ( const vec &val_up ) {
+	void store_data ( const vec &val_up ) {
 		if ( v2h_up.length() > 0 ) {
 			history.shift_right ( 0, v2h_up.length() );
@@ -818,5 +818,5 @@
 	void set_connection ( const RV &rv, const RV &rv_up ) {
 		// create link between up and down
-		datalink_part::set_connection ( rv, rv_up );
+		datalink_part::set_connection ( rv, rv_up.remove_time() ); // only non-delayed version
 
 		// create rvs of history
@@ -871,6 +871,6 @@
 	//! update buffer
 	void step ( const vec &dt, const vec &ut ) {
-		dl1.step ( dt );
-		dl2.step ( ut );
+		dl1.store_data ( dt );
+		dl2.store_data ( ut );
 	}
 };
Index: library/bdm/base/datasources.h
===================================================================
--- library/bdm/base/datasources.h (revision 700)
+++ library/bdm/base/datasources.h (revision 723)
@@ -163,9 +163,9 @@
 	public:
 		void step() {
-			yt2rgr.step(yt); // y is now history
+			yt2rgr.store_data(yt); // y is now history
 			ut2rgr.filldown ( ut,rgr );
 			yt2rgr.filldown ( yt,rgr );
 			yt=ipdf->samplecond ( rgr );
-			ut2rgr.step(ut); //u is now history
+			ut2rgr.store_data(ut); //u is now history
 		}
 		void getdata ( vec &dt_out ) const {
@@ -175,5 +175,6 @@
 		}
 		void write(const vec &ut0){ut=ut0;}
-
+		void write(const vec &ut0, const ivec &ind){set_subvector(ut,ind,ut0);}
+		
 		/*!
 		\code
Index: library/bdm/stat/discrete.cpp
===================================================================
--- library/bdm/stat/discrete.cpp (revision 643)
+++ library/bdm/stat/discrete.cpp (revision 723)
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 713)
+++ library/bdm/stat/exp_family.h (revision 723)
@@ -708,5 +708,7 @@
 		//!access function
 		mat _R() const { return this->iepdf._R().to_mat(); }
-
+		//!access function
+		sq_T __R() const { return this->iepdf._R(); }
+		
 		//! Debug stream
 		template<typename sq_M>
Index: library/bdm/itpp_ext.cpp
===================================================================
--- library/bdm/itpp_ext.cpp (revision 712)
+++ library/bdm/itpp_ext.cpp (revision 723)
@@ -108,5 +108,5 @@
 bvec operator| (const bvec &a, const bvec &b)
 {
-	it_assert_debug (b.size() != a.size(), "operator&(): Vectors of different lengths");
+	it_assert_debug (b.size() == a.size(), "operator|(): Vectors of different lengths");
 
 	bvec temp (a.size());
Index: library/tests/testsuite/epdf_test.cpp
===================================================================
--- library/tests/testsuite/epdf_test.cpp (revision 722)
+++ library/tests/testsuite/epdf_test.cpp (revision 723)
Index: library/tests/testsuite/dirfile-format.matrix
===================================================================
--- library/tests/testsuite/dirfile-format.matrix (revision 723)
+++ library/tests/testsuite/dirfile-format.matrix (revision 723)
@@ -0,0 +1,4 @@
+r_0 RAW d 1
+r_1 RAW d 1
+th_alog RAW d 1
+th_blog RAW d 1
Index: library/tests/testsuite/datasource_test.cpp
===================================================================
--- library/tests/testsuite/datasource_test.cpp (revision 722)
+++ library/tests/testsuite/datasource_test.cpp (revision 723)
Index: library/tests/testsuite/logger_test.cpp
===================================================================
--- library/tests/testsuite/logger_test.cpp (revision 722)
+++ library/tests/testsuite/logger_test.cpp (revision 723)
@@ -64,7 +64,7 @@
 	remove_all ( ls.c_str() );
 	makedir ( ls, false );
-	remove_all ( "exp/dirfile" );
+	remove_all ( "dirfilelog_files" );
 
-	dirfilelog L ( "exp/dirfile", 10 );
+	dirfilelog L ( "dirfilelog_files", 10 );
 
 	int rid = L.add ( r, "" );
@@ -82,5 +82,5 @@
 
 	std::string expected ( load_test_file ( "logger_test_dirfile_format.matrix" ) );
-	std::string actual ( load_test_file ( "exp/dirfile/format" ) );
+	std::string actual ( load_test_file ( "dirfilelog_files/format" ) );
 	CHECK_EQUAL ( expected, actual );
 }
Index: library/tests/testsuite/LQG_test.cpp
===================================================================
--- library/tests/testsuite/LQG_test.cpp (revision 722)
+++ library/tests/testsuite/LQG_test.cpp (revision 723)
@@ -5,82 +5,86 @@
 using namespace bdm;
 
-TEST ( LQG_test ) {
+TEST(LQG_test) {
 	LQG reg;
-	shared_ptr<StateSpace<fsqmat> > stsp = new StateSpace<fsqmat>;
+	shared_ptr<StateSpace<chmat> > stsp=new StateSpace<chmat>;
 	// 2 x 1 x 1
-	stsp-> set_parameters ( eye ( 2 ), ones ( 2, 1 ), ones ( 1, 2 ), ones ( 1, 1 ), /* Q,R */ eye ( 2 ), eye ( 1 ) );
-	reg.set_system ( stsp ); // A, B, C
-	reg.set_control_parameters ( eye ( 1 ), eye ( 1 ),  vec_1 ( 1.0 ), 6 ); //Qy, Qu, horizon
+	stsp-> set_parameters(eye(2), ones(2,1), ones(1,2), ones(1,1), /* Q,R */ eye(2), eye(1));
+	reg.set_system(stsp); // A, B, C
+	reg.set_control_parameters(eye(1), eye(1),  vec_1(1.0), 6); //Qy, Qu, horizon
+	reg.validate();
 
 	reg.redesign();
-	double reg_apply = reg.ctrlaction ( "0.5, 1.1", "0.0" ) ( 0 ); /*convert vec to double*/
-	CHECK_CLOSE ( reg_apply, -0.248528137234392, 0.0001 );
+	double reg_apply=reg.ctrlaction("0.5, 1.1","0.0")(0); /*convert vec to double*/
+	CHECK_CLOSE(reg_apply, -0.248528137234392, 0.0001);
 }
 
-TEST ( to_state_test ) {
+TEST(to_state_test) {
 	mlnorm<fsqmat> ml;
-	mat A = "1.1, 2.3";
-	ml.set_parameters ( A, vec_1 ( 1.3 ), eye ( 1 ) );
-	RV yr = RV ( "y", 1 );
-	RV ur = RV ( "u", 1 );
-	ml.set_rv ( yr );
-	yr.t_plus ( -1 );
-	ml.set_rvc ( concat ( yr, ur ) );
-
-	shared_ptr<StateCanonical > Stsp = new StateCanonical;
-	Stsp->connect_mlnorm ( ml );
-
-	/* results from
+	mat A="1.1, 2.3";
+	ml.set_parameters(A, vec_1(1.3), eye(1));
+	RV yr=RV("y",1);
+	RV ur=RV("u",1);
+	ml.set_rv(yr);
+	yr.t_plus(-1);
+	ml.set_rvc(concat(yr, ur));
+	
+	shared_ptr<StateCanonical > Stsp=new StateCanonical;
+	Stsp->connect_mlnorm(ml);
+	
+	/* results from 
 	[A,B,C,D]=tf2ss([2.3 0],[1 -1.1])
 	*/
-	CHECK_CLOSE_EX ( Stsp->_A().get_row ( 0 ), vec ( "1.1" ), 0.0001 );
-	CHECK_CLOSE_EX ( Stsp->_C().get_row ( 0 ), vec ( "2.53" ), 0.0001 );
-	CHECK_CLOSE_EX ( Stsp->_D().get_row ( 0 ), vec ( "2.30" ), 0.0001 );
+	CHECK_CLOSE_EX(Stsp->_A().get_row(0), vec("1.1"), 0.0001);
+	CHECK_CLOSE_EX(Stsp->_C().get_row(0), vec("2.53"), 0.0001);
+	CHECK_CLOSE_EX(Stsp->_D().get_row(0), vec("2.30"), 0.0001);
 }
 
-TEST ( to_state_arx_test ) {
-	mlnorm<fsqmat> ml;
-	mat A = "1.1, 2.3";
-	ml.set_parameters ( A, vec_1 ( 1.3 ), eye ( 1 ) );
-	RV yr = RV ( "y", 1 );
-	RV ur = RV ( "u", 1 );
-	ml.set_rv ( yr );
-	yr.t_plus ( -1 );
-	ml.set_rvc ( concat ( yr, ur ) );
-
-	shared_ptr<StateFromARX> Stsp = new StateFromARX;
-	Stsp->connect_mlnorm ( ml );
-
-	/* results from
+TEST(to_state_arx_test) {
+	mlnorm<chmat> ml;
+	mat A="1.1, 2.3, 3.4";
+	ml.set_parameters(A, vec_1(1.3), eye(1));
+	RV yr=RV("y",1);
+	RV ur=RV("u",1);
+	ml.set_rv(yr);
+	ml.set_rvc(concat(yr.copy_t(-1), concat(ur, ur.copy_t(-1))));
+	
+	shared_ptr<StateFromARX> Stsp=new StateFromARX;
+	RV xrv; RV urv;
+	Stsp->connect_mlnorm(ml,xrv,urv);
+	
+	/* results from 
 	[A,B,C,D]=tf2ss([2.3 0],[1 -1.1])
 	*/
-	cout << "---" << endl;
+	CHECK_CLOSE_EX(Stsp->_A().get_row(0), vec("1.1, 3.4, 1.3"), 0.0001);
+	CHECK_CLOSE_EX(Stsp->_B().get_row(0), vec("2.3"), 0.0001);
+	CHECK_CLOSE_EX(Stsp->_C().get_row(0), vec("1, 0, 0"), 0.0001);
 }
 
-TEST ( arx_LQG_test ) {
-	mlnorm<fsqmat> ml;
-	mat A = "1.81, -.8189, .00468, .00438";
-	ml.set_parameters ( A, vec_1 ( 0.0 ), 0.00001*eye ( 1 ) );
-	RV yr = RV ( "y", 1 );
-	RV ur = RV ( "u", 1 );
-	RV rgr = yr.copy_t ( -1 );
-	rgr.add ( yr.copy_t ( -2 ) );
-	rgr.add ( yr.copy_t ( -2 ) );
-	rgr.add ( ur.copy_t ( -2 ) );
-	rgr.add ( ur.copy_t ( -1 ) );
-
-	ml.set_rv ( yr );
-	ml.set_rvc ( rgr );
+TEST(arx_LQG_test){
+	mlnorm<chmat> ml;
+	mat A="1.81, -.81, .00468, .00438";
+	ml.set_parameters(A, vec_1(0.0), 0.00001*eye(1));
+	RV yr=RV("y",1);
+	RV ur=RV("u",1);
+	RV rgr = yr.copy_t(-1);
+	rgr.add(yr.copy_t(-2));
+	rgr.add(yr.copy_t(-2));
+	rgr.add(ur.copy_t(-1));
+	rgr.add(ur);
+	
+	ml.set_rv(yr);
+	ml.set_rvc(rgr);
 	ml.validate();
-
-	shared_ptr<StateCanonical > Stsp = new StateCanonical;
-	Stsp->connect_mlnorm ( ml );
-
+	
+	shared_ptr<StateFromARX> Stsp=new StateFromARX;
+	RV xrv; RV urv;
+	Stsp->connect_mlnorm(ml,xrv,urv);
+	
 	LQG L;
-	L.set_system ( Stsp );
-	L.set_control_parameters ( eye ( 1 ), eye ( 1 ), vec_1 ( 0.0 ), 100 );
+	L.set_system(Stsp);
+	L.set_control_parameters(eye(1), sqrt(1.0/1000)*eye(1), vec_1(0.0), 100);
 	L.validate();
-
+	
 	L.redesign();
-	cout << L.to_string() << endl;
+	cout << L.to_string()<<endl;
 }
Index: library/tests/testsuite/datalink_test.cpp
===================================================================
--- library/tests/testsuite/datalink_test.cpp (revision 722)
+++ library/tests/testsuite/datalink_test.cpp (revision 723)
@@ -101,5 +101,5 @@
 
 	vec val_up ( "1 2" );
-	dl.step ( val_up );
+	dl.store_data(val_up);
 	val_up = "3 4";
 
