Index: /bdm/estim/libKF.h
===================================================================
--- /bdm/estim/libKF.h (revision 83)
+++ /bdm/estim/libKF.h (revision 85)
@@ -262,5 +262,5 @@
 		A ( dimx,dimx ), B ( dimx,dimu ), C ( dimy,dimx ), D ( dimy,dimu ),
 		Q(dimx), R(dimy),
-		est ( rv ), fy ( rvy ), _mu(est._mu()), _P(est._R()), _yp(fy._mu()),_Ry(fy._R()) {
+		est ( rv ), fy ( rvy ), _yp(fy._mu()),_Ry(fy._R()), _mu(est._mu()), _P(est._R()) {
 
 	this->set_parameters ( K0.A, K0.B, K0.C, K0.D, K0.R, K0.Q );
@@ -279,5 +279,5 @@
 		A ( dimx,dimx ), B ( dimx,dimu ), C ( dimy,dimx ), D ( dimy,dimu ),
 		Q(dimx), R (dimy),
-		est ( rvx ), fy ( rvy ), _mu(est._mu()), _P(est._R()), _yp(fy._mu()),_Ry(fy._R()) {
+		est ( rvx ), fy ( rvy ),  _yp(fy._mu()),_Ry(fy._R()),_mu(est._mu()), _P(est._R()) {
 };
 
Index: /bdm/math/libDC.h
===================================================================
--- /bdm/math/libDC.h (revision 75)
+++ /bdm/math/libDC.h (revision 85)
@@ -149,4 +149,7 @@
 		vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;};
 
+		//! Add another matrix in fsq form with weight w
+		void add ( const fsqmat &fsq2, double w=1.0 ){M+=fsq2.M;};
+
 		//! Access functions
 		void setD (const vec &nD){M=diag(nD);}
@@ -155,4 +158,5 @@
 		//! Access functions
 		void setD (const vec &nD, int i){for(int j=i;j<nD.length();j++){M(j,j)=nD(j-i);}} //Fixme can be more general
+
 
 		//! add another fsqmat matrix
Index: /bdm/stat/libFN.h
===================================================================
--- /bdm/stat/libFN.h (revision 62)
+++ /bdm/stat/libFN.h (revision 85)
@@ -18,5 +18,5 @@
 using namespace itpp;
 
-//! class representing function $f(x) = a$, here rv is empty
+//! class representing function \f$f(x) = a\f$, here \c rv is empty
 class constfn : public fnc
 {
@@ -32,8 +32,8 @@
 };
 
-//! Class representing function $f(x) = Ax+B$
+//! Class representing function \f$f(x) = Ax+B\f$
 class linfn: public fnc
 {
-		//! Identification of $x$
+		//! Identification of \f$x\f$
 		RV rv;
 		//! Matrix A
@@ -53,5 +53,5 @@
 
 /*!
-\brief Class representing a differentiable function of two variables $f(x,u)$.
+\brief Class representing a differentiable function of two variables \f$f(x,u)\f$.
 
 Function of two variables.
@@ -73,5 +73,5 @@
 		int dimu;
 	public:
-		//! Evaluates $f(x0,u0)$ (VS: Do we really need common eval? )
+		//! Evaluates \f$f(x0,u0)\f$ (VS: Do we really need common eval? )
 		vec eval ( const vec &cond )
 		{
@@ -80,9 +80,9 @@
 		};
 
-		//! Evaluates $f(x0,u0)$
+		//! Evaluates \f$f(x0,u0)\f$
 		virtual vec eval ( const vec &x0, const vec &u0 ) {return zeros ( dimy );};
-		//! Evaluates \f$A=\frac{d}{dx}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed. @param x0 numeric value of $x$, @param u0 numeric value of $u$ @param A a place where the result will be stored.
+		//! Evaluates \f$A=\frac{d}{dx}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed. @param x0 numeric value of \f$x\f$, @param u0 numeric value of \f$u\f$ @param A a place where the result will be stored.
 		virtual void dfdx_cond ( const vec &x0, const vec &u0, mat &A , bool full=true ) {};
-		//! Evaluates \f$A=\frac{d}{du}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed.	@param x0 numeric value of $x$, @param u0 numeric value of $u$ @param A a place where the result will be stored.
+		//! Evaluates \f$A=\frac{d}{du}f(x,u)|_{x0,u0}\f$ and writes result into \c A . @param full denotes that even unchanged entries are to be rewritten. When, false only the changed elements are computed.	@param x0 numeric value of \f$x\f$, @param u0 numeric value of \f$u\f$ @param A a place where the result will be stored.
 		virtual void dfdu_cond ( const vec &x0, const vec &u0, mat &A, bool full=true ) {};
 		//!Default constructor (dimy is not set!)
@@ -94,5 +94,5 @@
 };
 
-//! Class representing function $f(x,u) = Ax+Bu$
+//! Class representing function \f$f(x,u) = Ax+Bu\f$
 //TODO can be generalized into multilinear form!
 class bilinfn: public diffbifn
Index: /bdm/stat/libDS.h
===================================================================
--- /bdm/stat/libDS.h (revision 33)
+++ /bdm/stat/libDS.h (revision 85)
@@ -22,5 +22,5 @@
 * \brief Class representing off-line data stored in memory
 
-The data are stored in an internal matrix \c Data . Each column of Data corresponds to one discrete time observation $t$. Access to this matrix is via indexes \c rowid and \c delays.
+The data are stored in an internal matrix \c Data . Each column of Data corresponds to one discrete time observation \f$t\f$. Access to this matrix is via indexes \c rowid and \c delays.
 
 The data can be loaded from a file.
Index: /bdm/stat/libBM.cpp
===================================================================
--- /bdm/stat/libBM.cpp (revision 32)
+++ /bdm/stat/libBM.cpp (revision 85)
@@ -24,6 +24,6 @@
 	sizes = in_sizes;
 	times = in_times;
-	size = 0;
-	for(i=0;i<len;i++){size+=sizes(i);}
+	tsize = 0;
+	for(i=0;i<len;i++){tsize+=sizes(i);}
 };
 
@@ -32,9 +32,9 @@
 }
 
-RV::RV () : size(0),len(0){};
+RV::RV () : tsize(0),len(0){};
 
 void RV::add (const RV &rv2) {
 	// TODO
-	size+=rv2.size;
+	tsize+=rv2.tsize;
 	len +=rv2.len;
 	ids=concat(ids,rv2.ids);
@@ -78,5 +78,5 @@
 
 ivec RV::indexlist(){
-	ivec indlist(size);
+	ivec indlist(tsize);
 	int i;
 	int pos = 0;
Index: /bdm/stat/libEF.h
===================================================================
--- /bdm/stat/libEF.h (revision 77)
+++ /bdm/stat/libEF.h (revision 85)
@@ -191,5 +191,5 @@
  \brief Normal distributed linear function with linear function of mean value;
 
- Mean value $mu=A*rvc$.
+ Mean value \f$mu=A*rvc\f$.
 */
 template<class sq_T>
@@ -197,6 +197,6 @@
 	//! Internal epdf that arise by conditioning on \c rvc
 	enorm<sq_T> epdf;
+	mat A;
 	vec& _mu; //cached epdf.mu;
-	mat A;
 public:
 	//! Constructor
@@ -216,5 +216,5 @@
 
 Mean value, \f$\mu\f$, of this density is given by \c rvc .
-Standard deviation of the random walk is proportional to one $k$-th the mean.
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
 This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
 
@@ -225,5 +225,5 @@
 	//! Internal epdf that arise by conditioning on \c rvc
 	egamma epdf;
-	//! Constant $k$
+	//! Constant \f$k\f$
 	double k;
 	//! cache of epdf.beta
@@ -245,8 +245,8 @@
 \brief  Gamma random walk around a fixed point
 
-Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, $p$. $k$ is the coefficient of the geometric combimation
+Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
 \f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
 
-Standard deviation of the random walk is proportional to one $k$-th the mean.
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
 This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
 
@@ -280,5 +280,5 @@
 	//! Number of particles
 	int n;
-	//! Sample weights $w$
+	//! Sample weights \f$w\f$
 	vec w;
 	//! Samples \f$x^{(i)}, i=1..n\f$
Index: /bdm/stat/libBM.h
===================================================================
--- /bdm/stat/libBM.h (revision 62)
+++ /bdm/stat/libBM.h (revision 85)
@@ -28,5 +28,5 @@
 protected:
 	//! size = sum of sizes
-	int size;
+	int tsize;
 	//! len = number of individual rvs
 	int len;
@@ -54,6 +54,8 @@
 	friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
 
-	//! Return length (number of scalars) of the RV.
-	int count() const {return size;} ;
+	//! Return number of scalars in the RV.
+	int count() const {return tsize;} ;
+	//! Return length (number of entries) of the RV.
+	int length() const {return len;} ;
 
 	//TODO why not inline and later??
@@ -62,5 +64,5 @@
 	ivec find ( RV rv2 );
 	//! Add (concat) another variable to the current one
-	void add (const RV &rv2 );
+ 	void add (const RV &rv2 );
 	//! Add (concat) another variable to the current one
 	friend RV concat (const RV &rv1, const RV &rv2 );
@@ -78,8 +80,17 @@
 	//!access function
 	Array<std::string>& _names(){return names;};
-};
-
-
-//! Class representing function $f(x)$ of variable $x$ represented by \c rv
+
+	//!access function
+	int id(int at){return ids(at);};
+	//!access function
+	int size(int at){return sizes(at);};
+	//!access function
+	int time(int at){return times(at);};
+	//!access function
+	std::string name(int at){return names(at);};
+};
+
+
+//! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv
 
 class fnc {
@@ -90,5 +101,5 @@
 	//!default constructor
 	fnc(int dy):dimy(dy){};
-	//! function evaluates numerical value of $f(x)$ at $x=cond$
+	//! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
 	virtual vec eval ( const vec &cond ) {
 		return vec ( 0 );
@@ -199,5 +210,5 @@
 	void linkrvs ( RV &drv, RV &urv );
 
-	//! Moves from $t$ to $t+1$, i.e. perfroms the actions and reads response of the system.
+	//! Moves from \f$t\f$ to \f$t+1\f$, i.e. perfroms the actions and reads response of the system.
 	void step();
 
@@ -242,5 +253,5 @@
 \brief Conditional Bayesian Filter
 
-Evaluates conditional filtering density $f(rv|rvc,data)$ for a given \c rvc which is specified in each step by calling function \c condition.
+Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition.
 
 This is an interface class used to assure that certain BM has operation \c condition .
