/*!
  \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
{

	//! Rectangular support
	//! Support ponits are located inbetween ranges! For example:
	//! For ranges=[0,1] and gridsizes=[1] the support point is 0.5
	class rectangular_support: public root {
		protected:
			//! Array of boundaries (2D vectors: [begining,end]) for each dimension
			Array<vec> ranges;
			//! Number of support points in each dimension
			ivec gridsizes;
			//! dimension
			int dim;
			//! Number of data points
			int Npoints;
			//! active vector for first_vec and next_vec
			vec actvec;
			//! indeces of active vector
			vec actvec_ind;
			//! length of steps in each dimension
			vec steps;
		public:
			//! default constructor
			rectangular_support(){};
			
			//! set parameters
			void set_parameters(const Array<vec> &ranges0, const ivec &gridsize0){
				ranges=ranges0;
				gridsizes=gridsize0;
				initialize();
			}
			//! Internal functio to set temporaries correctly
			void initialize() {
				dim = ranges.length();
				it_assert_debug(gridsizes.length()==dim,"Incompatible dimensions of input");
				Npoints = prod(gridsizes);
				it_assert_debug(Npoints>0,"Wrong input parameters");
				
				//precompute steps
				steps.set_size(dim);
				for ( int j = 0; j < dim; j++ ) {
					steps ( j ) = ( ranges ( j ) ( 1 ) - ranges(j)(0) ) / gridsizes ( j );
				}
				actvec.set_size(dim);
				actvec_ind.set_size(dim);
			}
			//! return vector at position given by vector of indeces
			vec get_vec(const ivec &inds){
				vec v;
				for ( int j = 0; j < dim; j++ ) {
					it_assert_debug(inds(j)<gridsizes(j), "Index out of bounds");
					v ( j ) = ranges(j)(0) + (0.5+inds(j))*steps(j);
				}
				return v;
			}

			//! convert dimension indeces into linear index, the indexing is in the same way as in \c next_vec()
			long linear_index (const ivec inds){
				long ind=0;
				it_assert_debug(inds.length()==dim,"Improper indeces inds");
				int dim_skips=1; // skips in active dimension, in the first dimension, the skips are 1 index per value
						
				for (int j=0; j<dim; j++){ 
					ind += dim_skips*(inds(j)); // add shift in linear index caused by this dimension
					dim_skips*=gridsizes(j);  // indeces in the next dimension are repeated with period gridsizes(j) times greater that in this dimesion
				}
				return ind;
			} 
			//! set the first corner to actvec
			const vec& first_vec(){
				for ( int j = 0; j < dim; j++ ) {
					actvec ( j ) = ranges(j)(0) + 0.5*steps(j);
					actvec_ind(j) = 0;
				}
				return actvec;
			}
			//! Get next active vector, call ONLY after first_vector()!
			const vec& next_vec() {
				// go through all dimensions
				int j=0;
				while (j<dim) {
					if ( actvec_ind ( j ) == gridsizes ( j ) - 1 ) { //j-th index is full
						actvec_ind ( j ) = 0; //shift back
						actvec ( j ) = ranges ( j ) ( 0 ) + 0.5*steps(j);
						j++;
					} else {
						actvec_ind ( j ) ++;
						actvec ( j ) += steps ( j );
						break;
					}
				}
				return actvec;
			} 
			
			ivec nearest_point(const vec &val){
				ivec inds;
				inds.set_size(dim);
				for(int j=0; j<dim; j++){
					if (val(j)<ranges(j)(0))
						inds(j) = 0;
					else {
						if (val(j)>ranges(j)(1))
							inds(j) = gridsizes(j)-1;
						else {
							inds(j) = ::round(val(j) - ranges(j)(0)/ steps(j));
						}
					}
				}
				return inds;
			}

			//! Access function
			int points() const {return Npoints;}
			
			//! access function
			const vec& _steps() const {return steps;}
			
			void from_setting (const Setting &set) {
				UI::get (ranges , set, "ranges", UI::compulsory);
				UI::get (gridsizes, set, "gridsizes", UI::compulsory);
				initialize();
			}
	};
	UIREGISTER(rectangular_support);
	
	class grid_fnc: public fnc{
		protected:
			rectangular_support sup;
			vec values;
		public:
			//! constructor function
			void set_support(rectangular_support &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&)){
				if (sup.points()>0){
					values(0) = (*evalptr)(sup.first_vec());
					for (int j=1; j<sup.points(); j++){ 
						values(j)=(*evalptr)(sup.next_vec());
					}
				}
			} 
			//! get value nearest to the given point
			double nearest_val(const vec &val){return values( sup.linear_index(sup.nearest_point(val)));}
			
			vec eval(const vec &val){return vec_1(nearest_val(val));}
	};
	UIREGISTER(grid_fnc);

	//! 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:
			rectangular_support sup;
			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