/*!
  \file
  \brief Probability distributions for discrete support densities
  \author Vaclav Smidl.

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

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

#ifndef DISCR_H
#define DISCR_H


#include "../shared_ptr.h"
#include "../base/bdmbase.h"
#include "../math/chmat.h"

namespace bdm {

	class support_base: public root {
	protected:
		//! dimension
		int dim;
		//! Number of data points
		int Npoints;
		//! active vector for first_vec and next_vec
		vec actvec;
		//! volume of the active point
		double actvolume;
	public:
		//! default constructor
		support_base() : dim ( 0 ), Npoints ( 0 ) {
		}
						
		//! set the first vector to corner and store result in actvec
		virtual const vec& first_vec() NOT_IMPLEMENTED(empty_vec);
		
		//! Get next active vector, call ONLY after first_vector()!
		virtual const vec& next_vec() NOT_IMPLEMENTED(empty_vec);;
		
		//! return active vector, call ONLY after first_vector()!
		virtual const vec& act_vec() {
			return actvec;
		};
		//! return active vector, call ONLY after first_vector()!
		virtual const double& act_volume() {
			return actvolume;
		};
		
		//! Access function
		int points() const {
			return Npoints;
		}
				
	};
	UIREGISTER ( support_base );
	
//! Rectangular support
//! Support points are located inbetween ranges! For example:
//! For ranges=[0,1] and gridsizes=[1] the support point is 0.5
class rectangular_support: public support_base {
protected:
    //! Array of boundaries (2D vectors: [begining,end]) for each dimension
    Array<vec> ranges;
    //! Number of support points in each dimension
    ivec gridsizes;
	//! indices of active vector
	vec actvec_ind;
	//
	vec steps;
public:
    //! default constructor
    rectangular_support() : support_base() {
    }

    //! set parameters
    void set_parameters ( const Array<vec> &ranges0, const ivec &gridsize0 );

    //! Internal functio to set temporaries correctly
    void initialize();

    //! return vector at position given by vector of indices
    vec get_vec ( const ivec &inds );

    //! convert dimension indices into linear index, the indexing is in the same way as in \c next_vec()
    long linear_index ( const ivec inds );

    //! set the first vector to corner and store result in actvec
    const vec& first_vec();

    //! Get next active vector, call ONLY after first_vector()!
    const vec& next_vec();

    //! return active vector, call ONLY after first_vector()!
    const vec& act_vec() {
        return actvec;
    };

    //! \todo to je asi navic .. v predkovi!
    ivec nearest_point ( const vec &val );

    //! Access function
    int points() const {
        return Npoints;
    }

    //! access function
    //! \todo opet pouze do potomka..
    const vec& _steps() const {
        return steps;
    }

    void from_setting ( const Setting &set );
};
UIREGISTER ( rectangular_support );

//! Discrete support with stored support points
class discrete_support: public support_base {
protected:
    //! storage of support points
    Array<vec> Spoints;
    //! index in iterators
    int idx;
	//! volumes for each point
	vec volumes;
public:
    //! Default constructor
    discrete_support() : Spoints ( 0 ), idx ( 0 ), volumes(0) {}
    //! Access function
    int points() const {
        return Spoints.length();
    }
    //! set the first vector to corner and store result in actvec
    const vec& first_vec() {

		bdm_assert_debug ( Spoints.length() > 0, "Empty support" );
        idx = 0;
		actvolume=volumes(0);
        return Spoints ( idx );
    }
    //! set next vector after calling first_vec()
    const vec& next_vec() {
        bdm_assert_debug ( Spoints.length() > idx - 1, "Out of support points" );
		idx++;
		actvolume=volumes(idx);
		actvec = Spoints ( idx );
        return actvec;
    }

    /*!
    \code
      class = "discrete_support";
      points = ( [1,2..], [2,2..], ...); // list of points
       === OR ===
      epdf = {class="epdf_offspring",...}; // epdf rfom which to sample
      npoints = 100;                     // number of samples
    \endcode
    */
    void from_setting ( const Setting &set );

	void validate(){
		bdm_assert(volumes.length()==Spoints.length(),"dimension mismatch");
		Npoints=Spoints.length();
		if (Npoints>0)
			dim = Spoints(0).length();
	}
    //! access function
    Array<vec> & _Spoints() {
        return Spoints;
    }
};
UIREGISTER ( discrete_support );

//! Function defined by values on a fixed grid and interpolated inbetween them
class grid_fnc: public fnc {
protected:
    //! grid - function support
    shared_ptr<support_base> sup;
    //! function values on the grid
    vec values;
public:
    //! constructor function
    void set_support ( shared_ptr<support_base> &sup0 ) {
        sup = sup0;
        values = zeros ( sup->points() );
    }
    //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
    void set_values ( double ( *evalptr ) ( const vec& ) );

    //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
    void set_values ( const epdf &ep );

    //! get value nearest to the given point
    double nearest_val ( const vec &val ) {
		NOT_IMPLEMENTED(0.0);
//        return values ( sup.linear_index ( sup.nearest_point ( val ) ) );
    }

	double integrate() ;

    vec eval ( const vec &val ) {
        return vec_1 ( nearest_val ( val ) );
    }
    const vec & _values() const {
        return values;
    }
};
UIREGISTER ( grid_fnc );

//! \brief Piecewise constant pdf on rectangular support
//!
//! Each point on the grid represents a centroid around which the density is constant.
//! This is a trivial point-mass density where all points have the same mass.
class egrid: public epdf {
protected:
    //! support of the pdf - grid
    rectangular_support sup;
    //! values at the grid
    vec values;
public:
    //! we assume that evallog is not called too often otherwise we should cache log(values)
    double evallog ( const vec &val ) {
        return log ( values ( sup.linear_index ( sup.nearest_point ( val ) ) ) );
    }
};
}
#endif //DISCR_H
