#include "example_fitter.h"

#include "iml_rand.h"
#include "mvl_mclk.h"

//===============================================================================
// main program:  example of how to use the generic fitter
//===============================================================================
int main(void) {

	// WE'LL GENERATE 10 VALUES FOR FITTING.
	IML_LINT DataCount = 10;
	
	// ALLOCATE A FITTER FOR MAX OF 10 PIECES OF DATA
	my_fitter DataFitter(DataCount);
	
	// GET A TRIAL COUNT FROM A HUMAN
	IML_LINT TrialCnt = 20;
	fprintf(stderr,"How many trials/category? "); fflush(stderr);
	scanf("%d",&TrialCnt);
	if(TrialCnt < 5) {
		fprintf(stderr,"%d trials is not enough. Using 5.\n",TrialCnt);
		TrialCnt = 5;
		} // end if
	
	// USE THESE FOR GENERATING DATA
	float Beta = 2.5;
	float Alpha = pow(10,-20/20);
	mvl_mac_clock Timer;
	uniform_rv UniformRV;
	UniformRV.seed(Timer.seed());
	
	// NOW LETS GENERATE THE DATA
	fprintf(stderr,"Generate %d trials/category using Alpha = %5.3f and Beta=%5.3f\n",
		TrialCnt, Alpha,Beta);
	for(long I=0; I<DataCount; I++) {
	
		// GET A HANDLE ON THIS DATA STRUCTURE FOR EASE OF ACCESS
		my_data &DP = DataFitter.Data[I];
		
		// CALCULATE THE CONTRAST
		DP.Contrast = pow(10,(-30.0+I*2.0)/20.0); // contrast in 2db steps
		
		// INITIALIZE THE TRIAL COUNT
		DP.TrialCount = TrialCnt;
		
		// CALCULATE THE EXPECTED PERCENT CORRECT FOR THE MODEL
		float ModelPc = weibull(DP.Contrast, Alpha, Beta);
		
		// GENERATE DATA USING BINOMIAL PROBABILITY
		DP.NumberCorrect = 0;
		for(long K=0; K<TrialCnt; K++)
			DP.NumberCorrect += (UniformRV.sample() <= ModelPc);
		
		fprintf(stderr,"For model Pc=%.3f NC=%.1f\n",ModelPc, DP.NumberCorrect);
		
		// THIS IS USED IN CALCULATING LOGLIKELIHOOD
		DP.NumberInCorrect = TrialCnt-DP.NumberCorrect;
		
		// WHAT WE WANT FOR THE NEXT STEP
		float DataPc = DP.NumberCorrect/(IML_DOUBLE)TrialCnt;
		DP.DataFractionCorrect = DataPc;
		
		// BINOMIAL VARIANCE FOR NumberCorrect IS TrialCnt*p*q
		// BECAUSE WE ARE USING % WE USE 1.0/(p*q/TrialCnt) = TrialCnt/(p*q) AS OUR WEIGHT
		float PQ = DataPc*(1-DataPc);
		
		// CHECK FOR POSSIBLE DIVIDE-BY-ZERO
		PQ = (PQ < 0.0001) ? 0.0001 : PQ;
		
		// WEIGHT IS 1/VARIANCE
		DP.ConfidenceWeight = TrialCnt/PQ;
		} // end for

	// SET CONVERGENCE CRITERIA:  TOLERANCE AND MAX-POWELL ITERATIONS
	DataFitter.set_convergence(1e-10, 30);
	
	// FAKE SOME GUESSES SO WE CAN WATCH THE FITTER CONVERGE.
	Alpha*=2.0;  Beta/=2.0;
	fprintf(stderr,"Begin fit guessing Alpha=%5.3f and Beta=%5.3f\n",Alpha,Beta);

	// CONFIGURE THE PARAMETERS - ID   NAME			MIN    	MAX   GUESS	OPTIMIZE-FLAG
	DataFitter.configure_parameter(0, "Weibull Alpha",	0.001,	0.99,	Alpha,	TRUE);
	DataFitter.configure_parameter(1, "Weibull Beta",	0.1,	10.0,	Beta,	TRUE);

	// CALCULATE THIS FIRST SO WE PERFORM THE FIT ON THE CHI-SQUARE SURFACE,
	// AND THE OPTIMUM VALUE WILL ACTUALLY BE CHI-SQUARE. THIS IS NOT
	// NECESSARY FOR THE FIT TO CONVERGE.
	IML_DOUBLE BestLogLikelihood = DataFitter.bestLL();
	
	// OPTIMUM VALUE IS ALSO CHI-SQUARE BECAUSE THATS OUR OBJECTIVE FUNCTION;
	IML_DOUBLE OptimumValue = DataFitter.fitdata();
	
	fprintf(stderr,"Evaluations =%d; LineMinimizations = %d; PowellIterations = %d;\n",
		DataFitter.Evaluations, 
		DataFitter.LineMinimizations, 
		DataFitter.PowellIterations
		);
	
	// PRINT THE CHISQUARE FIT
	fprintf(stderr,"Chisquare fit of %f with %d DOF\n",OptimumValue, DataCount-2);
	
	// GRAB THE OPTIMUM VALUES
	Alpha = DataFitter.parameter_value(0);
	Beta = DataFitter.parameter_value(1);
	
	// CALCULATE THE HESSIAN, COVARIANCE, AND CORRELATION MATRICES
	DataFitter.calculate_covariance_information();
	
	// GRAB THE PARAMETER STANDARD DEVIATIONS
	float AlphaSD = sqrt(DataFitter.covariance(0,0));
	float BetaSD  = sqrt(DataFitter.covariance(1,1));

	// PRINT THESE VITAL STATISTICS
	fprintf(stderr,"Alpha = %.3f +/- %.3f\n",Alpha, 2*AlphaSD);
	fprintf(stderr,"Beta  = %.3f +/- %.3f\n",Beta,  2*BetaSD);
	fprintf(stderr,"With correlation coefficient %.3f\n",DataFitter.correlation(0,1));
	
	
	// CALCULATE ALL THE FIT STATISTICS.
	DataFitter.calculate_fit_statistics();
	
	// PRINT MORE VITAL STATISTICS
	fprintf(stderr,"r^2: %5.3f   COD: %5.3f  MSC: %5.3f\n", 
		DataFitter.RSquared,
		DataFitter.CoeffOfDetermination,
		DataFitter.MSC
		);
	
	//==============================================================================
	// CALCULATE THE RIGOROUS 95% CHI-SQUARE LIMITS: UPPER/LOWER FOR BOTH PARAMETERS
	//==============================================================================
	IML_LINT MaxSteps = 15;
	
	// THIS IS THE CORRECT ERROR SURFACE CONTOUR FOR CHI-SQUARE.
	IML_DOUBLE Contour95 = OptimumValue + 3.84;
	float AMax = DataFitter.find_parameter_CI( 
		gml_generic_fitter::CI_HI,	// WHICH BOUNDARY TO FIND
		Contour95, 			// WHAT CONTOUR LEVEL TO FIND
		0, 				// PARAMETER ID
		Alpha, 				// LOWER GUESS BOUNDARY
		0.99, 				// UPPER GUESS BOUNDARY
		0.01, 				// ANSWER TOLLERANCE
		MaxSteps);			// MAX ITERATIONS TO TAKE
		
	float AMin = DataFitter.find_parameter_CI(
		gml_generic_fitter::CI_LO, Contour95, 0, 0.001, Alpha, 0.01, MaxSteps);

	float BMax = DataFitter.find_parameter_CI(
		gml_generic_fitter::CI_HI, Contour95, 1, Beta, Beta*3, 0.01, MaxSteps);

	float BMin = DataFitter.find_parameter_CI(
		gml_generic_fitter::CI_LO, Contour95, 1, 0.1, Beta/4, 0.01, MaxSteps);

	fprintf(stderr,"Rigorous Alpha Min,Max is %5.3f  %5.3f\n",AMin, AMax);
	fprintf(stderr,"Rigorous Beta  Min,Max is %5.3f  %5.3f\n",BMin, BMax);

return(0);
}