#ifndef __EXAMPLE_FITTER_H__
#define __EXAMPLE_FITTER_H__

#include <math.h>
#include "gml_gfit.h"


//=============================================================================
/** We use a Weibull as our psychometric function.
Our data go from 50% to 100% so we use 1.0-0.5*exp(-(C/A)^B) where
A is the 81% correct point and B determines the psychometric function slope.
*/
inline IML_DOUBLE weibull(IML_DOUBLE C, IML_DOUBLE A, IML_DOUBLE B) {
	return(1.0 - 0.5*exp( -1 * pow( C/A, B) ) );
	}

//=============================================================================
// FIRST WE DEFINE A DATA STRUCTURE WHICH WILL HOLD THE DATA AND 
// INTERMEDIATE MODEL VALUES.
//=============================================================================
/** The my_data class is used to contain example (psychometric) data.  
Calculated values are stored here as well.  The class knows how to
calculate error, in this case, negative log-likelihood, and how to
calculate the minimum log-likelihood so we can try to calculate
a chi-squared value once the best model is found.
*/
class my_data {

	public:
		/** The independent variable. */
		float Contrast;
		
		/** A parameter of the experiment. */
		float TrialCount;	
			
		/** The dependent variable. */
		float NumberCorrect; 		

		/** An alias for the dependent variable. */
		float NumberInCorrect; 
			
		/** Another alias for the dependent variable. */
		float DataFractionCorrect;
		
		/** The model estimate of the dependent variable. */
		float ModelFractionCorrect;

		/** Calculated as the binomial variance (estimate) for the data. */
		float ConfidenceWeight;
		
		/** This will hold the value of negative log-likelihood in the case
		where the model passes exactly through our DataFractionCorrect. */
		float BestLogLikelihood;

		/** This is the function that will be minimized for the model. 
		This will only be called by the evaluate() function which must
		be defined in my fitter below. */
		IML_DOUBLE logLikelihood(void)	 {
				return(-NumberCorrect*log(ModelFractionCorrect)
					- NumberInCorrect*log(1.0-ModelFractionCorrect) );  
			};  
		
		/** Calculate the "error" value in the case where the model
		passes right throught the data point. */
		void set_best_LL(void) { 
			if( !NumberCorrect || ! NumberInCorrect) 
				BestLogLikelihood= 0;
			else 
				BestLogLikelihood = 
					-NumberCorrect*log(DataFractionCorrect) 
					- NumberInCorrect*log(1.0-DataFractionCorrect ); 
			};  
		
};

//=============================================================================
//=============================================================================
/** The my_fitter class is an example of how to use the <em>gml_generic_fitter</em> class.
We will us it to fit a Weibull function to psychometric data.
The class is derived from the <em>gml_generic_fitter</em> class,
as described in the documentation.  All I need to do is tell 
<em>gml_generic_fitter</em> that I need two parameters, make sure 
to set the DataCount variable that I inherit from <em>gml_generic_fitter</em>, 
then initialize my own stuff.
*/
class my_fitter : public gml_generic_fitter {
	public:

		/** A pointer which holds an array of my_data structures. */
		my_data IML_PTR Data;
		
		/** A variable to hold an important quantity. */
		IML_DOUBLE BestLogLikelihood;
		
		/** The only constructor.  
		It tells <em>gml_generic_fitter</em> that there are two parameters, 
		and initializes the DataCount variable inherited from 
		<em>gml_generic_fitter</em>, and allocates the my_data array. */
		my_fitter(IML_USINT DCnt) : gml_generic_fitter(2) {

			DCnt = (DCnt > 0 ? DCnt : 1);
			DataCount = DCnt;
			Data = new my_data[DCnt];
			BestLogLikelihood = 0;
			}; // end

		/** The only destructor.  
		Just dallocates the my_data array. */
		~my_fitter() { if(Data) delete [] Data; }; 

		/** An interface function that must be defined.  
		Returns the model error for the point so that goodness-of-fit 
		statistics can be calculated. */
		virtual IML_DOUBLE model_error(register IML_ULINT I) 	
			{  return(Data[I].ModelFractionCorrect - Data[I].DataFractionCorrect); };
			
		/** An interface function that must be defined.  
		Returns the data value of the point so that goodness-of-fit 
		statistics can be calculated. */
		virtual IML_DOUBLE data(register IML_ULINT I) 	
			{  return(Data[I].DataFractionCorrect); };
			
		/** An interface function that must be defined.  
		Returns the model confidence weight for the point so that 
		goodness-of-fit statistics can be calculated. */
		virtual IML_DOUBLE weight(register IML_ULINT I) 
			{ return(Data[I].ConfidenceWeight); };

			
		/** The only interface function that must <b>absolutely</b> be defined.  
		Returns the model error upon which the fitter shall operate, and
		from which covariance and correlation matrices will eventually be
		calcualted. */
		virtual IML_DOUBLE evaluate( register gml_double_vec &PVector) { 
		
			Evaluations++;
		
			// LOOP OVER ALL THE DATA AND CALCULATE THE MODEL ESTIMATE
			for(register IML_LINT Index=0; Index < DataCount; Index++) {
				my_data &DP = Data[Index];
				DP.ModelFractionCorrect = weibull(DP.Contrast, PVector[0], PVector[1]);				
				DP.ModelFractionCorrect = GML_PFLOOR(0.001,DP.ModelFractionCorrect);
				DP.ModelFractionCorrect = GML_PCEIL(0.999,DP.ModelFractionCorrect);
				} // end for
				
			// CALCULATE THE LOGLILIHOOD
			register IML_DOUBLE LL=0;
			for(register IML_LINT Index=0; Index < DataCount; Index++)
				LL += Data[Index].logLikelihood();
				
			// AND RETURN THE (PUTATIVE) CHI-SQUARE VALUE
			return( 2*(LL-BestLogLikelihood) ); 
			}; // end evaluate

		/** Calculate the absolute minimum best fit. */
		IML_DOUBLE bestLL(void) { 
			BestLogLikelihood=0;
			for(register IML_LINT Index=0; Index < DataCount; Index++) {
				Data[Index].set_best_LL();
				BestLogLikelihood += Data[Index].BestLogLikelihood;
				} // END FOR
			return(BestLogLikelihood); 
			}; // end bestLL

};


#endif // __EXAMPLE_FITTER_H__