Index: /library/bdm/estim/arx.h
===================================================================
--- /library/bdm/estim/arx.h (revision 529)
+++ /library/bdm/estim/arx.h (revision 536)
@@ -110,7 +110,4 @@
 	//!\name Access attributes
 	//!@{
-	const egiw* _e() const {
-		return &est ;
-	};
 	const egiw& posterior() const {
 		return est;
Index: /library/bdm/estim/mixtures.h
===================================================================
--- /library/bdm/estim/mixtures.h (revision 477)
+++ /library/bdm/estim/mixtures.h (revision 536)
@@ -21,4 +21,5 @@
 namespace bdm {
 
+//! enum switch for internal approximation used in MixEF
 enum MixEF_METHOD { EM = 0, QB = 1};
 
@@ -118,7 +119,4 @@
 		return *est;
 	}
-	const eprod* _e() const {
-		return est;
-	}
 	emix* epredictor() const;
 	//! Flatten the density as if it was not estimated from the data
Index: /library/bdm/estim/kalman.h
===================================================================
--- /library/bdm/estim/kalman.h (revision 529)
+++ /library/bdm/estim/kalman.h (revision 536)
@@ -125,7 +125,4 @@
 		return est;
 	}
-	const enorm<sq_T>* _e() const {
-		return &est;
-	}
 	//!access function
 	mat& __K() {
@@ -214,7 +211,4 @@
 		return E;
 	};
-	const enorm<fsqmat>* _e() const {
-		return &E;
-	};
 	const mat _R() {
 		return P;
@@ -277,4 +271,5 @@
 	// TODO dodelat void to_setting( Setting &set ) const;
 
+	const enorm<chmat>& posterior() {return est;}
 };
 
@@ -423,5 +418,5 @@
 
 		est.set_rv ( RV ( "MM", A ( 0 )->posterior().dimension(), 0 ) );
-		est.set_parameters ( A ( 0 )->_e()->mean(), A ( 0 )->_e()->_R() );
+		est.set_parameters ( A ( 0 )->posterior().mean(), A ( 0 )->posterior()._R() );
 	}
 	void bayes ( const vec &dt ) {
@@ -439,6 +434,6 @@
 		case 1: {
 			int mi = max_index ( w );
-			const enorm<chmat>* st = ( Models ( mi )->_e() );
-			est.set_parameters ( st->mean(), st->_R() );
+			const enorm<chmat> &st = Models ( mi )->posterior() ;
+			est.set_parameters ( st.mean(), st._R() );
 		}
 		break;
@@ -451,9 +446,5 @@
 		}
 	}
-	//all posterior densities are equal => return the first one
-	const enorm<chmat>* _e() const {
-		return &est;
-	}
-	//all posterior densities are equal => return the first one
+	//! posterior density
 	const enorm<chmat>& posterior() const {
 		return est;
Index: /library/bdm/estim/particles.h
===================================================================
--- /library/bdm/estim/particles.h (revision 488)
+++ /library/bdm/estim/particles.h (revision 536)
@@ -111,5 +111,5 @@
 			dim = E.dimension() + A ( 0 )->posterior().dimension();
 			for ( int i = 0; i < _w.length() ; i++ ) {
-				Coms ( i ) = A ( i )->_e();
+				Coms ( i ) = &(A ( i )->posterior());
 			}
 		}
@@ -218,7 +218,4 @@
 		return jest;
 	}
-	const epdf* _e() const {
-		return &jest;    //Fixme: is it useful?
-	}
 	//! Set postrior of \c rvc to samples from epdf0. Statistics of BMs are not re-computed! Use only for initialization!
 	/*	void set_est ( const epdf& epdf0 ) {
@@ -234,5 +231,5 @@
 
 	//!Access function
-	BM* _BM ( int i ) {
+	const BM* _BM ( int i ) {
 		return BMs ( i );
 	}
Index: /library/bdm/shared_ptr.h
===================================================================
--- /library/bdm/shared_ptr.h (revision 529)
+++ /library/bdm/shared_ptr.h (revision 536)
@@ -171,4 +171,5 @@
 };
 
+//! Compare shared pointers 
 template<typename T, typename U>
 bool operator== ( shared_ptr<T> const &a, shared_ptr<U> const &b ) {
@@ -176,4 +177,5 @@
 }
 
+//! Compare shared pointers 
 template<typename T, typename U>
 bool operator!= ( shared_ptr<T> const &a, shared_ptr<U> const &b ) {
@@ -181,4 +183,5 @@
 }
 
+//! Compare shared pointers 
 template<typename T, typename U>
 bool operator< ( shared_ptr<T> const &a, shared_ptr<U> const &b ) {
Index: /library/bdm/math/square_mat.cpp
===================================================================
--- /library/bdm/math/square_mat.cpp (revision 495)
+++ /library/bdm/math/square_mat.cpp (revision 536)
@@ -64,13 +64,10 @@
 
 ldmat::ldmat ( const mat &V ) : sqmat ( V.cols() ) {
-//TODO check if correct!! Based on heuristic observation of lu()
 
 	it_assert_debug ( dim == V.rows(), "ldmat::ldmat matrix V is not square!" );
 
 	// L and D will be allocated by ldform()
-
 	//Chol is unstable
 	this->ldform ( chol ( V ), ones ( dim ) );
-//	this->ldform(ul(V),ones(dim));
 }
 
Index: /library/bdm/mex/config2mxstruct.h
===================================================================
--- /library/bdm/mex/config2mxstruct.h (revision 477)
+++ /library/bdm/mex/config2mxstruct.h (revision 536)
@@ -8,7 +8,10 @@
 using namespace libconfig;
 
+//! Reimplementation of libconfig's Config class for Matlab mxArray structures
 class UImxConfig : public Config {
 public:
+	//! Matlab structure where the info is stored
 	mxArray *mxconfig;
+	//! Load file in libconfig syntax to Matlab arrays
 	UImxConfig ( const char * filename ) {
 		Config config;
@@ -16,4 +19,5 @@
 		mxconfig = group2mxstruct ( config.getRoot() );
 	}
+	//! Convert existing Setting to Matlab arrays
 	UImxConfig ( const Setting &setting ) {
 		mxconfig = group2mxstruct ( setting );
@@ -21,5 +25,5 @@
 
 private:
-
+	//! Convert libconfig's array to Matlab vector
 	mxArray* array2mxvector ( const Setting &setting )  {
 		if ( !setting.isArray() ) mexErrMsgTxt ( "Given setting is not an array" );
@@ -40,4 +44,5 @@
 	}
 
+	//! Convert libconfig's array to Matlab matrix
 	mxArray* list2mxmatrix ( const Setting &setting )  {
 		if ( !setting.isList() || ( "matrix" != setting[0] ) )
@@ -58,5 +63,6 @@
 		return result;
 	}
-
+	
+	//! Convert libconfig's gourp to Matlab structure
 	mxArray* group2mxstruct ( const Setting &setting ) {
 		if ( !setting.isGroup() ) mexErrMsgTxt ( "Given setting is not a group." );
@@ -102,5 +108,5 @@
 
 	}
-
+	//! Convert libconfig's list  to Matlab cell
 	mxArray* list2mxcell ( const Setting &setting )  {
 		if ( !setting.isList() ) mexErrMsgTxt ( "Given setting is not a list." );
Index: /library/bdm/mex/mex_parser.h
===================================================================
--- /library/bdm/mex/mex_parser.h (revision 477)
+++ /library/bdm/mex/mex_parser.h (revision 536)
@@ -8,6 +8,8 @@
 using namespace libconfig;
 
+//! Class for writing Settings into Matlab's mxArray
 class UImxArray : public Config {
 public:
+	//! Build an instance of Config with fields filled from the given \a mxarray
 	UImxArray ( const mxArray *mxarray ) : Config() {
 		Setting & setting = this->getRoot(); //setting is a group
@@ -22,4 +24,5 @@
 		setAutoConvert ( true );
 	}
+	//! Add libconfig's \c list to the structure
 	void addList ( const mxArray *mxarray, const char* name ) {
 		Setting & setting = this->getRoot(); //setting is a group
@@ -27,4 +30,5 @@
 		fillList ( child, mxarray );
 	}
+	//! Add libconfig's \c group to the structure
 	void addGroup ( const mxArray *mxarray, const char* name ) {
 		Setting & setting = this->getRoot(); //setting is a group
@@ -32,4 +36,5 @@
 		fillGroup ( child, mxarray );
 	}
+	//! Operator for more convenient access to this Confic
 	operator Setting&() {
 		return getRoot();
Index: /library/bdm/mex/mex_logger.h
===================================================================
--- /library/bdm/mex/mex_logger.h (revision 529)
+++ /library/bdm/mex/mex_logger.h (revision 536)
@@ -37,4 +37,6 @@
 
 }
+
+//! Logger storing results into an mxArray
 class mexlog : public memlog {
 public:
Index: /library/bdm/base/bdmbase.h
===================================================================
--- /library/bdm/base/bdmbase.h (revision 532)
+++ /library/bdm/base/bdmbase.h (revision 536)
@@ -29,5 +29,7 @@
 
 typedef std::map<string, int> RVmap;
+//! Internal global variable storing sizes of RVs
 extern ivec RV_SIZES;
+//! Internal global variable storing names of RVs
 extern Array<string> RV_NAMES;
 
@@ -160,6 +162,4 @@
 		}
 		//!@}
-
-		//TODO why not inline and later??
 
 		//! \name Algebra on Random Variables
@@ -192,5 +192,5 @@
 		//!@{
 
-		//! generate \c str from rv, by expanding sizes TODO to_string..
+		//! generate \c str from rv, by expanding sizes 
 		str tostr() const;
 		//! when this rv is a part of bigger rv, this function returns indices of self in the data vector of the bigger crv.
@@ -223,5 +223,4 @@
 
 		// TODO dodelat void to_setting( Setting &set ) const;
-
 		//! Invalidate all named RVs. Use before initializing any RV instances, with care...
 		static void clear_all();
@@ -391,6 +390,5 @@
 SHAREDPTR(epdf);
 
-//! Conditional probability density, e.g. modeling some dependencies.
-//TODO Samplecond can be generalized
+//! Conditional probability density, e.g. modeling \f$ f( x | y) \f$, where \f$ x \f$ is random variable, \c rv, and \f$ y \f$ is conditioning variable, \c rvc.
 class mpdf : public root
 {
@@ -405,7 +403,9 @@
 
 	protected:
+		//! set internal pointer \c ep to point to given \c iepdf
 		void set_ep (epdf &iepdf) {
 			ep = &iepdf;
 		}
+		//! set internal pointer \c ep to point to given \c iepdf
 		void set_ep (epdf *iepdfp) {
 			ep = iepdfp;
@@ -487,8 +487,10 @@
 SHAREDPTR(mpdf);
 
+//! Mpdf with internal epdf that is modified by function \c condition
 template <class EPDF>
 class mpdf_internal: public mpdf
 {
 	protected :
+		//! Internal epdf used for sampling
 		EPDF iepdf;
 	public:
@@ -555,4 +557,5 @@
 		//! Constructor
 		datalink() : downsize (0), upsize (0) { }
+		//! Conevnience constructor
 		datalink (const RV &rv, const RV &rv_up) {
 			set_connection (rv, rv_up);
@@ -607,5 +610,5 @@
 		//! Constructor
 		datalink_m2e() : condsize (0) { }
-
+		//! Set connection between vectors
 		void set_connection (const RV &rv, const RV &rvc, const RV &rv_up) {
 			datalink::set_connection (rv, rv_up);
@@ -621,5 +624,5 @@
 			return tmp;
 		}
-
+		//! Copy corresponding values to Up.condition
 		void pushup_cond (vec &val_up, const vec &val, const vec &cond) {
 			it_assert_debug (downsize == val.length(), "Wrong val");
@@ -643,4 +646,5 @@
 		//! Constructor
 		datalink_m2m() {};
+		//! Set connection between the vectors
 		void set_connection (const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up) {
 			datalink_m2e::set_connection (rv, rvc, rv_up);
@@ -713,10 +717,10 @@
 class mepdf : public mpdf
 {
-
+		//! Internal shared pointer to epdf
 		shared_ptr<epdf> iepdf;
 	public:
 		//!Default constructor
 		mepdf() { }
-
+		//! Set internal shared pointer
 		mepdf (shared_ptr<epdf> em) {
 			iepdf = em;
@@ -893,5 +897,5 @@
 			return NULL;
 		};
-		//!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t})
+		//!Constructs conditional density of 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t}) \f$
 		virtual mpdf* predictor() const {
 			it_error ("Not implemented");
@@ -939,5 +943,4 @@
 		}
 		virtual const epdf& posterior() const = 0;
-		virtual const epdf* _e() const = 0;
 		//!@}
 
Index: /library/bdm/stat/merger.cpp
===================================================================
--- /library/bdm/stat/merger.cpp (revision 507)
+++ /library/bdm/stat/merger.cpp (revision 536)
@@ -122,5 +122,5 @@
 			if ( DBG ) {
 				cout << "Resampling =" << 1. / sum_sqr ( w ) << endl;
-				cout << Mix._e()->mean() << endl;
+				cout << Mix.posterior().mean() << endl;
 				cout << sum ( Smp_ex, 2 ) / Npoints << endl;
 				cout << Smp_ex*Smp_ex.T() / Npoints << endl;
Index: /library/bdm/stat/exp_family.h
===================================================================
--- /library/bdm/stat/exp_family.h (revision 535)
+++ /library/bdm/stat/exp_family.h (revision 536)
@@ -44,6 +44,4 @@
 		//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
 		virtual double lognc() const = 0;
-		//!TODO decide if it is really needed
-		virtual void dupdate (mat &v) {it_error ("Not implemented");};
 		//!Evaluate normalized log-probability
 		virtual double evallog_nn (const vec &val) const{it_error ("Not implemented");return 0.0;};
@@ -220,4 +218,5 @@
 		ldmat est_theta_cov() const;
 
+		//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
 		void mean_mat (mat &M, mat&R) const;
 		//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
@@ -340,6 +339,7 @@
 			if (evalll) {last_lognc = est.lognc();}
 		}
-		const epdf& posterior() const {return est;};
-		const eDirich* _e() const {return &est;};
+		//! reimplemnetation of BM::posterior()
+		const eDirich& posterior() const {return est;};
+		//! constructor function		
 		void set_parameters (const vec &beta0) {
 			est.set_parameters (beta0);
@@ -373,10 +373,9 @@
 
 		vec sample() const;
-		//! TODO: is it used anywhere?
-//	mat sample ( int N ) const;
 		double evallog (const vec &val) const;
 		double lognc () const;
-		//! Returns poiter to alpha and beta. Potentially dengerous: use with care!
+		//! Returns pointer to internal alpha. Potentially dengerous: use with care!
 		vec& _alpha() {return alpha;}
+		//! Returns pointer to internal beta. Potentially dengerous: use with care!
 		vec& _beta() {return beta;}
 		vec mean() const {return elem_div (alpha, beta);}
@@ -482,5 +481,4 @@
 		//!@}
 
-		double eval (const vec &val) const  {return nk;}
 		double evallog (const vec &val) const  {
 			if (any (val < low) && any (val > high)) {return inf;}
@@ -516,5 +514,5 @@
  \brief Normal distributed linear function with linear function of mean value;
 
- Mean value \f$mu=A*rvc+mu_0\f$.
+ Mean value \f$ \mu=A*\mbox{rvc}+\mu_0 \f$.
 */
 template < class sq_T, template <typename> class TEpdf = enorm >
@@ -524,4 +522,5 @@
 		//! Internal epdf that arise by conditioning on \c rvc
 		mat A;
+		//! Constant additive term 
 		vec mu_const;
 //			vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
@@ -558,4 +557,5 @@
 		mat _R() { return this->iepdf._R().to_mat(); }
 
+		//! Debug stream
 		template<typename sq_M>
 		friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);
@@ -649,10 +649,14 @@
 {
 	protected:
+		//! Variable \f$ \Lambda \f$ from theory
 		ldmat Lambda;
+		//! Reference to variable \f$ R \f$
 		ldmat &_R;
+		//! Variable \f$ R_e \f$
 		ldmat Re;
 	public:
 		mlstudent () : mlnorm<ldmat, enorm> (),
 				Lambda (),	_R (iepdf._R()) {}
+		//! constructor function
 		void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
 			it_assert_debug (A0.rows() == mu0.length(), "");
@@ -968,5 +972,7 @@
 		double delta;
 	public:
+		//! Set internal structures
 		void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
+		//! Sample matrix argument
 		mat sample_mat() const {
 			mat X = zeros (p, p);
@@ -998,11 +1004,18 @@
 };
 
+//! Inverse Wishart on Choleski decomposition
+/*! Being computed by conversion from `standard' Wishart
+*/
 class eiWishartCh: public epdf
 {
 	protected:
+		//! Internal instance of Wishart density
 		eWishartCh W;
+		//! size of Ch
 		int p;
+		//! parameter delta
 		double delta;
 	public:
+		//! constructor function
 		void set_parameters (const mat &Y0, const double delta0) {
 			delta = delta0;
@@ -1011,4 +1024,5 @@
 		}
 		vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
+		//! access function
 		void _setY (const vec &y0) {
 			mat Ch (p, p);
@@ -1044,4 +1058,5 @@
 };
 
+//! Random Walk on inverse Wishart
 class rwiWishartCh : public mpdf_internal<eiWishartCh>
 {
@@ -1049,12 +1064,14 @@
 		//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
 		double sqd;
-		//reference point for diagonal
+		//!reference point for diagonal
 		vec refl;
+		//! power of the reference
 		double l;
+		//! dimension
 		int p;
 
 	public:
 		rwiWishartCh() : sqd (0), l (0), p (0) {}
-
+		//! constructor function
 		void set_parameters (int p0, double k, vec ref0, double l0) {
 			p = p0;
@@ -1326,4 +1343,5 @@
 void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};
 
+//! \todo unify this stuff with to_string()
 template<class sq_T>
 std::ostream &operator<< (std::ostream &os,  mlnorm<sq_T> &ml)
Index: /library/bdm/stat/merger.h
===================================================================
--- /library/bdm/stat/merger.h (revision 529)
+++ /library/bdm/stat/merger.h (revision 536)
@@ -318,4 +318,5 @@
 SHAREDPTR ( merger_base );
 
+//! Merger using importance sampling with mixture proposal density
 class merger_mix : public merger_base {
 protected:
Index: /library/bdm/stat/emix.h
===================================================================
--- /library/bdm/stat/emix.h (revision 529)
+++ /library/bdm/stat/emix.h (revision 536)
@@ -165,4 +165,5 @@
 
 	shared_ptr<epdf> marginal ( const RV &rv ) const;
+	//! Update already existing marginal density  \c target
 	void marginal ( const RV &rv, emix &target ) const;
 	shared_ptr<mpdf> condition ( const RV &rv ) const;
@@ -290,5 +291,5 @@
 		set_elements ( mFacs );
 	}
-
+	//! Set internal \c mpdfs from given values
 	void set_elements (const Array<shared_ptr<mpdf> > &mFacs );
 
@@ -339,11 +340,4 @@
 		return smp;
 	}
-	mat samplecond ( const vec &cond,  int N ) {
-		mat Smp ( dimension(), N );
-		for ( int i = 0; i < N; i++ ) {
-			Smp.set_col ( i, samplecond ( cond ) );
-		}
-		return Smp;
-	}
 
 	//! Load from structure with elements:
@@ -372,5 +366,7 @@
 	Array<datalink*> dls;
 public:
+	//! Default constructor
 	eprod () : epdfs ( 0 ), dls ( 0 ) {};
+	//! Set internal 
 	void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
 		epdfs = epdfs0;//.set_length ( epdfs0.length() );
@@ -470,5 +466,5 @@
 	//! Set weights \c w and components \c R
 	void set_parameters ( const vec &w0, const Array<shared_ptr<mpdf> > &Coms0 ) {
-		//!\TODO check if all components are OK
+		//!\todo check if all components are OK
 		Coms = Coms0;
 		w=w0;	
Index: /library/tests/epdf_harness.cpp
===================================================================
--- /library/tests/epdf_harness.cpp (revision 529)
+++ /library/tests/epdf_harness.cpp (revision 536)
@@ -210,5 +210,5 @@
 	Array<vec> actual(CurrentContext::max_trial_count);
 	do {
-		mat smp = mep.samplecond ( vec ( 0 ), nsamples );
+		mat smp = mep.samplecond_m ( vec ( 0 ), nsamples );
 		vec emu = sum ( smp, 2 ) / nsamples;
 		actual( tc ) = emu;
@@ -231,5 +231,5 @@
 	Array<mat> actual(CurrentContext::max_trial_count);
 	do {
-		mat smp = mep.samplecond ( vec ( 0 ), nsamples );
+		mat smp = mep.samplecond_m ( vec ( 0 ), nsamples );
 		vec emu = sum ( smp, 2 ) / nsamples;
 		mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
Index: /library/tests/arx_elem_test.cpp
===================================================================
--- /library/tests/arx_elem_test.cpp (revision 477)
+++ /library/tests/arx_elem_test.cpp (revision 536)
@@ -12,5 +12,5 @@
 	mat mu ( 1, 1 );
 	mat R ( 1, 1 );
-	Ar._e()->mean_mat ( mu, R );
+	Ar.posterior().mean_mat ( mu, R );
 	cout << "Prior moments: mu=" << mu << ", R=" << R << endl;
 
@@ -24,5 +24,5 @@
 	// Ar is now filled with estimates of N(0,1);
 	cout << "Empirical moments: mu=" << sum ( smp ) / ndat << ", R=" << sum_sqr ( smp ) / ndat - pow ( sum ( smp ) / ndat, 2 ) << endl;
-	Ar._e()->mean_mat ( mu, R );
+	Ar.posterior().mean_mat ( mu, R );
 	cout << "Posterior moments: mu=" << mu << ", R=" << R << endl;
 
Index: /library/doc/tutorial/tut_arx.dox
===================================================================
--- /library/doc/tutorial/tut_arx.dox (revision 290)
+++ /library/doc/tutorial/tut_arx.dox (revision 536)
@@ -36,8 +36,10 @@
 \begin{array}{c} [y_{t}',\,\psi_{t}']\\ \\\end{array}
 +(1-\phi) V_0
-\f]</dd>
-<dt>"Degree of freedom"</dd>  <dd>which is an accumulator of number of data records \f[
+\f]
+</dd>
+<dt>"Degree of freedom"</dt>  <dd>which is an accumulator of number of data records \f[
 \nu_t = \phi \nu_{t-1} + 1 + (1-\phi) \nu_0
-\f]</dd>
+\f]
+</dd>
 </dl>
 where \f$ \phi \f$ is the forgetting factor, typically \f$ \phi \in [0,1]\f$ roughly corresponding to the effective length of the exponential window by relation:\f[
