#ifndef __TREFSIM_HPP__
#define __TREFSIM_HPP__

#include <cmath>
#include <set>
#include <utility>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>
#include "tref.hpp"

namespace tref
{

template <class TREF> class sim {
protected:
	TREF *tref;
	gsl_rng *r;
	bool verbose;
	std::vector<double> objectQuality;
	std::vector<double> userAbility;
	unsigned long int objects;
	unsigned long int users;
	double links;
	double linksMax;
	double qualityScale;
	double abilityScale;
	double ratingScale;
	long unsigned int iter;

	virtual void generateObjectQuality()
	{
		for(std::vector<double>::iterator Q=objectQuality.begin();Q!=objectQuality.end();++Q)
			*Q = qualityScale * gsl_rng_uniform_pos(r);
	}

	virtual void generateUserAbility()
	{
		for(std::vector<double>::iterator sigma2=userAbility.begin();sigma2!=userAbility.end();++sigma2)
			*sigma2 = abilityScale * gsl_rng_uniform_pos(r);
	}

	virtual double rateObject(objectID o, userID u)
	{
		return ( ( ratingScale * userAbility[u] * ((2*gsl_rng_uniform(r))-1) ) + objectQuality[o] );
	}

	void generateRatings()
	{
		if(links!=linksMax) {
			std::vector<long unsigned int> objectLinks(objects,0);
			std::vector<long unsigned int> userLinks(users,0);
			std::set< std::pair<objectID,userID> > L;
			unsigned long int unlinkedObjects = 0, unlinkedUsers = 0;
			for(long unsigned int l=0;l!=links;++l) {
				std::pair<std::set< std::pair<objectID,userID> >::iterator,bool> ret;
				objectID o;
				userID u;
				do {
					o = gsl_rng_uniform_int(r,objects);
					u = gsl_rng_uniform_int(r,users);
					ret = L.insert(std::pair<objectID,userID>(o,u));
				} while(!(ret.second));
				tref->addRating(o,u,rateObject(o,u));
				++(objectLinks[o]);
				++(userLinks[u]);
				//std::cout << "Added a link between object " << o << " and user " << u << std::endl;
			}
			for(objectID o=0;o!=objects;++o) {
				if(objectLinks[o]==0) {
					//std::cout << "Object " << o << " has no links!" << std::endl;
					++unlinkedObjects;
				}
			}
			for(userID u=0;u!=users;++u) {
				if(userLinks[u]==0) {
					//std::cout << "User " << u << " has no links!" << std::endl;
					++unlinkedUsers;
				}
			}
			if(verbose && ((unlinkedObjects==0) || (unlinkedUsers==0))) {
				std::cout << unlinkedObjects << " objects and " << unlinkedUsers << " users ";
				std::cout << "have no links!" << std::endl;
			}
		} else {
			for(userID u=0;u!=users;++u)
				for(objectID o=0;o!=objects;++o)
					tref->addRating(o,u,rateObject(o,u));
		}
	}
public:
	sim(unsigned long int _objects, unsigned long int _users,
	    double sparsity, double beta, double epsilon, double convergence,
	    double _qualityScale, double _abilityScale, double _ratingScale,
	    gsl_rng *_r, bool _verbose)
	{
		verbose = _verbose;
		objects = _objects;
		users = _users;
		qualityScale = _qualityScale;
		abilityScale = _abilityScale;
		ratingScale = _ratingScale;
		objectQuality.resize(objects);
		userAbility.resize(users);
		r = _r;
		linksMax = ((double) objects)*((double) users);
		if(verbose)
			std::cout << "Maximum possible links: " << linksMax << std::endl;
		links = sparsity*linksMax;
		if(verbose)
			std::cout << "Actual number of links: " << links << std::endl;
		tref = new TREF(objects,users,links,beta,epsilon,convergence,verbose);
	}

	~sim()
	{
		delete tref;
	}

	void simulation()
	{
		generateObjectQuality();
		generateUserAbility();
		generateRatings();
		iter = tref->iterationCycle();
	}

	double qualityDelta()
	{
		double Delta = 0;
		for(objectID o=0;o!=objects;++o)
			Delta += pow((tref->objectQuality(o)-objectQuality[o]),2);
		return sqrt(Delta/objects);
	}

	long unsigned int iterations() { return iter; }
};


template <class TREF> class paretosim : public sim<TREF> {
	virtual double rateObject(objectID o, userID u)
	{
		double Cf = 2;
		if(this->ratingScale > 3)
			Cf = (this->ratingScale-2)*(this->ratingScale-3);
		double value = gsl_ran_pareto(this->r,this->ratingScale-1,sqrt(0.5*Cf*this->userAbility[u]));
		return (((gsl_rng_uniform(this->r)<0.5)?-value:value) + this->objectQuality[o]);
	}
public:
	paretosim(unsigned long int _objects, unsigned long int _users,
	          double sparsity, double beta, double epsilon, double convergence,
	          double _qualityScale, double _abilityScale, double _ratingScale,
	          gsl_rng *_r, bool _verbose) : sim<TREF> (_objects, _users,
                                                           sparsity, beta, epsilon, convergence,
                                                           _qualityScale, _abilityScale, _ratingScale,
                                                           _r, _verbose)
	{
	}
};


template <class TREF> class llsim : public sim<TREF> {
protected:
	virtual double rateObject(objectID o, userID u)
	{
		double value;
		do {
			value = sim<TREF>::rateObject(o,u);
		} while (value<0 || value > this->qualityScale);
		return value;
	}
public:
	llsim(unsigned long int _objects, unsigned long int _users,
	      double sparsity, double beta, double epsilon, double convergence,
	      double _qualityScale, double _abilityScale, double _ratingScale,
	      gsl_rng *_r, bool _verbose) : sim<TREF> (_objects, _users,
	                                               sparsity, beta, epsilon, convergence,
	                                               _qualityScale, _abilityScale, _ratingScale,
	                                               _r, _verbose)
	{
	}
};


template <class TREF> class intsim : public sim<TREF> {
protected:
	virtual double rateObject(objectID o, userID u)
	{
		double value = sim<TREF>::rateObject(o,u);
		value = floor(value+0.5);
		return value;
	}
public:
	intsim(unsigned long int _objects, unsigned long int _users,
	       double sparsity, double beta, double epsilon, double convergence,
	       double _qualityScale, double _abilityScale, double _ratingScale,
	       gsl_rng *_r, bool _verbose) : sim<TREF> (_objects, _users,
	                                                sparsity, beta, epsilon, convergence,
	                                                _qualityScale, _abilityScale, _ratingScale,
	                                                _r, _verbose)
	{
	}
};


template <template <class T> class TREFSIM, class TREF> class multisim {
protected:
	std::vector<double> qualityDelta;
	std::vector<double> iterations;
	double mean(std::vector<double> &v) { return gsl_stats_mean(&(v[0]),1,v.size()); }
	double median(std::vector<double> &v)
	{
		std::vector<double> w(v.begin(),v.end());
		gsl_sort(&(w[0]),1,w.size());
		return gsl_stats_median_from_sorted_data(&(w[0]),1,w.size());
	}
	double sd(std::vector<double> &v) { return gsl_stats_sd(&(v[0]),1,v.size()); }
	double min(std::vector<double> &v) { return gsl_stats_min(&(v[0]),1,v.size()); }
	double max(std::vector<double> &v) { return gsl_stats_max(&(v[0]),1,v.size()); }
public:
	multisim(unsigned long int sims,
	         unsigned long int objects, unsigned long int users,
	         double sparsity, double beta, double epsilon, double convergence,
	         double qualityScale, double abilityScale, double ratingScale,
	         gsl_rng *r, bool verbose)
	{
		long unsigned int iter_total = 0;
		qualityDelta.resize(sims);
		iterations.resize(sims);

		for(unsigned long int s=0;s!=sims;++s) {
			TREFSIM<TREF> sim(objects,users,
			                  sparsity,beta,epsilon,convergence,
			                  qualityScale,abilityScale,ratingScale,
			                  r,verbose);
			sim.simulation();
			qualityDelta[s] = sim.qualityDelta();
			iter_total += iterations[s] = sim.iterations();
			if(verbose)
				std::cout << "[" << s << "]\t" << qualityDelta[s] << "\t" << iterations[s] << std::endl;
		}
		if(verbose)
			std::cout << "Total iterations: " << iter_total << std::endl;
	}

	double qualityDeltaMean() { return mean(qualityDelta); }
	double qualityDeltaMedian() { return median(qualityDelta); }
	double qualityDeltaSD() { return sd(qualityDelta); }
	double qualityDeltaMin() { return min(qualityDelta); }
	double qualityDeltaMax() { return max(qualityDelta); }
	double iterationsMean() { return mean(iterations); }
	double iterationsMedian() { return median(iterations); }
	double iterationsSD() { return sd(iterations); }
	double iterationsMin() { return min(iterations); }
	double iterationsMax() { return max(iterations); }

	void stats()
	{
		std::cout << "\t" << qualityDeltaMean();
		std::cout << "\t" << qualityDeltaMin();
		std::cout << "\t" << qualityDeltaMax();
		std::cout << "\t" << iterationsMean();
		std::cout << "\t" << iterationsMin();
		std::cout << "\t" << iterationsMax() << std::endl;
	}
};

}

#endif /* __TREFSIM_HPP__ */