/*!
  \file
  \brief Basic structures of probability calculus: random variables, probability densities, Bayes rule
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef BDMBASE_H
#define BDMBASE_H

#include <map>

#include "../itpp_ext.h"
#include "../bdmroot.h"
#include "../shared_ptr.h"
#include "user_info.h"

using namespace libconfig;
using namespace itpp;
using namespace std;

namespace bdm {
//! Structure of RV, i.e. RVs expanded into a flat list of IDs, used for debugging.
class str {
public:
	//! vector id ids (non-unique!)
	ivec ids;
	//! vector of times
	ivec times;
	//!Default constructor
	str ( ivec ids0, ivec times0 ) : ids ( ids0 ), times ( times0 ) {
		bdm_assert ( times0.length() == ids0.length(), "Incompatible input" );
	};
};

/*!
* \brief Class representing variables, most often random variables

The purpose of this class is to decribe a vector of data. Such description is used for connecting various vectors between each other, see class datalink.

The class is implemented using global variables to assure uniqueness of description:

 In is a vector
\dot
digraph datalink {
rankdir=LR;
subgraph cluster0 {
node [shape=record];
label = "MAP \n std::map<string,int>";
map [label="{{\"a\"| \"b\" | \"c\"} | {<3> 3 |<1> 1|<2> 2}}"];
color = "white"
}
subgraph cluster1{
node [shape=record];
label = "NAMES";
names [label="{<1> \"b\" | <2> \"c\" | <3>\"a\" }"];
color = "white"
}
subgraph cluster2{
node [shape=record];
label = "SIZES";
labelloc = b;
sizes [label="{<1>1 |<2> 4 |<3> 1}"];
color = "white"
}
map:1 -> names:1;
map:1 -> sizes:1;
map:3 -> names:3;
map:3 -> sizes:3;
}
\enddot
*/

class RV : public root {
private:
	typedef std::map<string, int> str2int_map;

	//! Internal global variable storing sizes of RVs
	static ivec SIZES;
	//! Internal global variable storing names of RVs
	static Array<string> NAMES;

	//! TODO
	const static int BUFFER_STEP;

	//! TODO
	static str2int_map MAP;

public:

protected:
	//! size of the data vector
	int dsize;
	//! number of individual rvs
	int len;
	//! Vector of unique IDs
	ivec ids;
	//! Vector of shifts from current time
	ivec times;

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
	//!@{

	//! Full constructor
	RV ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) {
		init ( in_names, in_sizes, in_times );
	}

	//! Constructor with times=0
	RV ( const Array<std::string> &in_names, const ivec &in_sizes ) {
		init ( in_names, in_sizes, zeros_i ( in_names.length() ) );
	}

	//! Constructor with sizes=1, times=0
	RV ( const Array<std::string> &in_names ) {
		init ( in_names, ones_i ( in_names.length() ), zeros_i ( in_names.length() ) );
	}

	//! Constructor of empty RV
	RV() : dsize ( 0 ), len ( 0 ), ids ( 0 ), times ( 0 ) {}

	//! Constructor of a single RV with given id
	RV ( string name, int sz, int tm = 0 );

	// compiler-generated copy constructor is used
	//!@}

	//! \name Access functions
	//!@{

	//! State output, e.g. for debugging.
	friend std::ostream &operator<< ( std::ostream &os, const RV &rv );

	string to_string() const {ostringstream o; o << *this; return o.str();}
	
	//! total size of a random variable
	int _dsize() const {
		return dsize;
	}

	//! access function
	const ivec& _ids() const {
		return ids;
	}

	//! Recount size of the corresponding data vector
	int countsize() const;
	//! Vector of cumulative sizes of RV
	ivec cumsizes() const;
	//! Number of named parts
	int length() const {
		return len;
	}
	int id ( int at ) const {
		return ids ( at );
	}
	int size ( int at ) const {
		return SIZES ( ids ( at ) );
	}
	int time ( int at ) const {
		return times ( at );
	}
	std::string name ( int at ) const {
		return NAMES ( ids ( at ) );
	}
	//! returns name of a scalar at position scalat, i.e. it can be in the middle of vector name, in that case it adds "_%d" to it
	std::string scalarname ( int scalat ) const {
		bdm_assert(scalat<dsize,"Wrong input index");
		int id=0;
		int scalid=0;
		while (scalid+SIZES(ids(id))<=scalat)  { scalid+=SIZES(ids(id)); id++;};
		//now id is the id of variable of interest
		if (size(id)==1)
			return  NAMES ( ids ( id ) );
		else
			return  NAMES ( ids ( id ) )+ "_" + num2str(scalat-scalid);

	}
	void set_time ( int at, int time0 ) {
		times ( at ) = time0;
	}
	//!@}

	//! \name Algebra on Random Variables
	//!@{

	//! 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;
	//! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
	bool add ( const RV &rv2 );
	//! Subtract  another variable from the current one
	RV subt ( const RV &rv2 ) const;
	//! Select only variables at indices ind
	RV subselect ( const ivec &ind ) const;

	//! Select only variables at indices ind
	RV operator() ( const ivec &ind ) const {
		return subselect ( ind );
	}

	//! Select from data vector starting at di1 to di2
	RV operator() ( int di1, int di2 ) const;

	//! Shift \c time by delta.
	void t_plus ( int delta );
	//!@}

	//! @{ \name Time manipulation functions
	//! returns rvs with time set to 0 and removed duplicates
	RV remove_time() const {
		return RV ( unique ( ids ), dummy );
	}
	//! create new RV from the current one with time shifted by given value
	RV copy_t(int dt) const {
		RV tmp=*this;
		tmp.t_plus(dt);
		return tmp;
	}
	//! return rvs with expanded delayes and sorted in the order of: \f$ [ rv_{0}, rv_{-1},\ldots  rv_{max_delay}]\f$
	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_plus ( -1 );
			tmp.add ( rvt ); //shift u1
		}
		return tmp;
	}
	//!@}

	//!\name Relation to vectors
	//!@{

	//! 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.
	//! 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)
	void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
	//! Minimum time-offset
	int mint() const {
			return times.length()>0 ? min (times) : 0;
	}
	//!@}

	/*! \brief UI for class RV (description of data vectors)

	\code
	class = 'RV';                   
	names = {'a', 'b', 'c', ...};   // UNIQUE IDENTIFIER same names = same variable
									// names are also used when storing results 
	--- optional ---
	sizes = [1, 2, 3, ...];         // size of each name. default = ones()
									// if size = -1, it is found out from previous instances of the same name
	times = [-1, -2, 0, ...];       // time shifts with respect to current time, default = zeros()
	\endcode
	*/
	void from_setting ( const Setting &set );

	// TODO dodelat void to_setting( Setting &set ) const;

	//! Invalidate all named RVs. Use before initializing any RV instances, with care...
	static void clear_all();
	//! function for debugging RV related stuff
	string show_all();

};
UIREGISTER ( RV );
SHAREDPTR ( RV );

//! Concat two random variables
RV concat ( const RV &rv1, const RV &rv2 );

/*!
@brief Class for storing results (and semi-results) of an experiment

This class abstracts logging of results from implementation. This class replaces direct logging of results (e.g. to files or to global variables) by calling methods of a logger. Specializations of this abstract class for specific storage method are designed.
*/
class logger : public root {
	protected:
		//! RVs of all logged variables.
		Array<RV> entries;
		//! Names of logged quantities, e.g. names of algorithm variants
		Array<string> names;
		//!separator of prefixes of entries
		const string separator;
	public:
		//!Default constructor
		logger(const string separator0) : entries ( 0 ), names ( 0 ), separator(separator0) {}
		
		//! returns an identifier which will be later needed for calling the \c logit() function
		//! For empty RV it returns -1, this entry will be ignored by \c logit().
		virtual int add ( const RV &rv, string prefix = "" ) {
			int id;
			if ( rv._dsize() > 0 ) {
				id = entries.length();
				names = concat ( names, prefix ); // diff
				entries.set_length ( id + 1, true );
				entries ( id ) = rv;
			} else {
				id = -1;
			}
			return id; // identifier of the last entry
		}
		
		//! log this vector
		virtual void logit ( int id, const vec &v ) {
			bdm_error ( "Not implemented" );
		};
		//! log this double
		virtual void logit ( int id, const double &d ) {
			bdm_error ( "Not implemented" );
		};
		
		//! Shifts storage position for another time step.
		virtual void step() {
			bdm_error ( "Not implemneted" );
		};
		
		//! Finalize storing information
		virtual void finalize() {};
		
		//! Initialize the storage
		virtual void init() {};
		
		//!separator of prefixes for this logger
		const string& prefix_sep() {return separator;}
};


//! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv
class fnc : public root {
protected:
	//! Length of the output vector
	int dimy;
	//! Length of the input vector
	int dimc;
public:
	//!default constructor
	fnc() {};
	//! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
	virtual vec eval ( const vec &cond ) {
		return vec ( 0 );
	};

	//! function substitutes given value into an appropriate position
	virtual void condition ( const vec &val ) {};

	//! access function
	int dimension() const {
		return dimy;
	}
	//! access function
	int dimensionc() const {
		return dimc;
	}
};

class epdf;

//! 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 {
protected:
	//!dimension of the condition
	int dimc;
	//! random variable in condition
	RV rvc;

	//! TODO
	int dim;

	//! TODO
	RV rv;

public:
	//! \name Constructors
	//! @{

	mpdf() : dimc ( 0 ), rvc(), dim(0), rv() { }

	mpdf ( const mpdf &m ) : dimc ( m.dimc ), rvc ( m.rvc ), rv( m.rv ), dim( m.dim) { }
	//!@}

	//! \name Matematical operations
	//!@{

	//! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
	virtual vec samplecond ( const vec &cond ) {
		bdm_error ( "Not implemented" );
		return vec();
	}

	//! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
	virtual mat samplecond_m ( const vec &cond, int N );

	//! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
	virtual double evallogcond ( const vec &dt, const vec &cond ) {
		bdm_error ( "Not implemented" );
		return 0.0;
	}

	//! Matrix version of evallogcond
	virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {
		vec v ( Dt.cols() );
		for ( int i = 0; i < Dt.cols(); i++ ) {
			v ( i ) = evallogcond ( Dt.get_col ( i ), cond );
		}
		return v;
	}

	//! Array<vec> version of evallogcond
	virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
		bdm_error ( "Not implemented" );
		return vec();
	}

	//! \name Access to attributes
	//! @{

	const RV& _rv() const {
		return rv;
	}
	const RV& _rvc() const {
		return rvc;
	}

	int dimension() const {
		return dim;
	}
	int dimensionc() {
		return dimc;
	}

	//! Load from structure with elements:
	//!  \code
	//! { class = "mpdf_offspring",
	//!   rv = {class="RV", names=(...),}; // RV describing meaning of random variable
	//!   rvc= {class="RV", names=(...),}; // RV describing meaning of random variable in condition
	//!   // elements of offsprings
	//! }
	//! \endcode
	//!@}
	void from_setting ( const Setting &set );
	//!@}

	//! \name Connection to other objects
	//!@{
	void set_rvc ( const RV &rvc0 ) {
		rvc = rvc0;
	}
	void set_rv ( const RV &rv0 ) {
		rv = rv0;
	}

	bool isnamed() {
		return ( dim == rv._dsize() ) && ( dimc == rvc._dsize() );
	}
	//!@}
};
SHAREDPTR ( mpdf );

//! Probability density function with numerical statistics, e.g. posterior density.
class epdf : public mpdf {
protected:
	//! dimension of the random variable
	int dim;
	//! Description of the random variable
	RV rv;

public:
	/*! \name Constructors
	 Construction of each epdf should support two types of constructors:
	\li empty constructor,
	\li copy constructor,

	The following constructors should be supported for convenience:
	\li constructor followed by calling \c set_parameters()
	\li constructor accepting random variables calling \c set_rv()

	 All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors.
	@{*/
	epdf() : dim ( 0 ), rv() {};
	epdf ( const epdf &e ) : dim ( e.dim ), rv ( e.rv ) {};
	epdf ( const RV &rv0 ) : dim ( rv0._dsize() ) {
		set_rv ( rv0 );
	};
	void set_parameters ( int dim0 ) {
		dim = dim0;
	}
	//!@}

	//! \name Matematical Operations
	//!@{

	//! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
	virtual vec sample() const {
		bdm_error ( "not implemented" );
		return vec();
	}

	//! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$
	virtual mat sample_m ( int N ) const;

	//! Compute log-probability of argument \c val
	//! In case the argument is out of suport return -Infinity
	virtual double evallog ( const vec &val ) const {
		bdm_error ( "not implemented" );
		return 0.0;
	}

	//! Compute log-probability of multiple values argument \c val
	virtual vec evallog_m ( const mat &Val ) const;

	//! Compute log-probability of multiple values argument \c val
	virtual vec evallog_m ( const Array<vec> &Avec ) const;

	//! Return conditional density on the given RV, the remaining rvs will be in conditioning
	virtual shared_ptr<mpdf> condition ( const RV &rv ) const;

	//! Return marginal density on the given RV, the remainig rvs are intergrated out
	virtual shared_ptr<epdf> marginal ( const RV &rv ) const;

	//! return expected value
	virtual vec mean() const {
		bdm_error ( "not implemneted" );
		return vec();
	}

	//! return expected variance (not covariance!)
	virtual vec variance() const {
		bdm_error ( "not implemneted" );
		return vec();
	}

	//! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
	virtual void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {
		vec mea = mean();
		vec std = sqrt ( variance() );
		lb = mea - 2 * std;
		ub = mea + 2 * std;
	};
	//!@}

	//! \name Connection to other classes
	//! Description of the random quantity via attribute \c rv is optional.
	//! For operations such as sampling \c rv does not need to be set. However, for \c marginalization
	//! and \c conditioning \c rv has to be set. NB:
	//! @{

	//!Name its rv
	void set_rv ( const RV &rv0 ) {
		rv = rv0;
	}

	//! True if rv is assigned
	bool isnamed() const {
		bool b = ( dim == rv._dsize() );
		return b;
	}

	//! Return name (fails when isnamed is false)
	const RV& _rv() const {
		//bdm_assert_debug ( isnamed(), "" );
		return rv;
	}
	//! store values of the epdf on the following levels:
	//!  #1 mean
	//!  #2 mean + lower & upper bound
	void log_register(logger &L, const string &prefix){
		RV r;
		if ( isnamed() ) {
			r = _rv();
		} else {
			r = RV ( "", dimension() );
		};
		root::log_register(L,prefix);
		logrec->ids.set_size(3);
		if (log_level >0){
			logrec->ids(0) = logrec->L.add ( r, prefix + logrec->L.prefix_sep()+ "mean" );
		}
		if (log_level >1){
			logrec->ids(1) = logrec->L.add ( r, prefix + logrec->L.prefix_sep()+ "lb" );
			logrec->ids(2) = logrec->L.add ( r, prefix + logrec->L.prefix_sep()+ "ub" );
		}
	}
	//!@}

	//! \name Access to attributes
	//! @{

	//! Size of the random variable
	int dimension() const {
		return dim;
	}
	//! Load from structure with elements:
	//!  \code
	//! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable
	//!   // elements of offsprings
	//! }
	//! \endcode
	//!@}
	void from_setting ( const Setting &set ) {
		shared_ptr<RV> r = UI::build<RV> ( set, "rv", UI::optional );
		if ( r ) {
			set_rv ( *r );
		}
	}


/*! \brief Unconditional mpdf, allows using epdf in the role of mpdf.

*/


/// MEPDF BEGINS HERE TODO
	//! empty
	vec samplecond ( const vec &cond ) {
		return sample();
	}
	double evallogcond ( const vec &val, const vec &cond ) {
		return evallog ( val );
	}
};
SHAREDPTR ( epdf );

//! 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:
	//! constructor
	mpdf_internal() : mpdf(), iepdf() {
//		set_ep ( iepdf ); TODO!
	}

	//! Update \c iepdf so that it represents this mpdf conditioned on \c rvc = cond
	//! This function provides convenient reimplementation in offsprings
	virtual void condition ( const vec &cond ) {
		bdm_error ( "Not implemented" );
	}

	//!access function to iepdf
	EPDF& e() {
		return iepdf;
	}

	//! Reimplements samplecond using \c condition()
	vec samplecond ( const vec &cond );
	//! Reimplements evallogcond using \c condition()
	double evallogcond ( const vec &val, const vec &cond );
	//! Efficient version of evallogcond for matrices
	virtual vec evallogcond_m ( const mat &Dt, const vec &cond );
	//! Efficient version of evallogcond for Array<vec>
	virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond );
	//! Efficient version of samplecond
	virtual mat samplecond_m ( const vec &cond, int N );
};

/*! \brief DataLink is a connection between two data vectors Up and Down

Up can be longer than Down. Down must be fully present in Up (TODO optional)
See chart:
\dot
digraph datalink {
  node [shape=record];
  subgraph cluster0 {
    label = "Up";
      up [label="<1>|<2>|<3>|<4>|<5>"];
    color = "white"
}
  subgraph cluster1{
    label = "Down";
    labelloc = b;
      down [label="<1>|<2>|<3>"];
    color = "white"
}
    up:1 -> down:1;
    up:3 -> down:2;
    up:5 -> down:3;
}
\enddot

*/
class datalink {
protected:
	//! Remember how long val should be
	int downsize;

	//! Remember how long val of "Up" should be
	int upsize;

	//! val-to-val link, indices of the upper val
	ivec v2v_up;

public:
	//! Constructor
	datalink() : downsize ( 0 ), upsize ( 0 ) { }

	//! Convenience constructor
	datalink ( const RV &rv, const RV &rv_up ) {
		set_connection ( rv, rv_up );
	}

	//! set connection, rv must be fully present in rv_up
	virtual void set_connection ( const RV &rv, const RV &rv_up );

	//! set connection using indices
	virtual void set_connection ( int ds, int us, const ivec &upind );

	//! 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"
	virtual void filldown ( const vec &val_up, vec &val_down ) {
		bdm_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
		val_down = val_up ( v2v_up );
	}
	//! Fill val of "Up" by my pieces
	virtual void pushup ( vec &val_up, const vec &val ) {
		bdm_assert_debug ( downsize == val.length(), "Wrong val" );
		bdm_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
		set_subvector ( val_up, v2v_up, val );
	}
	//! access functions
	int _upsize() {
		return upsize;
	}
	//! access functions
	int _downsize() {
		return downsize;
	}
	//! for future use
	virtual ~datalink(){}
};

/*! 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;
	//! v2h: indeces of up too be pushed to h
	ivec v2h_up;
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 ( v2h_up.length() > 0 ) {
			history.shift_right ( 0, v2h_up.length() );
		  history.set_subvector ( 0, val_up(v2h_up) ); 
		}
	}
	//! 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 rv_hist0 =rv_hist.remove_time();  // these RVs will form history at time =0
		// now we need to know what is needed from Up
		rv_hist = rv_hist.expand_delayes(); // full regressor - including time 0
		Hrv=rv_hist.subt(rv_hist0);   // remove time 0
		history = zeros ( Hrv._dsize() );

		// decide if we need to copy val to history
		if (Hrv._dsize() >0) {
			v2h_up = rv_hist0.dataind(rv_up); // indeces of elements of rv_up to be copied
		} // else v2h_up is empty
		
		Hrv.dataind ( rv, h2v_hist, h2v_down );

		downsize = v2v_down.length() + h2v_down.length();
		upsize = v2v_up.length();
	}
	//! set history of variable given by \c rv1 to values of \c hist.
	void set_history(const RV& rv1, const vec &hist0){
		bdm_assert(rv1._dsize()==hist0.length(),"hist is not compatible with given rv1");
		ivec ind_H;
		ivec ind_h0;
		Hrv.dataind(rv1, ind_H, ind_h0); // find indeces of rv in 
		set_subvector(history, ind_H, hist0(ind_h0)); // copy given hist to appropriate places
	}
};

//! buffered datalink from 2 vectors to 1
class datalink_2to1_buffered {
protected:
	//! link 1st vector to down
	datalink_buffered dl1;
	//! link 2nd vector to down
	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 );
	}
	//! fill values of down from the values of the two up vectors
	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 );
	}
	//! update buffer
	void step ( const vec &dt, const vec &ut ) {
		dl1.step ( dt );
		dl2.step ( ut );
	}
};



//! Data link with a condition.
class datalink_m2e: public datalink {
protected:
	//! Remember how long cond should be
	int condsize;

	//!upper_val-to-local_cond link, indices of the upper val
	ivec v2c_up;

	//!upper_val-to-local_cond link, indices of the local cond
	ivec v2c_lo;

public:
	//! Constructor
	datalink_m2e() : condsize ( 0 ) { }

	//! Set connection between vectors
	void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up );

	//!Construct condition
	vec get_cond ( const vec &val_up );

	//! Copy corresponding values to Up.condition
	void pushup_cond ( vec &val_up, const vec &val, const vec &cond );
};

//!DataLink is a connection between mpdf and its superordinate (Up)
//! This class links
class datalink_m2m: public datalink_m2e {
protected:
	//!cond-to-cond link, indices of the upper cond
	ivec c2c_up;
	//!cond-to-cond link, indices of the local cond
	ivec c2c_lo;

public:
	//! 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 );
		//establish c2c connection
		rvc.dataind ( rvc_up, c2c_lo, c2c_up );
		bdm_assert_debug ( c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given" );
	}

	//! Get cond for myself from val and cond of "Up"
	vec get_cond ( const vec &val_up, const vec &cond_up ) {
		vec tmp ( condsize );
		set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) );
		set_subvector ( tmp, c2c_lo, cond_up ( c2c_up ) );
		return tmp;
	}
	//! Fill

};


//! \brief Combines RVs from a list of mpdfs to a single one.
RV get_composite_rv ( const Array<shared_ptr<mpdf> > &mpdfs, bool checkoverlap = false );

/*! \brief Abstract class for discrete-time sources of data.

The class abstracts operations of: 
\li  data aquisition, 
\li  data-preprocessing, such as  scaling of data, 
\li  data resampling from the task of estimation and control.
Moreover, for controlled systems, it is able to receive the desired control action and perform it in the next step. (Or as soon as possible).

The DataSource has three main data interaction structures:
\li input, \f$ u_t \f$, 
\li output \f$ y_t \f$,
\li data, \f$ d_t=[y_t,u_t, \ldots ]\f$ a collection of all inputs and outputs and possibly some internal variables too.

*/

class DS : public root {
protected:
		//! size of data returned by \c getdata()
	int dtsize;
		//! size of data
	int utsize;
		//!size of output
		int ytsize;
	//!Description of data returned by \c getdata().
	RV Drv;
	//!Description of data witten by by \c write().
	RV Urv; //
		//!Description of output data
		RV Yrv; //
public:
	//! default constructors
		DS() : Drv(), Urv(),Yrv() {log_level=1;};

	//! Returns maximum number of provided data, by default it is set to maximum allowed length, shorter DS should overload this method! See, MemDS.max_length().
	virtual int max_length() {return std::numeric_limits< int >::max();}
	//! Returns full vector of observed data=[output, input]
	virtual void getdata ( vec &dt ) const {
		bdm_error ( "abstract class" );
	}
	//! Returns data records at indeces.
	virtual void getdata ( vec &dt, const ivec &indeces ) {
		bdm_error ( "abstract class" );
	}

	//! Accepts action variable and schedule it for application.
	virtual void write (const vec &ut ) {
		bdm_error ( "abstract class" );
	}

	//! Accepts action variables at specific indeces
	virtual void write (const vec &ut, const ivec &indeces ) {
		bdm_error ( "abstract class" );
	}

	//! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system.
	virtual void step() 
	{
		bdm_error ( "abstract class" );
	}

	//! Register DS for logging into logger L
	virtual void log_register (logger &L,  const string &prefix ) {
		bdm_assert ( ytsize == Yrv._dsize(), "invalid DS: ytsize (" + num2str ( ytsize ) + ") different from Drv " + num2str ( Yrv._dsize() ) );
		bdm_assert ( utsize == Urv._dsize(), "invalid DS: utsize (" + num2str ( utsize ) + ") different from Urv " + num2str ( Urv._dsize() ) );

		root::log_register(L,prefix);
		//we know that 
		if (log_level >0){
			logrec->ids.set_size(2);
			logrec->ids(0) = logrec->L.add ( Yrv, "" );
			logrec->ids(1) = logrec->L.add ( Urv, "" );
		}
	}
	//! Register DS for logging into logger L
	virtual void log_write ( ) const {
		if (log_level >0) {
			vec tmp ( Yrv._dsize() + Urv._dsize());
			getdata ( tmp );
			// d is first in getdata
			logrec->L.logit ( logrec->ids(0), tmp.left ( Yrv._dsize() ) );
			// u follows after d in getdata
			logrec->L.logit ( logrec->ids(1), tmp.mid ( Yrv._dsize(), Urv._dsize() ) );
		}
	}
	//!access function
	virtual const RV& _drv() const {
		return Drv;
	}
	//!access function
	const RV& _urv() const {
		return Urv;
	}
	//!access function
	const RV& _yrv() const {
		return Yrv;
	}
	//! set random variables
	virtual void set_drv (const  RV &yrv, const RV &urv) {
		Yrv = yrv;
		Drv = concat(yrv,urv);
		Urv = urv;
	}
};

/*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.

This object represents exact or approximate evaluation of the Bayes rule:
\f[
f(\theta_t | d_1,\ldots,d_t) = \frac{f(y_t|\theta_t,\cdot) f(\theta_t|d_1,\ldots,d_{t-1})}{f(y_t|d_1,\ldots,d_{t-1})}
\f]

Access to the resulting posterior density is via function \c posterior().

As a "side-effect" it also evaluates log-likelihood of the data, which can be accessed via function _ll().
It can also evaluate predictors of future values of \f$y_t\f$, see functions epredictor() and predictor().

Alternatively, it can evaluate posterior density conditioned by a known constant, \f$ c_t \f$:
\f[
f(\theta_t | c_t, d_1,\ldots,d_t) \propto  f(y_t,\theta_t|c_t,\cdot, d_1,\ldots,d_{t-1})
\f]

The value of \f$ c_t \f$ is set by function condition().

*/

class BM : public root {
protected:
	//! Random variable of the data (optional)
	RV drv;
	//!Logarithm of marginalized data likelihood.
	double ll;
	//!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
	bool evalll;
	
public:
	//! \name Constructors
	//! @{

	BM() : ll ( 0 ), evalll ( true ) { };
	BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
	//! \brief Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
	//! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
	virtual BM* _copy_() const {
		return NULL;
	};
	//!@}

	//! \name Mathematical operations
	//!@{

	/*! \brief Incremental Bayes rule
	@param dt vector of input data
	*/
	virtual void bayes ( const vec &dt ) = 0;
	//! Batch Bayes rule (columns of Dt are observations)
	virtual void bayesB ( const mat &Dt );
	//! Evaluates predictive log-likelihood of the given data record
	//! I.e. marginal likelihood of the data with the posterior integrated out.
	virtual double logpred ( const vec &dt ) const {
		bdm_error ( "Not implemented" );
		return 0.0;
	}

	//! Matrix version of logpred
	vec logpred_m ( const mat &dt ) const {
		vec tmp ( dt.cols() );
		for ( int i = 0; i < dt.cols(); i++ ) {
			tmp ( i ) = logpred ( dt.get_col ( i ) );
		}
		return tmp;
	}

	//!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$
	virtual epdf* epredictor() const {
		bdm_error ( "Not implemented" );
		return NULL;
	};
	//!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 {
		bdm_error ( "Not implemented" );
		return NULL;
	};
	//!@}

	//! \name Extension to conditional BM
	//! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
	//! Alternatively, it can be used for automated connection to DS when the condition is observed
	//!@{

	//! Name of extension variable
	RV rvc;
	//! access function
	const RV& _rvc() const {
		return rvc;
	}

	//! Substitute \c val for \c rvc.
	virtual void condition ( const vec &val ) {
		bdm_error ( "Not implemented!" );
	}

	//!@}


	//! \name Access to attributes
	//!@{
	//! access function
	const RV& _drv() const {
		return drv;
	}
	//! access function
	void set_drv ( const RV &rv ) {
		drv = rv;
	}
	//! access to rv of the posterior
	void set_rv ( const RV &rv ) {
		const_cast<epdf&> ( posterior() ).set_rv ( rv );
	}
	//! return internal log-likelihood of the last data vector
	double _ll() const {
		return ll;
	}
	//! switch evaluation of log-likelihood on/off
	void set_evalll ( bool evl0 ) {
		evalll = evl0;
	}
	//! return posterior density
	virtual const epdf& posterior() const = 0;
	//!@}

	//! \name Logging of results
	//!@{

	//! Set boolean options from a string, recognized are: "logbounds,logll"
	virtual void set_options ( const string &opt ) {
		if ( opt.find ( "logbounds" ) != string::npos ) {
			const_cast<epdf&>(posterior()).set_log_level(2) ;
		} else {
			const_cast<epdf&>(posterior()).set_log_level(1) ;
		}
		if ( opt.find ( "logll" ) != string::npos ) {
			log_level = 1;
		}
	}

	//! Add all logged variables to a logger
	//! Log levels two digits: xy where
	//!  * y = 0/1 log-likelihood is to be logged
	//!  * x = level of the posterior (typically 0/1/2 for nothing/mean/bounds)
	virtual void log_register ( logger &L, const string &prefix = "" ) {
		root::log_register(L,prefix);		

		const_cast<epdf&>(posterior()).log_register(L, prefix+L.prefix_sep()+"apost"); 
		
		if ((log_level) > 0){
			logrec->ids.set_size(1);
			logrec->ids(0) = L.add(RV("ll",1), prefix+L.prefix_sep()+"ll");
		}
	}
	//! Save results to the given logger, details of what is stored is configured by \c LIDs and \c options
	virtual void log_write ( ) const {
		posterior().log_write();
		if (log_level >0) { logrec->L.logit ( logrec->ids ( 0 ), ll );}
	}
	//!@}
	void from_setting ( const Setting &set ) {
		shared_ptr<RV> r = UI::build<RV> ( set, "drv", UI::optional );
		if ( r ) {
			set_drv ( *r );
		}
		string opt;
		if ( UI::get ( opt, set, "options", UI::optional ) ) {
			set_options ( opt );
		}
	}

};
//! array of pointers to epdf
typedef Array<shared_ptr<epdf> > epdf_array;
//! array of pointers to mpdf
typedef Array<shared_ptr<mpdf> > mpdf_array;

template<class EPDF>
vec mpdf_internal<EPDF>::samplecond ( const vec &cond ) {
	condition ( cond );
	vec temp = iepdf.sample();
	return temp;
}

template<class EPDF>
mat mpdf_internal<EPDF>::samplecond_m ( const vec &cond, int N ) {
	condition ( cond );
	mat temp ( dimension(), N );
	vec smp ( dimension() );
	for ( int i = 0; i < N; i++ ) {
		smp = iepdf.sample();
		temp.set_col ( i, smp );
	}

	return temp;
}

template<class EPDF>
double mpdf_internal<EPDF>::evallogcond ( const vec &dt, const vec &cond ) {
	double tmp;
	condition ( cond );
	tmp = iepdf.evallog ( dt );
	return tmp;
}

template<class EPDF>
vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Dt, const vec &cond ) {
	condition ( cond );
	return iepdf.evallog_m ( Dt );
}

template<class EPDF>
vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
	condition ( cond );
	return iepdf.evallog_m ( Dt );
}

}; //namespace
#endif // BDMBASE_H