/* 
* 
* Copyright (C) 2013 UMPC
* 
* This file is part of CLaRiNet.
*  
* CLaRiNet is free software: you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* CLaRiNet is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with CLaRiNet.  If not, see <http://www.gnu.org/licenses/>.
*/
// Library of functions for Finite-Size Scaling analysis

#include <stdlib.h>
#include <stddef.h>
#include "Erreur.h"
#include "FiniteSizeScaling.h"
#include "math.h"
#include "float.h"
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

using namespace std;

int RecursionNoe ( int n, double * p , int * g , int * h , double * _Pn )
{
	int offg = 0;
	int offh = -1;
	int offp = -1;
	int offQ = 0;
	
	int m,i,k ;
	double * Q , * Qprev ;
	double logCik , logp , logTermk ;
	
	if ( !( Q = (double *) malloc( (n+1) * sizeof(double) ) ) )
	{
		    SetErreur("RecursionNoe : could not allocate Q\n");
            return(-1);
	}

	if ( !( Qprev = (double *) malloc( (n+1) * sizeof(double) ) ) )
	{
		    SetErreur("RecursionNoe : could not allocate Qprev\n");
            return(-1);
	}
	
	for ( i = 0 ; i < (n+1) ; i++ )
		Q[i] = 0 ;
	
	Q[0+offQ] = 1;
			
	for ( m = 1 ; m < (2*n+1) ; m++ )
	{
	
		// Mise à jour Qprev <- Q
		
		for ( i = 0 ; i < (n+1) ; i++ )
			Qprev[i] = Q[i] ;

		if ( p[m+offp] > 0 )
			logp = log(p[m+offp]) ;
			
		for ( i = h[m+1+offh]-1 ; i < (g[m-1+offg]+1) ; i++ )
		{
			Q[i+offQ] = Qprev[i+offQ] ; // = 0;
		
			if ( p[m+offp] > 0 )
			{		
				for ( k = h[m+offh]-1 ; k < i ; k++ ) //(i+1) ; k++ ) k = i déjà traité via l'initialisation
				{
					if ( Qprev[k+offQ] > 0 )
					{
						logCik = logCoeffBinom ( i , k ) ;						
						logTermk = logCik + log(Qprev[k+offQ]) + ((double)(i-k)) * logp ;					
						Q[i+offQ] += exp(logTermk) ;

					}
				}
			}				
		}	
	}
	
	(*_Pn) = Q[n+offQ];
	
	free(Q);
	free(Qprev);
	
	return(0);
}

int RecursionBohmMohanty ( int n, double * w , int * c , double * _Qn )
{
	// Même calcul que précédemment mais avec la méthode de Böhm et Mohanty (1997)
	
	int l,i ;
	int cl,ci;
	double res1,res2,logTermi,logC;
	double wl,xi;
	double * R;
	
	
	if ( !( R = (double *) malloc( (2*n+2) * sizeof(double) ) ) )
	{
		    SetErreur("RecursionBohmMohanty : could not allocate R\n");
            return(-1);
	}
	
	for ( i = 0 ; i<(2*n+2) ; i++)
		R[i] = 0 ;
			
	R[0] = -1;
	
	for ( l = 1 ; l < (2*n+2) ; l++ )
	{
		cl = c[l] ;
		wl = w[l];
		
		for ( i = 0 ; i < l ; i++ )
		{
			ci = c[i] ;

			if ( ci <= cl )
			{	
				if ( wl > 0 )
				{
					xi = w[i]/wl;
				}
				else
				{
					xi = 1.0;		// convention	(0/0) = 1
				}						
			
				if ( ! nLogx ( cl-ci , 1.0-xi , & res1 ) )
				{						
					if ( ! nLogx ( ci , xi , & res2 ) )
					{
						logC = logCoeffBinom ( cl , ci ) ;
						logTermi = logC + res1 + res2 ;
						R[l] -= exp(logTermi)*R[i] ;
					}
				}
			}
		}
		
		
	}
	
	(*_Qn) = R[2*n+1];
	
	free(R);
	
	return(0);
}

int RecursionBohmMohanty2 ( int n, int i0 , double L , double H , double * _Qn )
{
	// Version modifiée pour calculer la probabilité que la statistique d'ordre U{i0:n} soit comprise entre L et H
	// Ne calcule que les termes non-triviaux de la suite, les bornes a et b ne changeant de valeur qu'en i0
    /* 
           |------------ 1  	
	       |         **
	H -----|        *
	       |  ******
	      ****  
	     * | 
	   *   |------- L
	 **    | 
	0 ---- i0	
	
	*/
	
	int l,i ;
	double LogL,Log1_L,LogH,Log1_H,LogTermi,LogC;
	double * R;
	
	if ( !( R = (double *) malloc( (n+1) * sizeof(double) ) ) )
	{
		    SetErreur("RecursionBohmMohanty2 : could not allocate R\n");
            return(-1);
	}
	
	for ( i = 0 ; i < (n+1) ; i++)
		R[i] = 1 ;		

	LogH = log(H);		
	Log1_H = log(1.0-H);
	LogL = log(L);
	Log1_L = log(1.0-L);
		
	for ( l = i0 ; l < (n+1) ; l++ )
	{
	
		LogC = 0;
		
		for ( i = 0 ; i < i0 ; i++ )
		{
			LogTermi = LogC + (l-i)*Log1_H + i*LogH ;
			R[l] -= exp(LogTermi)*R[i] ;
			
			LogC += ( log(l-i) - log(i+1) );
		}
		
		for ( i = i0 ; i < (l+1) ; i++ )
		{
			LogTermi = LogC + (l-i)*Log1_L + i*LogL ;
			R[l] -= exp(LogTermi) ;
			
			LogC += ( log(l-i) - log(i+1) );
		}
		
	}
	
	(*_Qn) = R[n];
	
	free(R);
	
	return(0);
}


int nLogx ( int n , double x , double * _res )
{
	if (x<0 | n<0)
		return(-1);

	if (x>0)
	{
		*_res = n*log(x);
		return(0);
	}
	else
	{
		if (n==0)
		{	
			*_res = 0;
			return(0);
		}
		else
		{
			return(-1);
		}
	}
}

double logCoeffBinom ( int n , int k ) 
{
	// renvoie le logarithme naturel du coefficient binomial C(n,k)

	double r ;
	int i;
	
	r = 0;
	
	for ( i = k+1 ; i<(n+1) ; i++ )
		r += log((double)i);
		
	for ( i = 1 ; i<(n-k+1) ; i++ )
		r -= log((double)i);	
	
	return(r);
}

double BerkJonesLike ( double phat, double p )
{
	double K ; 
	
	if (p>0.0 && p<1.0)
	{
		K = phat * log(phat/p) + (1-phat)*log((1-phat)/(1-p)) ;
	}
	else
	{
		K = DBL_MAX ;
	}
	
	return(K);
}

int ZeroBJ ( double Fn , double lambda , double x1, double x2, double tol , int ITMAX , double * _res )
{
	int iter;
	double a=x1,b=x2,c=x2,d,e,min1,min2;
	double fa,fb,fc,p,q,r,s,tol1,xm;
	double EPS = 2.0*DBL_EPSILON ;

	fa = BerkJonesLike ( Fn , x1 ) - lambda ;
	fb = BerkJonesLike ( Fn , x2 ) - lambda ;
	
	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
	{
		SetErreur("ZeroBJ : x1 and x2 must bracket the solution\n");
		return(-1);
	}
		
	fc=fb;
	
	for (iter=1;iter<=ITMAX;iter++) {
		if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (fabs(fc) < fabs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*fabs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (fabs(xm) <= tol1 || fb == 0.0) 
		{
			(*_res)= b ;
			return(0);
		}
		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0) q = -q;
			p=fabs(p);
			min1=3.0*xm*q-fabs(tol1*q);
			min2=fabs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (fabs(d) > tol1)
			b += d;
		else
			b += SIGN(tol1,xm);
			
		fb = BerkJonesLike ( Fn , b ) - lambda ;
	}
	
	SetErreur("ZeroBJ : max number of iterations exceeded (%d)\n",ITMAX);
	return(-1);
}

int CoverageLevel ( int n , int i0 , double Fn , double lambda , int ITMAX , double tol , double * _L , double * _H , double * _cl )
{
	// Calcule la proba de couverture de la stat U{i0:n}, de fréquence empirique Fn, connaissant la valeur lambda
	// Renvoie également les bornes L,H associées
	
	double L,H,cl;
	
	// On commence par calculer les bornes L et H
	
/*	if ( Fn < 0.5)
		tol = (1.0D-3) * Fn ;
	else
		tol = (1.0D-3) * (1.0-Fn) ;*/
		
	
	if ( ZeroBJ ( Fn , lambda , 0.0 , Fn , tol , ITMAX , & L ) )
	{
		SetErreur("CoverageLevel : could not compute lower bound\n");
		return(-1);
	}
	
	if ( ZeroBJ ( Fn , lambda , Fn , 1.0 , tol , ITMAX , & H ) )
	{
		SetErreur("CoverageLevel : could not compute upper bound\n");
		return(-1);
	}
	
	// Calcul de la probabilité de couverture associée à ces bornes
	
	if ( RecursionBohmMohanty2 ( n , i0 ,  L , H , & cl ) )
	{	
		SetErreur("CoverageLevel : could not compute coverage probability\n");
		return(-1);
	}
	
	(*_L) = L ;
	(*_H) = H ;
	(*_cl) = cl;
	
	return(0);
}

int PointWiseConfidence ( int n , int i0 , double Fn , double alpha , double UBlambda , int ITMAX_BM , int ITMAX_BJ , double * _L , double * _H , double * _cl )
{
	// Calcule l'intervalle de confiance de la statistique d'ordre X{i0:n}, de fréquence empirique Fn, au niveau de couverture 1-alpha
	// Autorise le test de ITMAX_BM valeurs de lambda, et 2*ITMAX_BJ appels à ZeroBJ pour chaque lambda testé
	
	int iter;
	double x1,x2,a,b,c,d,e,min1,min2;
	double fa,fb,fc,p,q,r,s,tol,tol1,xm,TOL_BJ;
	double EPS = DBL_EPSILON ;

	x1 = 0;
	x2 = UBlambda ;	// majorant à fournir (?)
	
	a=x1;
	b=x2;
	c=x2;
		
	fa = alpha - 1.0 ;
	tol = EPS*alpha;
	
	if (Fn<0.5)
		TOL_BJ = EPS*Fn ;
	else
		TOL_BJ = EPS*(1.0-Fn) ;
	
	if ( CoverageLevel ( n , i0 , Fn , b , ITMAX_BJ , TOL_BJ , _L , _H , _cl ) )
	{
		SetErreur("PointWiseConfidence : could not compute coverage level, iter = %d, lambda = %e, TOL_BJ = %e\n",0,b,TOL_BJ);
		return(-1);
	}
	
	fb = (*_cl) + alpha - 1.0 ;
	
	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
	{
		SetErreur("PointWiseConfidence : UBlambda must be larger than the solution\n");
		return(-1);
	}
		
	fc=fb;
	
	for (iter=1;iter<=ITMAX_BM;iter++) {
		if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (fabs(fc) < fabs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*fabs(b)+0.5*tol;
		xm=0.5*(c-b);
		
		if (fabs(xm) <= tol1 || fb == 0.0) 
		{
			/* Retour de la solution ICI */
			// vide car on ne récupère pas la valeur de lambda, juste *_L , *_H et *_cl
			return(0);
			
		}
		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0) q = -q;
			p=fabs(p);
			min1=3.0*xm*q-fabs(tol1*q);
			min2=fabs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (fabs(d) > tol1)
			b += d;
		else
			b += SIGN(tol1,xm);
			
		/* Ré-évalutation à b */
		
		if ( CoverageLevel ( n , i0 , Fn , b , ITMAX_BJ , TOL_BJ , _L , _H , _cl ) )
		{
			SetErreur("PointWiseConfidence : could not compute coverage level, iter = %d, lambda = %e, TOL_BJ = %e\n",iter,b,TOL_BJ);
			return(-1);
		}
	
		fb = (*_cl) + alpha - 1.0 ;
		
	}
	
	SetErreur("PointWiseConfidence : max number of iterations exceeded (%d)\n",ITMAX_BM);
	return(-1);
}


