<commit_before>// SVN date:    $Date$
// SVN author:  $Author$
// SVN rev.:    $Revision$
// SVN URL:     $HeadURL$
// SVN ID:      $Id$

/*	GenCurveFit.c -- An XOP for curvefitting via Differential Evolution.
 See:
 Wormington, et. al., "Characterisation of structures from X-ray Scattering
 Data using Genetic Algorithms", Phil. Trans. R. Soc. Lond. A (1999) 357, 2827-2848
 
 And
 The Motofit packages: http://motofit.sourceforge.net/
 Nelson, "Co-refinement of multiple-contrast neutron/X-ray reflectivity data using Motofit",
 J. Appl. Cryst. (2006). 39,273-276.
 
 @copyright: <NAME> and the Australian Nuclear Science and Technology Organisation 2007.
 
 */

/*
 put all the standard headers in GenCurveFit.h, including those from wavemetrics
 */
#include "XOPStandardHeaders.h"
#include "GenCurveFitXOP.h"
#include "gencurvefit.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "updateXOP.h"


//gTheWindow is a window created to show the latest position of the fit
XOP_WINDOW_REF gTheWindow = NULL;

/*
 a variable which sees if the fit is to be aborted.
 */
int Abort_the_fit = 0;

/*
 some memory utilities
 */

//create a platform independent routine for continuous reallocation of memory, appending data to it
void *myrealloc(void *src_ptr, size_t size)
{
    /* There might be a realloc() out there that doesn't like reallocing
	 NULL pointers, so we take care of it here */
    if(src_ptr)
		return realloc(src_ptr, size);
    else
		return malloc(size);
}

size_t
WriteMemoryCallback(void *ptr, size_t size, size_t nmemb, void *data)
{
    size_t realsize = size * nmemb;
    struct MemoryStruct *mem = (struct MemoryStruct *)data;
	
    mem->memory = (char *)myrealloc(mem->memory, mem->size + realsize + 1);
    if (mem->memory) {
		memcpy(&(mem->memory[mem->size]), ptr, realsize);
		mem->size += realsize;
		mem->memory[mem->size] = 0;
    }
    return realsize;
}

/*
 the fitfunction for the genetic optimisation
 */
int lgencurvefit_fitfunction(void *userdata, const double *coefs, unsigned int numcoefs, double *model, const double **xdata, long datapoints, unsigned int numDataDims){
	int err = 0;
	
	CountInt ii;
	unsigned int jj;
	int requiredParameterTypes[MAX_MDFIT_SIZE + 2];
	allFitFunc allParameters;
	fitFunc parameters;
	double result;
	double *dp;
	GenCurveFitInternals *goiP = (GenCurveFitInternals*) userdata;

	memset(&parameters, 0, sizeof(parameters));
	memset(&allParameters, 0, sizeof(allParameters));
	

	//cmd-dot or abort button
	if(CheckAbort(0) == -1){
		err = FIT_ABORTED;
		goto done;
	}

	/*
	 check if all the input isn't NULL
	 */
	if(goiP->coefs == NULL || goiP->dataCalc == NULL){
		err = UNSPECIFIED_ERROR;
		goto done;
	}
	
	/*
	 Place the parameters into the coefficient wave
	 */
	if(err = MDStoreDPDataInNumericWave(goiP->GenCurveFitCoefs, (double*) coefs))
		goto done;
			
	switch(goiP->isAAO){
		case 0:
			parameters.waveH = goiP->GenCurveFitCoefs;
			requiredParameterTypes[0] = WAVE_TYPE;
			
			for(ii = 0 ; ii < datapoints ; ii += 1){
				for(jj = 0 ; jj < numDataDims ; jj += 1)
					parameters.x[jj] = xdata[jj][ii];
				
				// call the users fit function and put the result in the output array
				if (err = CallFunction(&goiP->fi, (void*) &parameters, &result))
					goto done;
				model[ii] = result;
			}
						
			break;
		case 1:
			allParameters.waveC = goiP->GenCurveFitCoefs;
			allParameters.waveY = goiP->dataCalc;
			requiredParameterTypes[0] = WAVE_TYPE;
			requiredParameterTypes[1] = WAVE_TYPE;
			
			for(ii = 0 ; ii < numDataDims ; ii += 1){
				requiredParameterTypes[ii + 2] = WAVE_TYPE;
				if(goiP->xcalc[ii] == NULL){
					err = UNSPECIFIED_ERROR;
					goto done;
				}
				allParameters.waveX[ii] = goiP->xcalc[ii];
			}

			// call the users fit function and put the result in the output wave
			if (err = CallFunction(&goiP->fi, (void*)&allParameters, &result))
				goto done;

			// the user may have changed the number of points in the output wave
			if(goiP->dataCalc == NULL || WavePoints(goiP->dataCalc ) != datapoints){
				err = USER_CHANGED_FITWAVE;
				goto done;
			}
			dp = WaveData(goiP->dataCalc);
			memcpy(model, dp, sizeof(double) * datapoints);
			
			break;
		case 2:			
			if(goiP->sp == NULL){
				err = NULL_STRUCTURE;
				goto done;
			}
			goiP->sp->w = goiP->GenCurveFitCoefs;
			goiP->sp->yy = goiP->dataCalc;
			
			
			for(ii = 0 ; ii < numDataDims ; ii+=1){
				if(goiP->xcalc[ii] == NULL){
					err = UNSPECIFIED_ERROR;
					goto done;
				}
				goiP->sp->xx[ii] = goiP->xcalc[ii];
			}

			if (err = CallFunction(&goiP->fi, (fitfuncStruct*)&goiP->sp, &result))
				goto done;
			//don't want any dangling references to waves
			goiP->sp->yy = NULL;
			goiP->sp->w = NULL;
			memset(goiP->sp->xx, 0, sizeof(goiP->sp->xx));
			
			// the user may have changed the number of points in the output wave
			if(goiP->dataCalc == NULL || WavePoints(goiP->dataCalc) != datapoints){
				err = USER_CHANGED_FITWAVE;
				goto done;
			}
			
			dp = WaveData(goiP->dataCalc);
			memcpy(model, dp, sizeof(double) * datapoints);
			
			break;
		default:
			err = UNSPECIFIED_ERROR;
			goto done;
			break;
	}
	
	// check that the fitfunction didn't return any NaN or INF
	if(err = checkNanInfArray(model, datapoints)){
		err = FITFUNC_RETURNED_NANINF;
		goto done;
	}


done:		

	return err;
}

/*
 the cost function for the genetic optimisation
 */
double lgencurvefit_costfunction(void *userdata, const double *coefs, unsigned int numcoefs, const double *data, const double *model, const double *errors, long datapoints){
	double val;

	GenCurveFitInternals *goiP = (GenCurveFitInternals*) userdata;
	costFunc userCostFunc;
	double *dp;
	int err;
	
	switch(goiP->METH){//which cost function
		case 0:
			val = chisquared(userdata, coefs, numcoefs, data, model, errors, datapoints);			
			break;
		case 1:
			val = robust(userdata, coefs, numcoefs, data, model, errors, datapoints);
			break;
		case 2:
			userCostFunc.coefs = goiP->GenCurveFitCoefs;
			userCostFunc.yobs =	 goiP->yobs;
			userCostFunc.ycalc = goiP->dataCalc;
			userCostFunc.sobs =  goiP->sobs;
			
			//copy the original data into the yobs wave created for the purpose
			//some sneaky users probably try to change it.
			dp = WaveData(goiP->yobs);
			memcpy(dp, data, datapoints * sizeof(double));
			
			//copy the original data into the yobs wave created for the purpose
			//some sneaky users probably try to change it.
			dp = WaveData(goiP->sobs);
			memcpy(dp, errors, datapoints * sizeof(double));

			if(err = CallFunction(&goiP->minf, (void*) &userCostFunc, &val))
				goto done;
			
			if(WavePoints(goiP->yobs) != datapoints || WavePoints(goiP->sobs) != datapoints || WavePoints(goiP->GenCurveFitCoefs)!= numcoefs){
				err = COSTFUNC_WAVES_CHANGED;
				goto done;
			}
			
			if(IsNaN64(&val) || IsINF64(&val)){
				err = COSTFUNC_DOESNT_RETURN_NUMBER;
				goto done;
			}
			
			break;
		default:
			val = chisquared(userdata, coefs, numcoefs, data, model, errors, datapoints);			
			break;
	}
	done:
			
	return val;
};


/*
 the update function for the genetic optimisation
 */
int lgencurvefit_updatefunction(void *userdata,
								const double *coefs,
								 unsigned int numcoefs,
								 unsigned int iterations,
								 double cost,
								 unsigned int updatetime,
								 double convergenceNumber,
								 const double** population,
								 const unsigned int *varparams,
								 unsigned int numvarparams,
								 unsigned int totalpopsize,
								 const double *costmap){
	int err = 0;
	double val = 0;
	
	GenCurveFitInternals *goiP = (GenCurveFitInternals*) userdata;
	updtFunc userUpdateFunc;
	
	goiP->V_numfititers = iterations;
	
	/*
	 Display the coefficients so far.
	 */
	if(!goiP->noupdate && (updatetime == 1 || updatetime == 8))
		DisplayWindowXOP1Message(gTheWindow, numcoefs, coefs, cost, goiP->fi.name, iterations, updatetime, convergenceNumber);

	// perhaps the user wants to abort the fit using gui panel?
	if(Abort_the_fit){
		err = FIT_ABORTED;
		goto done;
	}

	//We want a record of all the improvements in the fit
	if(goiP->dump && (updatetime == 2 || updatetime == 1)){
		WriteMemoryCallback((void*) coefs, sizeof(double), goiP->totalnumparams, &(goiP->dumpRecord));
		if(goiP->dumpRecord.memory == NULL){
			err = NOMEM; goto done;
		}
	}
	
	//store the best fit coefficients in the original coef wave and update the output data.
	if((updatetime == 1 && !goiP->noupdate) || updatetime == 16){
		goiP->cost = cost;

		if(err = MDStoreDPDataInNumericWave(goiP->OUT_coefs, (void*) coefs))
			return err;
		WaveHandleModified(goiP->OUT_coefs);
		
		if(err = calcModelXY(&goiP->fi, goiP->OUT_coefs, goiP->OUT_data, goiP->OUT_x, goiP->numVarMD, goiP->isAAO, goiP->sp))
			return err;
		
		WaveHandleModified(goiP->OUT_data);
		
		if(goiP->OUT_res){
			if(err = calcModelXY(&goiP->fi, goiP->OUT_coefs, goiP->OUT_res, goiP->fullExtentOfData, goiP->numVarMD, goiP->isAAO, goiP->sp))
				return err;
			if(err = subtractTwoWaves(goiP->OUT_res, goiP->originalYobs))
				return err;
			if(err = scalarMultiply(goiP->OUT_res, -1))
				return err;
			WaveHandleModified(goiP->OUT_res);
		}
		DoUpdate();		
	}
	
	//the user defined a user IGOR update function so call that.
	if(!goiP->noupdate && goiP->useIgorUpdateFunction && population && costmap){
		userUpdateFunc.currentbestfit = goiP->OUT_coefs;
		userUpdateFunc.population = goiP->M_population;
		userUpdateFunc.costmap = goiP->W_costmap;
		userUpdateFunc.updatetime = updatetime;

		memcpy(WaveData(goiP->M_population), &(population[0][0]), totalpopsize * numvarparams * sizeof(double));
		memcpy(WaveData(goiP->W_costmap), costmap, totalpopsize * sizeof(double) );

		WaveHandleModified(goiP->M_population);
		WaveHandleModified(goiP->W_costmap);

		if(err = CallFunction(&goiP->igorUpdateFunction, (void*) &userUpdateFunc, &val))
			goto done;
		if(val)
			err = (int)val;
	}
	
	
done:
	return err;
}

/* 
 function obtains the covariance matrix for the fit.
 */
int getGCovarianceMatrix(GenCurveFitRuntimeParamsPtr p, GenCurveFitInternalsPtr goiP){
	int err = 0;
	double hessianDeterminant = 0;
	
	if(err = getCovarianceMatrix(&goiP->covarianceMatrix,
								 &hessianDeterminant,
								 goiP,
								 lgencurvefit_fitfunction,
								 goiP->cost,
								 goiP->coefs,
								 (int) goiP->totalnumparams,
								 goiP->holdvector,
								 goiP->dataObs,
								 goiP->dataSig,
								 (const double**) goiP->independentVariable,
								 (long) goiP->unMaskedPoints,
								 (int) goiP->numVarMD,
								 !p->WFlagEncountered))
		return err;
				
	goiP->V_logBayes = exp(-0.5 * (goiP->cost) / (double)(goiP->unMaskedPoints - goiP->numvarparams));// * pow(4*3.14159,(double) goiP->numvarparams) ;//* factorial((double)goiP->numvarparams);
	goiP->V_logBayes = goiP->V_logBayes / (sqrt(hessianDeterminant));
	goiP->V_logBayes = log(goiP->V_logBayes);

	return err;
}

/*
 converts a libgencurvefit error to an XOP error
 All libgencurvefit errors are negative, all XOP errors are positive
 */

int convertlgencurvefitErrorToXOP(int lgencurvefiterror){
	int err = lgencurvefiterror;
	
	if(lgencurvefiterror > -1)
		return err;
	
	switch (lgencurvefiterror) {
		case NO_MEMORY:
			err = NOMEM;
			break;
		case INCORRECT_LIMITS:
			err = LIMITS_INVALID;
			break;
		case HOLDVECTOR_COEFS_MISMATCH:
			err = HOLDSTRING_NOT_RIGHT_SIZE;
			break;
		case NO_VARYING_PARAMS:
			err = ALL_COEFS_BEING_HELD;
			break;
		case WRONG_NUMBER_OF_PARAMS:
			break;
		case COEFS_MUST_BE_WITHIN_LIMITS:
			err = LIMITS_INVALID;
			break;
		case PROBLEM_CALCULATING_COVARIANCE:
			err = SINGULAR_MATRIX;
			break;
		case NO_FIT_FUNCTION_SPECIFIED:
			err = FITFUNC_NOT_SPECIFIED;
			break;
		case NO_Y_ARRAY:
			break;
		case NO_X_ARRAY:
			break;
		case NO_E_ARRAY:
			break;
		case NO_COEFS_ARRAY:
			break;
		case NO_LIMITS_ARRAY:
			err = LIMITS_INVALID;
			break;
		case SINGULAR_MATRIX_ERROR:
			err = SINGULAR_MATRIX;
			break;
		default:
			break;
	}
	
	
	return err;
}


/*
 ExecuteGenCurveFit performs the genetic curvefitting routines
 returns 0 if no error
 returns errorcode otherwise
 */
int
ExecuteGenCurveFit(GenCurveFitRuntimeParamsPtr p)
{
	/* the function that is called by IGOR. Here's where we farm the work out to different places.*/
	/*
	 GenCurveFitInternals carries the internal data structures for doing the fitting.
	 */
	GenCurveFitInternals goi;
	//initialise the structure to be zero
	
	/*
	 err carries the errors for all the operations
	 err2 carries the error code if we still want to return err, but we need to finish off
	 something else first.	
	 */
	int err = 0, err2 = 0;
	double value[2];
	CountInt indices[MAX_DIMENSIONS];
	
	//libgencurvefit options
	gencurvefitOptions gco;
	
	//variables listed below are purely for outputs sake.
	char varname[MAX_OBJ_NAME + 1];
	CountInt dimensionSizes[MAX_DIMENSIONS + 1];
	double t1,t2, chi2;
	double *wP;
	long lt1 = 0;
	char note[200], note_buffer1[MAX_WAVE_NAME + 1], note_buffer2[MAX_WAVE_NAME + 1], cmd[MAXCMDLEN + 1];
	int output, ii, jj, isDisplayed, quiet;
	
	//initialise all the internal data structures to NULL
	memset(&goi, 0, sizeof(goi));
	
	//you have to use IGOR > 6.10
	if( igorVersion < 610 )
		return REQUIRES_IGOR_610;
	
	//reset the options for libgencurvefit
	memset(&gco, 0, sizeof(gencurvefitOptions));
	
	//reset the abort condition
	Abort_the_fit = 0;
	
	//can't do this in a threadsafe manner if we're not running IGOR >6.2
	if (igorVersion < 620 && !RunningInMainThread())
		return IGOR_OBSOLETE;
	
	strncpy(varname, "V_Fiterror", MAX_OBJ_NAME);
	if(FetchNumVar(varname, &t1, &t2)!=-1){
		if(!err)
			lt1 = 0;
		else 
			lt1 = 1;
		
		err = 0;
		if(err2 = SetIgorIntVar(varname, lt1, 1)){err = err2;goto done;};
	}	
	
	/*
	 Genetic Optimisation uses a lot of random numbers, we are seeding the generator here.
	 If you want to seed the generator we can 
	 */
	if(p->SEEDFlagEncountered)
		gco.seed = (int) p->SEEDFlag_seed;
	else 
		gco.seed = -1;
	
	/*
	 checkInput checks the input that IGOR sends to the XOP.  If everything is correct it returns 0, 
	 else it returns an error.  Errors can be caused by, e.g. different x and y wave lengths, etc.
	 */
	if(err = checkInput(p, &goi))
		goto done;
	
	/*
	 init_GenCurveFitInternals sets up the internal data arrays in the GenCurveFitInternals structure.  It holds
	 copies of the datapoints being fitted, arrays of the fitting structures, etc.  It calls malloc
	 several times to set aside memory for all this.  If this procedure works without a hitch then the function
	 returns 0.
	 */
	if(err = init_GenCurveFitInternals(p, &goi))
		goto done;
	
	if(p->OPTFlagEncountered && (((long)p->OPTFlag_opt) & (long)pow(2, 0)))
		gco.useinitialguesses = 1;
	
	
	gco.updatefun = (updatefunction) &lgencurvefit_updatefunction;
	gco.k_m = goi.k_m;
	gco.recomb = goi.recomb;
	gco.strategy = goi.STGY;
	gco.tolerance = goi.tolerance;
	gco.popsizeMultiplier = goi.popsize;
	gco.iterations = goi.iterations;
	gco.temp = goi.temperature;
	gco.updatefrequency = 31;
	
	/*
	 you want to polish the fit (at the end) using LevenbergMarquardt
	 */
	if(p->POLFlagEncountered)
		gco.polishUsingLM = 1;
	
	/*
	 if you want to do a Monte Carlo fit.
	 */
	if(p->MCFlagEncountered && p->WFlagEncountered)
		gco.monteCarlo = 1;
	
	//make an error wave before we start the fit, to make sure it's always there
	dimensionSizes[0] = goi.totalnumparams;
	dimensionSizes[1] = 0;
	if(err = MDMakeWave(&goi.W_sigma, "W_sigma", goi.cDF, dimensionSizes, NT_FP64, 1))
			goto done;
	for(ii = 0 ; ii < goi.totalnumparams ; ii+=1){
		indices[0] = ii;
		value[0] = 0;
		if(err = MDSetNumericWavePointValue(goi.W_sigma, indices, value))
			goto done;
	}
	WaveHandleModified(goi.W_sigma);
		
	if(p->MATFlagEncountered && !(p->MATFlagParamsSet[0] && (int) p->MATFlag_mat == 0)){
		//make the covariance matrix
		dimensionSizes[0] = goi.totalnumparams;
		dimensionSizes[1] = goi.totalnumparams;
		dimensionSizes[2] = 0;
		if(err = MDMakeWave(&goi.M_covariance, "M_Covar", goi.cDF, dimensionSizes, NT_FP64, 1))
				goto done;
			
		wP = WaveData(goi.M_covariance);
			
		memset(wP, 0, sizeof(double) * goi.totalnumparams * goi.totalnumparams);
			
		WaveHandleModified(goi.M_covariance);
	}
	
	/*
	 optimiseloop does the Differential Evolution, according to Storn and Price.  When this returns 0, then 
	 you have the best fit in the GenCurveFitInternals structure.  Otherwise it returns an error code.  If the user aborts
	 then the FIT_ABORTED error code is returned, but it is still possible to retrieve the best solution so far
	 */
	if(err = genetic_optimisation(lgencurvefit_fitfunction,
							   lgencurvefit_costfunction,
							   (unsigned int) goi.totalnumparams,
							   goi.coefs,
							   goi.holdvector,
							   (const double**) goi.limits,
							   (long) goi.unMaskedPoints,
							   goi.dataObs,
							   (const double**) goi.independentVariable,
							   goi.dataSig,
							   goi.numVarMD,
							   &chi2,
							   &gco,
							   &goi))
		goto done;
	
	//make sure the coefficients for the fit are updated, by calling the update function
	Abort_the_fit = 0;
	if(err = lgencurvefit_updatefunction(&goi,
										 goi.coefs,
										 (unsigned int) goi.totalnumparams,
										 goi.V_numfititers,
										 chi2,
										 16,
										 1,
										 NULL,
										 NULL,
										 0,
										 0,
										 NULL))
	   goto done;
	
	//output the dumprecord
	if(!err && goi.dump && goi.dumpRecord.memory)
		if(err2 = dumpRecordToWave(&goi, &goi.dumpRecord)){
			err = err2;
			goto done;
		}
	
	/*
	 if there are no errors, or if the user aborted, then return the best fit.
	 If the data is displayed in the top graph append the fitcurve to the top graph
	 */
	if((err == 0 || err == FIT_ABORTED) && !getGCovarianceMatrix(p, &goi) && goi.covarianceMatrix){
		//set the error wave
		jj = 0;
		for(ii = 0; ii < goi.totalnumparams ; ii += 1){
			indices[0] = ii;
			if(!goi.holdvector[ii]){
				value[0] = sqrt(goi.covarianceMatrix[goi.varParams[jj]][goi.varParams[jj]]);
				jj++;
			} else {
				value[0] = 0;
			}

			if(err = MDSetNumericWavePointValue(goi.W_sigma, indices, value))
				goto done;
			
		}
		WaveHandleModified(goi.W_sigma);
		
		if(p->MATFlagEncountered && !(p->MATFlagParamsSet[0] && (int) p->MATFlag_mat == 0)){
			//make the covariance matrix
			dimensionSizes[0] = goi.totalnumparams;
			dimensionSizes[1] = goi.totalnumparams;
			dimensionSizes[2] = 0;
			wP = WaveData(goi.M_covariance);			
			memcpy(wP, *goi.covarianceMatrix, sizeof(double) * goi.totalnumparams * goi.totalnumparams);
			WaveHandleModified(goi.M_covariance);
		}
	}
	//append the fit to the graph
	if(!err || err == FIT_ABORTED){
		if(err = isWaveDisplayed(p->dataWave.waveH, &isDisplayed))
			goto done;
		
		if(isDisplayed && goi.numVarMD == 1){
			if(err = isWaveDisplayed(goi.OUT_data, &isDisplayed))
				goto done;
			
			if(!isDisplayed){
				strncpy(cmd, "appendtograph/w=$(winname(0,1)) ", MAXCMDLEN);
				WaveName(goi.OUT_data,&note_buffer1[0]);
				strncat(cmd, &note_buffer1[0], MAXCMDLEN - strlen(note_buffer1));
				
				if(p->DFlagEncountered && p->XFlagEncountered){
					WaveName(p->XFlag_xx,&note_buffer2[0]);
					strncat(cmd, " vs ", MAXCMDLEN - strlen(cmd) - strlen(" vs "));
					strncat(cmd, note_buffer2, MAXCMDLEN - strlen(cmd) - strlen(note_buffer2) );
				}
				if(err = XOPSilentCommand(&cmd[0]))
					goto done;
			}
		}
	}
	
	/*this section sets history and global variables
	 V_Chisq
	 V_Fiterror
	 If there are no errors returned above, then V_Chisq is produced in the current datafolder.
	 V_fiterror is for if there is a problem with the fitting.  
	 If this global variable exists in IGOR at runtime, and there is no XOPerror then V_fiterror = 0.  
	 If there is an XOPerror, bit 0 of V_fiterror is set. However, the XOP will return gracefully to IGOR,
	 allowing the user that called it to detect this and carry on.
	 In normal situations V_fiterror may not exist, so the GenCurveFit XOP will return an error message.
	 */
	
	if(!err){
		SetOperationNumVar("V_Chisq", chi2);
		SetOperationNumVar("V_fitIters", (int)(goi.V_numfititers));
		SetOperationNumVar("V_npnts", (int)(goi.unMaskedPoints));
		SetOperationNumVar("V_nterms", (int)WavePoints(p->coefs));
		SetOperationNumVar("V_nheld", (int)(WavePoints(p->coefs) - goi.numvarparams));
		SetOperationNumVar("V_logBayes", goi.V_logBayes);
	}		
	
	strcpy(varname, "V_Fiterror");
	if(FetchNumVar(varname, &t1, &t2)!=-1){
		if(!err){
			lt1 = 0;
		} else {
			lt1 = 1;
		}
		err = 0;
		if(err2 = SetIgorIntVar(varname, lt1, 1))
		{err = err2;goto done;}
	}
	/*
	 This section puts a copy of the fit parameters into the history area, unless one sets quiet mode.
	 */
	quiet = 0;
	if(p->QFlagEncountered){
		quiet = 1;
		if(p->QFlagParamsSet[0] && (int) p->QFlag_quiet == 0)
			quiet = 0;
	}
	
	if(!quiet && (!err || err == FIT_ABORTED) && lt1==0 ){		
		if(!err)
			{output = snprintf(note, 199, "_______________________________\rGenetic Optimisation Successful\r");XOPNotice(note);}
		if(err == FIT_ABORTED)
			{output = snprintf(note, 199, "_______________________________\rGenetic Optimisation ABORTED\r");XOPNotice(note);}
		WaveName(p->dataWave.waveH, note_buffer1);
		
		output = snprintf(note, 
						  199,
						  "Fitting: %s to %s\r",
						  note_buffer1,
						  goi.fi.name);
		
		XOPNotice(note);
		
		output = snprintf(note,
						  199,
						  "V_fitIters = %li; V_Chisq = %g; V_npnts= %li; V_nterms= %li; V_nheld= %li; V_logBayes = %g\r",
						  goi.V_numfititers,
						  chi2,
						  goi.unMaskedPoints,
						  WavePoints(p->coefs),
						  WavePoints(p->coefs) - goi.numvarparams,
						  goi.V_logBayes);
		
		XOPNotice(note);
		
		for(ii = 0; ii < WavePoints(p->coefs) ; ii += 1){
			indices[0] = ii;
			indices[1] = 0;
			if(err = MDGetNumericWavePointValue(goi.W_sigma, indices, value))
				goto done;
			
			output = snprintf(note,<commit_msg>changing over to toolkit 6<commit_after>// SVN date:    $Date$
// SVN author:  $Author$
// SVN rev.:    $Revision$
// SVN URL:     $HeadURL$
// SVN ID:      $Id$

/*	GenCurveFit.c -- An XOP for curvefitting via Differential Evolution.
 See:
 Wormington, et. al., "Characterisation of structures from X-ray Scattering
 Data using Genetic Algorithms", Phil. Trans. R. Soc. Lond. A (1999) 357, 2827-2848
 
 And
 The Motofit packages: http://motofit.sourceforge.net/
 Nelson, "Co-refinement of multiple-contrast neutron/X-ray reflectivity data using Motofit",
 J. Appl. Cryst. (2006). 39,273-276.
 
 @copyright: <NAME> and the Australian Nuclear Science and Technology Organisation 2007.
 
 */

/*
 put all the standard headers in GenCurveFit.h, including those from wavemetrics
 */
#include "XOPStandardHeaders.h"
#include "GenCurveFitXOP.h"
#include "gencurvefit.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "updateXOP.h"


//gTheWindow is a window created to show the latest position of the fit
XOP_WINDOW_REF gTheWindow = NULL;

/*
 a variable which sees if the fit is to be aborted.
 */
int Abort_the_fit = 0;

/*
 some memory utilities
 */

//create a platform independent routine for continuous reallocation of memory, appending data to it
void *myrealloc(void *src_ptr, size_t size)
{
    /* There might be a realloc() out there that doesn't like reallocing
	 NULL pointers, so we take care of it here */
    if(src_ptr)
		return realloc(src_ptr, size);
    else
		return malloc(size);
}

size_t
WriteMemoryCallback(void *ptr, size_t size, size_t nmemb, void *data)
{
    size_t realsize = size * nmemb;
    struct MemoryStruct *mem = (struct MemoryStruct *)data;
	
    mem->memory = (char *)myrealloc(mem->memory, mem->size + realsize + 1);
    if (mem->memory) {
		memcpy(&(mem->memory[mem->size]), ptr, realsize);
		mem->size += realsize;
		mem->memory[mem->size] = 0;
    }
    return realsize;
}

/*
 the fitfunction for the genetic optimisation
 */
int lgencurvefit_fitfunction(void *userdata, const double *coefs, unsigned int numcoefs, double *model, const double **xdata, long datapoints, unsigned int numDataDims){
	int err = 0;
	
	CountInt ii;
	unsigned int jj;
	int requiredParameterTypes[MAX_MDFIT_SIZE + 2];
	allFitFunc allParameters;
	fitFunc parameters;
	double result;
	double *dp;
	GenCurveFitInternals *goiP = (GenCurveFitInternals*) userdata;

	memset(&parameters, 0, sizeof(parameters));
	memset(&allParameters, 0, sizeof(allParameters));
	

	//cmd-dot or abort button
	if(CheckAbort(0) == -1){
		err = FIT_ABORTED;
		goto done;
	}

	/*
	 check if all the input isn't NULL
	 */
	if(goiP->coefs == NULL || goiP->dataCalc == NULL){
		err = UNSPECIFIED_ERROR;
		goto done;
	}
	
	/*
	 Place the parameters into the coefficient wave
	 */
	if(err = MDStoreDPDataInNumericWave(goiP->GenCurveFitCoefs, (double*) coefs))
		goto done;
			
	switch(goiP->isAAO){
		case 0:
			parameters.waveH = goiP->GenCurveFitCoefs;
			requiredParameterTypes[0] = WAVE_TYPE;
			
			for(ii = 0 ; ii < datapoints ; ii += 1){
				for(jj = 0 ; jj < numDataDims ; jj += 1)
					parameters.x[jj] = xdata[jj][ii];
				
				// call the users fit function and put the result in the output array
				if (err = CallFunction(&goiP->fi, (void*) &parameters, &result))
					goto done;
				model[ii] = result;
			}
						
			break;
		case 1:
			allParameters.waveC = goiP->GenCurveFitCoefs;
			allParameters.waveY = goiP->dataCalc;
			requiredParameterTypes[0] = WAVE_TYPE;
			requiredParameterTypes[1] = WAVE_TYPE;
			
			for(ii = 0 ; ii < numDataDims ; ii += 1){
				requiredParameterTypes[ii + 2] = WAVE_TYPE;
				if(goiP->xcalc[ii] == NULL){
					err = UNSPECIFIED_ERROR;
					goto done;
				}
				allParameters.waveX[ii] = goiP->xcalc[ii];
			}

			// call the users fit function and put the result in the output wave
			if (err = CallFunction(&goiP->fi, (void*)&allParameters, &result))
				goto done;

			// the user may have changed the number of points in the output wave
			if(goiP->dataCalc == NULL || WavePoints(goiP->dataCalc ) != datapoints){
				err = USER_CHANGED_FITWAVE;
				goto done;
			}
			dp = WaveData(goiP->dataCalc);
			memcpy(model, dp, sizeof(double) * datapoints);
			
			break;
		case 2:			
			if(goiP->sp == NULL){
				err = NULL_STRUCTURE;
				goto done;
			}
			goiP->sp->w = goiP->GenCurveFitCoefs;
			goiP->sp->yy = goiP->dataCalc;
			
			
			for(ii = 0 ; ii < numDataDims ; ii+=1){
				if(goiP->xcalc[ii] == NULL){
					err = UNSPECIFIED_ERROR;
					goto done;
				}
				goiP->sp->xx[ii] = goiP->xcalc[ii];
			}

			if (err = CallFunction(&goiP->fi, (fitfuncStruct*)&goiP->sp, &result))
				goto done;
			//don't want any dangling references to waves
			goiP->sp->yy = NULL;
			goiP->sp->w = NULL;
			memset(goiP->sp->xx, 0, sizeof(goiP->sp->xx));
			
			// the user may have changed the number of points in the output wave
			if(goiP->dataCalc == NULL || WavePoints(goiP->dataCalc) != datapoints){
				err = USER_CHANGED_FITWAVE;
				goto done;
			}
			
			dp = WaveData(goiP->dataCalc);
			memcpy(model, dp, sizeof(double) * datapoints);
			
			break;
		default:
			err = UNSPECIFIED_ERROR;
			goto done;
			break;
	}
	
	// check that the fitfunction didn't return any NaN or INF
	if(err = checkNanInfArray(model, datapoints)){
		err = FITFUNC_RETURNED_NANINF;
		goto done;
	}


done:		

	return err;
}

/*
 the cost function for the genetic optimisation
 */
double lgencurvefit_costfunction(void *userdata, const double *coefs, unsigned int numcoefs, const double *data, const double *model, const double *errors, long datapoints){
	double val;

	GenCurveFitInternals *goiP = (GenCurveFitInternals*) userdata;
	costFunc userCostFunc;
	double *dp;
	int err;
	
	switch(goiP->METH){//which cost function
		case 0:
			val = chisquared(userdata, coefs, numcoefs, data, model, errors, datapoints);			
			break;
		case 1:
			val = robust(userdata, coefs, numcoefs, data, model, errors, datapoints);
			break;
		case 2:
			userCostFunc.coefs = goiP->GenCurveFitCoefs;
			userCostFunc.yobs =	 goiP->yobs;
			userCostFunc.ycalc = goiP->dataCalc;
			userCostFunc.sobs =  goiP->sobs;
			
			//copy the original data into the yobs wave created for the purpose
			//some sneaky users probably try to change it.
			dp = WaveData(goiP->yobs);
			memcpy(dp, data, datapoints * sizeof(double));
			
			//copy the original data into the yobs wave created for the purpose
			//some sneaky users probably try to change it.
			dp = WaveData(goiP->sobs);
			memcpy(dp, errors, datapoints * sizeof(double));

			if(err = CallFunction(&goiP->minf, (void*) &userCostFunc, &val))
				goto done;
			
			if(WavePoints(goiP->yobs) != datapoints || WavePoints(goiP->sobs) != datapoints || WavePoints(goiP->GenCurveFitCoefs)!= numcoefs){
				err = COSTFUNC_WAVES_CHANGED;
				goto done;
			}
			
			if(IsNaN64(&val) || IsINF64(&val)){
				err = COSTFUNC_DOESNT_RETURN_NUMBER;
				goto done;
			}
			
			break;
		default:
			val = chisquared(userdata, coefs, numcoefs, data, model, errors, datapoints);			
			break;
	}
	done:
			
	return val;
};


/*
 the update function for the genetic optimisation
 */
int lgencurvefit_updatefunction(void *userdata,
								const double *coefs,
								 unsigned int numcoefs,
								 unsigned int iterations,
								 double cost,
								 unsigned int updatetime,
								 double convergenceNumber,
								 const double** population,
								 const unsigned int *varparams,
								 unsigned int numvarparams,
								 unsigned int totalpopsize,
								 const double *costmap){
	int err = 0;
	double val = 0;
	
	GenCurveFitInternals *goiP = (GenCurveFitInternals*) userdata;
	updtFunc userUpdateFunc;
	
	goiP->V_numfititers = iterations;
	
	/*
	 Display the coefficients so far.
	 */
	if(!goiP->noupdate && (updatetime == 1 || updatetime == 8))
		DisplayWindowXOP1Message(gTheWindow, numcoefs, coefs, cost, goiP->fi.name, iterations, updatetime, convergenceNumber);

	// perhaps the user wants to abort the fit using gui panel?
	if(Abort_the_fit){
		err = FIT_ABORTED;
		goto done;
	}

	//We want a record of all the improvements in the fit
	if(goiP->dump && (updatetime == 2 || updatetime == 1)){
		WriteMemoryCallback((void*) coefs, sizeof(double), goiP->totalnumparams, &(goiP->dumpRecord));
		if(goiP->dumpRecord.memory == NULL){
			err = NOMEM; goto done;
		}
	}
	
	//store the best fit coefficients in the original coef wave and update the output data.
	if((updatetime == 1 && !goiP->noupdate) || updatetime == 16){
		goiP->cost = cost;

		if(err = MDStoreDPDataInNumericWave(goiP->OUT_coefs, (void*) coefs))
			return err;
		WaveHandleModified(goiP->OUT_coefs);
		
		if(err = calcModelXY(&goiP->fi, goiP->OUT_coefs, goiP->OUT_data, goiP->OUT_x, goiP->numVarMD, goiP->isAAO, goiP->sp))
			return err;
		
		WaveHandleModified(goiP->OUT_data);
		
		if(goiP->OUT_res){
			if(err = calcModelXY(&goiP->fi, goiP->OUT_coefs, goiP->OUT_res, goiP->fullExtentOfData, goiP->numVarMD, goiP->isAAO, goiP->sp))
				return err;
			if(err = subtractTwoWaves(goiP->OUT_res, goiP->originalYobs))
				return err;
			if(err = scalarMultiply(goiP->OUT_res, -1))
				return err;
			WaveHandleModified(goiP->OUT_res);
		}
		DoUpdate();		
	}
	
	//the user defined a user IGOR update function so call that.
	if(!goiP->noupdate && goiP->useIgorUpdateFunction && population && costmap){
		userUpdateFunc.currentbestfit = goiP->OUT_coefs;
		userUpdateFunc.population = goiP->M_population;
		userUpdateFunc.costmap = goiP->W_costmap;
		userUpdateFunc.updatetime = updatetime;

		memcpy(WaveData(goiP->M_population), &(population[0][0]), totalpopsize * numvarparams * sizeof(double));
		memcpy(WaveData(goiP->W_costmap), costmap, totalpopsize * sizeof(double) );

		WaveHandleModified(goiP->M_population);
		WaveHandleModified(goiP->W_costmap);

		if(err = CallFunction(&goiP->igorUpdateFunction, (void*) &userUpdateFunc, &val))
			goto done;
		if(val)
			err = (int)val;
	}
	
	
done:
	return err;
}

/* 
 function obtains the covariance matrix for the fit.
 */
int getGCovarianceMatrix(GenCurveFitRuntimeParamsPtr p, GenCurveFitInternalsPtr goiP){
	int err = 0;
	double hessianDeterminant = 0;
	
	if(err = getCovarianceMatrix(&goiP->covarianceMatrix,
								 &hessianDeterminant,
								 goiP,
								 lgencurvefit_fitfunction,
								 goiP->cost,
								 goiP->coefs,
								 (int) goiP->totalnumparams,
								 goiP->holdvector,
								 goiP->dataObs,
								 goiP->dataSig,
								 (const double**) goiP->independentVariable,
								 (long) goiP->unMaskedPoints,
								 (int) goiP->numVarMD,
								 !p->WFlagEncountered))
		return err;
				
	goiP->V_logBayes = exp(-0.5 * (goiP->cost) / (double)(goiP->unMaskedPoints - goiP->numvarparams));// * pow(4*3.14159,(double) goiP->numvarparams) ;//* factorial((double)goiP->numvarparams);
	goiP->V_logBayes = goiP->V_logBayes / (sqrt(hessianDeterminant));
	goiP->V_logBayes = log(goiP->V_logBayes);

	return err;
}

/*
 converts a libgencurvefit error to an XOP error
 All libgencurvefit errors are negative, all XOP errors are positive
 */

int convertlgencurvefitErrorToXOP(int lgencurvefiterror){
	int err = lgencurvefiterror;
	
	if(lgencurvefiterror > -1)
		return err;
	
	switch (lgencurvefiterror) {
		case NO_MEMORY:
			err = NOMEM;
			break;
		case INCORRECT_LIMITS:
			err = LIMITS_INVALID;
			break;
		case HOLDVECTOR_COEFS_MISMATCH:
			err = HOLDSTRING_NOT_RIGHT_SIZE;
			break;
		case NO_VARYING_PARAMS:
			err = ALL_COEFS_BEING_HELD;
			break;
		case WRONG_NUMBER_OF_PARAMS:
			break;
		case COEFS_MUST_BE_WITHIN_LIMITS:
			err = LIMITS_INVALID;
			break;
		case PROBLEM_CALCULATING_COVARIANCE:
			err = SINGULAR_MATRIX;
			break;
		case NO_FIT_FUNCTION_SPECIFIED:
			err = FITFUNC_NOT_SPECIFIED;
			break;
		case NO_Y_ARRAY:
			break;
		case NO_X_ARRAY:
			break;
		case NO_E_ARRAY:
			break;
		case NO_COEFS_ARRAY:
			break;
		case NO_LIMITS_ARRAY:
			err = LIMITS_INVALID;
			break;
		case SINGULAR_MATRIX_ERROR:
			err = SINGULAR_MATRIX;
			break;
		default:
			break;
	}
	
	
	return err;
}


/*
 ExecuteGenCurveFit performs the genetic curvefitting routines
 returns 0 if no error
 returns errorcode otherwise
 */
int
ExecuteGenCurveFit(GenCurveFitRuntimeParamsPtr p)
{
	/* the function that is called by IGOR. Here's where we farm the work out to different places.*/
	/*
	 GenCurveFitInternals carries the internal data structures for doing the fitting.
	 */
	GenCurveFitInternals goi;
	//initialise the structure to be zero
	
	/*
	 err carries the errors for all the operations
	 err2 carries the error code if we still want to return err, but we need to finish off
	 something else first.	
	 */
	int err = 0, err2 = 0;
	double value[2];
	CountInt indices[MAX_DIMENSIONS];
	
	//libgencurvefit options
	gencurvefitOptions gco;
	
	//variables listed below are purely for outputs sake.
	char varname[MAX_OBJ_NAME + 1];
	CountInt dimensionSizes[MAX_DIMENSIONS + 1];
	double t1,t2, chi2;
	double *wP;
	long lt1 = 0;
	char note[200], note_buffer1[MAX_WAVE_NAME + 1], note_buffer2[MAX_WAVE_NAME + 1], cmd[MAXCMDLEN + 1];
	int output, ii, jj, isDisplayed, quiet;
	
	//initialise all the internal data structures to NULL
	memset(&goi, 0, sizeof(goi));
	
	//you have to use IGOR > 6.10
	if( igorVersion < 610 )
		return REQUIRES_IGOR_610;
	
	//reset the options for libgencurvefit
	memset(&gco, 0, sizeof(gencurvefitOptions));
	
	//reset the abort condition
	Abort_the_fit = 0;
	
	//can't do this in a threadsafe manner if we're not running IGOR >6.2
	if (igorVersion < 620 && !RunningInMainThread())
		return IGOR_OBSOLETE;
	
	strncpy(varname, "V_Fiterror", MAX_OBJ_NAME);
	if(FetchNumVar(varname, &t1, &t2)!=-1){
		if(!err)
			lt1 = 0;
		else 
			lt1 = 1;
		
		err = 0;
		if(err2 = SetIgorIntVar(varname, lt1, 1)){err = err2;goto done;};
	}	
	
	/*
	 Genetic Optimisation uses a lot of random numbers, we are seeding the generator here.
	 If you want to seed the generator we can 
	 */
	if(p->SEEDFlagEncountered)
		gco.seed = (int) p->SEEDFlag_seed;
	else 
		gco.seed = -1;
	
	/*
	 checkInput checks the input that IGOR sends to the XOP.  If everything is correct it returns 0, 
	 else it returns an error.  Errors can be caused by, e.g. different x and y wave lengths, etc.
	 */
	if(err = checkInput(p, &goi))
		goto done;
	
	/*
	 init_GenCurveFitInternals sets up the internal data arrays in the GenCurveFitInternals structure.  It holds
	 copies of the datapoints being fitted, arrays of the fitting structures, etc.  It calls malloc
	 several times to set aside memory for all this.  If this procedure works without a hitch then the function
	 returns 0.
	 */
	if(err = init_GenCurveFitInternals(p, &goi))
		goto done;
	
	if(p->OPTFlagEncountered && (((long)p->OPTFlag_opt) & (long)pow(2, 0)))
		gco.useinitialguesses = 1;
	
	
	gco.updatefun = (updatefunction) &lgencurvefit_updatefunction;
	gco.k_m = goi.k_m;
	gco.recomb = goi.recomb;
	gco.strategy = goi.STGY;
	gco.tolerance = goi.tolerance;
	gco.popsizeMultiplier = goi.popsize;
	gco.iterations = goi.iterations;
	gco.temp = goi.temperature;
	gco.updatefrequency = 31;
	
	/*
	 you want to polish the fit (at the end) using LevenbergMarquardt
	 */
	if(p->POLFlagEncountered)
		gco.polishUsingLM = 1;
	
	/*
	 if you want to do a Monte Carlo fit.
	 */
	if(p->MCFlagEncountered && p->WFlagEncountered)
		gco.monteCarlo = 1;
	
	//make an error wave before we start the fit, to make sure it's always there
	dimensionSizes[0] = goi.totalnumparams;
	dimensionSizes[1] = 0;
	if(err = MDMakeWave(&goi.W_sigma, "W_sigma", goi.cDF, dimensionSizes, NT_FP64, 1))
			goto done;
	for(ii = 0 ; ii < goi.totalnumparams ; ii+=1){
		indices[0] = ii;
		value[0] = 0;
		if(err = MDSetNumericWavePointValue(goi.W_sigma, indices, value))
			goto done;
	}
	WaveHandleModified(goi.W_sigma);
		
	if(p->MATFlagEncountered && !(p->MATFlagParamsSet[0] && (int) p->MATFlag_mat == 0)){
		//make the covariance matrix
		dimensionSizes[0] = goi.totalnumparams;
		dimensionSizes[1] = goi.totalnumparams;
		dimensionSizes[2] = 0;
		if(err = MDMakeWave(&goi.M_covariance, "M_Covar", goi.cDF, dimensionSizes, NT_FP64, 1))
				goto done;
			
		wP = WaveData(goi.M_covariance);
			
		memset(wP, 0, sizeof(double) * goi.totalnumparams * goi.totalnumparams);
			
		WaveHandleModified(goi.M_covariance);
	}
	
	/*
	 optimiseloop does the Differential Evolution, according to <NAME>.  When this returns 0, then 
	 you have the best fit in the GenCurveFitInternals structure.  Otherwise it returns an error code.  If the user aborts
	 then the FIT_ABORTED error code is returned, but it is still possible to retrieve the best solution so far
	 */
	if(err = genetic_optimisation(lgencurvefit_fitfunction,
							   lgencurvefit_costfunction,
							   (unsigned int) goi.totalnumparams,
							   goi.coefs,
							   goi.holdvector,
							   (const double**) goi.limits,
							   (long) goi.unMaskedPoints,
							   goi.dataObs,
							   (const double**) goi.independentVariable,
							   goi.dataSig,
							   goi.numVarMD,
							   &chi2,
							   &gco,
							   &goi))
		goto done;
	
	//make sure the coefficients for the fit are updated, by calling the update function
	Abort_the_fit = 0;
	if(err = lgencurvefit_updatefunction(&goi,
										 goi.coefs,
										 (unsigned int) goi.totalnumparams,
										 goi.V_numfititers,
										 chi2,
										 16,
										 1,
										 NULL,
										 NULL,
										 0,
										 0,
										 NULL))
	   goto done;
	
	//output the dumprecord
	if(!err && goi.dump && goi.dumpRecord.memory)
		if(err2 = dumpRecordToWave(&goi, &goi.dumpRecord)){
			err = err2;
			goto done;
		}
	
	/*
	 if there are no errors, or if the user aborted, then return the best fit.
	 If the data is displayed in the top graph append the fitcurve to the top graph
	 */
	if((err == 0 || err == FIT_ABORTED) && !getGCovarianceMatrix(p, &goi) && goi.covarianceMatrix){
		//set the error wave
		jj = 0;
		for(ii = 0; ii < goi.totalnumparams ; ii += 1){
			indices[0] = ii;
			if(!goi.holdvector[ii]){
				value[0] = sqrt(goi.covarianceMatrix[goi.varParams[jj]][goi.varParams[jj]]);
				jj++;
			} else {
				value[0] = 0;
			}

			if(err = MDSetNumericWavePointValue(goi.W_sigma, indices, value))
				goto done;
			
		}
		WaveHandleModified(goi.W_sigma);
		
		if(p->MATFlagEncountered && !(p->MATFlagParamsSet[0] && (int) p->MATFlag_mat == 0)){
			//make the covariance matrix
			dimensionSizes[0] = goi.totalnumparams;
			dimensionSizes[1] = goi.totalnumparams;
			dimensionSizes[2] = 0;
			wP = WaveData(goi.M_covariance);			
			memcpy(wP, *goi.covarianceMatrix, sizeof(double) * goi.totalnumparams * goi.totalnumparams);
			WaveHandleModified(goi.M_covariance);
		}
	}
	//append the fit to the graph
	if(!err || err == FIT_ABORTED){
		if(err = isWaveDisplayed(p->dataWave.waveH, &isDisplayed))
			goto done;
		
		if(isDisplayed && goi.numVarMD == 1){
			if(err = isWaveDisplayed(goi.OUT_data, &isDisplayed))
				goto done;
			
			if(!isDisplayed){
				strncpy(cmd, "appendtograph/w=$(winname(0,1)) ", MAXCMDLEN);
				WaveName(goi.OUT_data,&note_buffer1[0]);
				strncat(cmd, &note_buffer1[0], MAXCMDLEN - strlen(note_buffer1));
				
				if(p->DFlagEncountered && p->XFlagEncountered){
					WaveName(p->XFlag_xx,&note_buffer2[0]);
					strncat(cmd, " vs ", MAXCMDLEN - strlen(cmd) - strlen(" vs "));
					strncat(cmd, note_buffer2, MAXCMDLEN - strlen(cmd) - strlen(note_buffer2) );
				}
				if(err = XOPSilentCommand(&cmd[0]))
					goto done;
			}
		}
	}
	
	/*this section sets history and global variables
	 V_Chisq
	 V_Fiterror
	 If there are no errors returned above, then V_Chisq is produced in the current datafolder.
	 V_fiterror is for if there is a problem with the fitting.  
	 If this global variable exists in IGOR at runtime, and there is no XOPerror then V_fiterror = 0.  
	 If there is an XOPerror, bit 0 of V_fiterror is set. However, the XOP will return gracefully to IGOR,
	 allowing the user that called it to detect this and carry on.
	 In normal situations V_fiterror may not exist, so the GenCurveFit XOP will return an error message.
	 */
	
	if(!err){
		SetOperationNumVar("V_Chisq", chi2);
		SetOperationNumVar("V_fitIters", (int)(goi.V_numfititers));
		SetOperationNumVar("V_npnts", (int)(goi.unMaskedPoints));
		SetOperationNumVar("V_nterms", (int)WavePoints(p->coefs));
		SetOperationNumVar("V_nheld", (int)(WavePoints(p->coefs) - goi.numvarparams));
		SetOperationNumVar("V_logBayes", goi.V_logBayes);
	}		
	
	strcpy(varname, "V_Fiterror");
	if(FetchNumVar(varname, &t1, &t2)!=-1){
		if(!err){
			lt1 = 0;
		} else {
			lt1 = 1;
		}
		err = 0;
		if(err2 = SetIgorIntVar(varname, lt1, 1))
		{err = err2;goto done;}
	}
	/*
	 This section puts a copy of the fit parameters into the history area, unless one sets quiet mode.
	 */
	quiet = 0;
	if(p->QFlagEncountered){
		quiet = 1;
		if(p->QFlagParamsSet[0] && (int) p->QFlag_quiet == 0)
			quiet = 0;
	}
	
	if(!quiet && (!err || err == FIT_ABORTED) && lt1==0 ){		
		if(!err)
			{output = snprintf(note, 199, "_______________________________\rGenetic Optimisation Successful\r");XOPNotice(note);}
		if(err == FIT_ABORTED)
			{output = snprintf(note, 199, "_______________________________\rGenetic Optimisation ABORTED\r");XOPNotice(note);}
		WaveName(p->dataWave.waveH, note_buffer1);
		
		output = snprintf(note, 
						  199,
						  "Fitting: %s to %s\r",
						  note_buffer1,
						  goi.fi.name);
		
		XOPNotice(note);
		
		output = snprintf(note,
						  199,
						  "V_fitIters = %li; V_Chisq = %g; V_npnts= %lld; V_nterms= %lld; V_nheld= %lld; V_logBayes = %g\r",
						  goi.V_numfititers,
						  chi2,
						  (SInt64) goi.unMaskedPoints,
						  (SInt64) WavePoints(p->coefs),
						  (SInt64) WavePoints(p->coefs) - goi.numvarparams,
						  goi.V_logBayes);
		
		XOPNotice(note);
		
		for(ii = 0; ii < WavePoints(p->coefs) ; ii += 1){
			indices[0] = ii;
			indices[1] = 0;
			if(err = MDGetNumericWavePointValue(goi.W_sigma, indices, value))
				goto done;
			
			output = snprintf(note,
{'fraction_non_alphanumeric': 0.11954389965792474, 'fraction_numerical': 0.011516533637400229, 'mean_word_length': 3.4240314769975786, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 310, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5791186', 'n_tokens_mistral': 19000, 'n_tokens_neox': 15719, 'n_words': 4048}