/*------------------------------------------------------------------------------*
 * File Name: CompareFits.c														*
 * Creation: 2001.10.11															*
 * Purpose: OriginC Source C file												*
 * Copyright (c) OriginLab Corp. 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007	*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *	ML 2/6/2002 PROPER_CONCATENATION											*
 *------------------------------------------------------------------------------*/
 
#include <stdio.h> // Most usful to include this, has printf and etc.
//#include <wksheet.h> // worksheet matrixlayer and column classes
#include <data.h>	 // Dataset, vector, matrix etc.
#include <File.h>	// File and StdioFile classes for basic file access
#include <Page.h> //File containing declaration of the class Page
#include <utilities.h> // basic utilities, including Win32 Kernal functions through mswin.h
//#include <graph.h> // graph, plot and graphic objects
//#include <complex.h>
#include <math.h> // all math function prototypes and constants
#include <stat.h>

#include "App_Utils.h"
#include "$local.h"

// Prototypes
int CompareFits(string strCurve1, string strCurve2, string strResWks, string strCategory, string strFunc);
static void FitCmpOut(double dfval, double dpval, Dataset &AData, Dataset &BData, string strCurve1, string strCurve2, string strFunc);

// perform fits on two datasets with same function, and compare the fit results
int CompareFits(string strCurve1, string strCurve2, string strResWks, string strCategory, string strFunc)
{

	// these are the curve names for the two fit datasets	
	Curve curve1(strCurve1);
	Curve curve2(strCurve2);

	// define datasets to hold x,y colums for the combined dataset
	Dataset xxData(strResWks,0);
	Dataset yyData(strResWks,1);
	// and the fit results
	Dataset AResults(strResWks,2);
	Dataset BResults(strResWks,3);

	// get size of fit curves
	int isize1 = curve1.GetSize();
	int isize2 = curve2.GetSize();

	// set size of combined datasets
	xxData.SetSize(isize1+isize2);
	yyData.SetSize(isize1+isize2);
	string strCombined;
	yyData.GetName(strCombined);

	// now combine the two y datasets of the curves
	if(curve1.HasX())
	{
		Dataset sx1;
		curve1.AttachX(sx1);
    	for (int ii=0; ii < isize1; ii++)
		{
			xxData[ii] = sx1[ii];
		}
	}
	else
	{
		for (int ii=0; ii < isize1; ii++)
		{
			xxData[ii] = ii+1;
		}

	}
	for (int ii=0; ii < isize1; ii++)
	{
		yyData[ii] = curve1[ii];
	}
	
	if(curve2.HasX())
	{
		Dataset sx2;
		curve2.AttachX(sx2);
    	for (int ii=0; ii < isize2; ii++)
		{
			xxData[ii+isize1] = sx2[ii];
		}
	}
	else
	{
		for (int ii=0; ii < isize2; ii++)
		{
			xxData[ii+isize1] = ii+1;
		}
	}
	for (ii=0; ii < isize2; ii++)
	{
		yyData[ii+isize1] = curve2[ii];
	}

	// now start fitting - first we fit the two data sets and the the combined dataset
	double ssr_combined, dof_combined;
	int jj, FitCmpNPara, iLineFunc;
	char pname[32];
		
	_LT_Obj
	{
		nlsf.category(strCategory);
		nlsf.init();
		nlsf.func$=strFunc;		
		nlsf.output(0);	
		nlsf.fitdata$=strCurve1;
//		nlsf.output(0);	
		nlsf.execute("parainit");
		nlsf.iterateex(1000);
		nlsf.output();
		FitCmpNPara=nlsf.npara;
		AResults.SetSize(4 + 3 * FitCmpNPara);
		AResults[0]=nlsf.dof;						// DOF
		AResults[1]=nlsf.ssr;						// SSR
		AResults[2]=nlsf.chisqr;					// CHISQR
		AResults[3]=FitCmpNPara;					// NPARA
		for(ii=1, jj=4; ii <= FitCmpNPara; ii++)
		{
			LT_set_var("ii", ii);
			LT_execute("%L = nlsf.n$(ii)$;");
			LT_get_str("%L", pname, 32);
			AResults.SetText(jj,pname);				// PARAMETER NAME
			jj++;
			AResults[jj] = nlsf.p$(ii);				// PARAMETER VALUE
			jj++;
			AResults[jj] = nlsf.e$(ii);				// PARAMETER ERROR
			jj++;
		}
		
		nlsf.init();
		nlsf.func$=strFunc;	
		nlsf.output(0);	
		nlsf.fitdata$=strCurve2;	
//		nlsf.output(0);	
		nlsf.execute("parainit");
		nlsf.iterateex(1000);
		nlsf.output();	
		BResults.SetSize(4 + 3 * FitCmpNPara);
		BResults[0]=nlsf.dof;
		BResults[1]=nlsf.ssr;
		BResults[2]=nlsf.chisqr;
		BResults[3]=FitCmpNPara;
		for(ii=1, jj=4; ii <= FitCmpNPara; ii++)
		{
			LT_set_var("ii", ii);
			LT_execute("%L = nlsf.n$(ii)$;");
			LT_get_str("%L", pname, 32);
			BResults.SetText(jj,pname);
			jj++;
			BResults[jj] = nlsf.p$(ii);
			jj++;
			BResults[jj] = nlsf.e$(ii);
			jj++;
		}
		
		nlsf.init();
		nlsf.func$=strFunc;	
		nlsf.output(0);	
		nlsf.fitdata$=strCombined;	
//		nlsf.output(0);	
		nlsf.execute("parainit");
		nlsf.iterateex(1000);
		nlsf.output();	
		dof_combined = nlsf.dof;
		ssr_combined = nlsf.ssr;
	}	

	// compute f statistics
	double ssr_sep = AResults[1] + BResults[1];
	double dof_sep = AResults[0] + BResults[0];
	double f1 = (ssr_combined - ssr_sep) / (dof_combined - dof_sep );
	double f2 = ssr_sep;
	double f  =f1 /f2;
	if(f < 0) f = 0 - f;
	if(f < 1E-8) f = 0.0;
	double p = 1 - invf(f, (dof_combined - dof_sep), dof_sep);
	
	FitCmpOut(f, p, AResults, BResults, strCurve1, strCurve2, strFunc);

	return 0;
}


static void FitCmpOut(double dfval, double dpval, Dataset &AData, Dataset &BData, string strCurve1, string strCurve2, string strFunc)
{
	Worksheet		wks;
	string			sWksName, sStr, sStr1, sStr2;
	unsigned int	iRow = 0;
	unsigned int	iCol = 0;
	unsigned int	iLoop = 0;
	unsigned int	iIdx = 0;
	int				iTemp;
	double			dTemp, dTemp1, dTemp2, dTemp3, dTemp4, dFitCmpCL;
	BOOL			bRes;
	
	wks.Create("fitcmp1", CREATE_TEMP);
	sWksName = wks.GetPage().GetName();
	LT_set_str("%M", sWksName);
	LT_get_var("FitCmpCL", &dFitCmpCL);
	wks.AppendRows(35 + 2 * AData[3]);	// include space for Line function test
	LT_execute("type.NumCols = 6;type.BeginResults();");

// The home position
	sStr1.Format(FC_HEADER1, strCurve1, strCurve2, strFunc);
	wks.SetCell(iRow, iCol, sStr1);
	iRow++;
// Blank everything in this row
	for(iCol = 0;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, " ");
	iRow++;

// First Fit
// Indent and type a Fit results message
	sStr1.Format(FC_HEADER4, strCurve1);
	wks.SetCell(iRow, 0, " ");
	wks.SetCell(iRow, 1, sStr1);
	for(iCol=2 ; iCol <= 5; iCol++)
		wks.SetCell(iRow,iCol," ");
	iRow++;
// Ruler
	wks.SetCell(iRow, 0, " ");
	for(iCol = 1;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, "---");
	iRow++;
// Headers
	wks.SetCell(iRow,0," ");
	sStr1 = FC_HEADER5;
	for(iCol = 1; iCol <= 5; iCol++)
	{
		sStr2 = sStr1.GetToken(iCol - 1, ' ');
		wks.SetCell(iRow, iCol, sStr2);
	}
	iRow++;
// Parameters
	for(iLoop = 1,iIdx = 4;iLoop <= AData[3];iLoop++,iRow++)
	{
		wks.SetCell(iRow,0," ");
		bRes = AData.GetText(iIdx, sStr1);
		iIdx++;
		wks.SetCell(iRow, 1, sStr1);
		dTemp1 = AData[iIdx];
		iIdx++;
		wks.SetCell(iRow, 2, dTemp1);
		dTemp2 = AData[iIdx];
		iIdx++;
		wks.SetCell(iRow, 3, dTemp2);
		dTemp = dTemp2 * tTable(dFitCmpCL, AData[0]);
		wks.SetCell(iRow, 4, dTemp1 - dTemp);
		wks.SetCell(iRow, 5, dTemp1 + dTemp);
	}

// Blank row
	for(iCol = 0;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, " ");
	iRow++;

// Special quantities
	wks.SetCell(iRow, 0, " ");
	sStr1 = "SSR = ";
	if(AData[1] == NANUM) sStr2 = "--"; else sStr2 = LocalizeDouble(AData[1], "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;
	wks.SetCell(iRow, 0, " ");
	sStr1 = "DOF = ";
	if(AData[0] == NANUM) sStr2 = "--"; else sStr2 = LocalizeDouble(AData[0], "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;
	wks.SetCell(iRow, 0, " ");
	sStr1 = "ChiSqr / DOF = ";
	if(AData[2] == NANUM) sStr2 = "--"; else sStr2 = LocalizeDouble(AData[2], "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
//	sStr1.Format("ChiSqr / DOF = %f", AData[2] / AData[0]);
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;

// Ruler
	wks.SetCell(iRow, 0, " ");
	for(iCol = 1;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, "---");
	iRow++;

// Blank row
	for(iCol = 0;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, " ");
	iRow++;

// Second Fit
// Indent and type a Fit results message
	sStr1.Format(FC_HEADER4, strCurve2);
	wks.SetCell(iRow, 0, " ");
	wks.SetCell(iRow, 1, sStr1);
	for(iCol=2 ; iCol <= 5; iCol++)
		wks.SetCell(iRow,iCol," ");
	iRow++;

// Ruler
	wks.SetCell(iRow, 0, " ");
	for(iCol = 1;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, "---");
	iRow++;

// Headers
	wks.SetCell(iRow,0," ");
	sStr1 = FC_HEADER5;
	for(iCol = 1; iCol <= 5; iCol++)
	{
		sStr2 = sStr1.GetToken(iCol - 1, ' ');
		wks.SetCell(iRow, iCol, sStr2);
	}
	iRow++;
// Parameters
	for(iLoop = 1,iIdx = 4;iLoop <= (int) BData[3];iLoop++,iRow++)
	{
		wks.SetCell(iRow,0," ");
		bRes = BData.GetText(iIdx, sStr1);
		iIdx++;
		wks.SetCell(iRow, 1, sStr1);
		dTemp1 = BData[iIdx];
		iIdx++;
		wks.SetCell(iRow, 2, dTemp1);
		dTemp2 = BData[iIdx];
		iIdx++;
		wks.SetCell(iRow, 3, dTemp2);
		dTemp = dTemp2 * tTable(dFitCmpCL, BData[0]);
		wks.SetCell(iRow, 4, dTemp1 - dTemp);
		wks.SetCell(iRow, 5, dTemp1 + dTemp);
	}

// Blank row
	for(iCol = 0;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, " ");
	iRow++;

// Special quantities
	wks.SetCell(iRow, 0, " ");
	sStr1 = "SSR = ";
	if(BData[1] == NANUM) sStr2 = "--"; else sStr2 = LocalizeDouble(BData[1], "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
//	sStr1.Format("SSR = %f", BData[1]);
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;
	wks.SetCell(iRow, 0, " ");
	sStr1 = "DOF = ";
	if(BData[0] == NANUM) sStr2 = "--"; else sStr2 = LocalizeDouble(BData[0], "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
//	sStr1.Format("DOF = %f", BData[0]);
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;
	wks.SetCell(iRow, 0, " ");
	sStr1 = "ChiSqr / DOF = ";
	if(BData[2] == NANUM) sStr2 = "--"; else sStr2 = LocalizeDouble(BData[2], "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
//	sStr1.Format("ChiSqr / DOF = %f", BData[2] / BData[0]);
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;

// Ruler
	wks.SetCell(iRow, 0, " ");
	for(iCol = 1;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, "---");
	iRow++;

// Blank row
	for(iCol = 0;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, " ");
	iRow++;

// F-Test Results
	wks.SetCell(iRow, 0, " ");
	wks.SetCell(iRow, 1, FC_HEADER6);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;

// Ruler
	wks.SetCell(iRow, 0, " ");
	for(iCol = 1;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, "---");
	iRow++;

// F and P Values
	wks.SetCell(iRow, 0, " ");
	sStr1 = FC_HEADER5B;
	if(dfval == NANUM) sStr2 = "NAN"; else sStr2 = LocalizeDouble(dfval, "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
	sStr2 = ", P = ";
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
	if(dpval == NANUM) sStr2 = "NAN"; else sStr2 = LocalizeDouble(dpval, "%G");
	/// ML 2/6/2002 PROPER_CONCATENATION
	//lstrcat(sStr1,sStr2);
	sStr1 += sStr2;
	/// end PROPER_CONCATENATION
//	sStr1.Format("F-Value = %f, P = %f", dfval, dpval);
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;

// Ruler
	wks.SetCell(iRow, 0, " ");
	for(iCol = 1;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, "---");
	iRow++;

// Conclusion
	LT_get_var("OFitCmpCL", &dFitCmpCL);
	dTemp = 1.0 - dFitCmpCL / 100.0;
	if(dpval > dTemp)
	{
		sStr1.Format(FC_FOOTER1, dTemp);
	}
	else
	{
		sStr1.Format(FC_FOOTER2, dTemp);
	}
	wks.SetCell(iRow, 0, " ");
	wks.SetCell(iRow, 1, sStr1);
	for(iLoop = 2; iLoop <= 5; iLoop++)
	{
		wks.SetCell(iRow, iLoop, " ");
	}
	iRow++;

// Blank row
	for(iCol = 0;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, " ");
	iRow++;

// Line Fit comparison
	if(strFunc.CompareNoCase("Line") == 0)
	{
		wks.SetCell(iRow, 0, " ");
		wks.SetCell(iRow, 1, FC_LINEARHEADER1);
		for(iCol = 2;iCol <= 5;iCol++)
			wks.SetCell(iRow, iCol, " ");
		iRow++;
		// Ruler
		wks.SetCell(iRow, 0, " ");
		for(iCol = 1;iCol <= 5;iCol++)
			wks.SetCell(iRow, iCol, "---");
		iRow++;
		// Intercepts
		dTemp1 = AData[6] * AData[6] / AData[2];
		dTemp2 = BData[6] * BData[6] / BData[2];
		dTemp3 = (AData[0] * AData[2] + BData[0] * BData[2]) / (AData[0] + BData[0]);
		LT_get_var("FitCmpCL", &dFitCmpCL);
		dTemp = tTable(dFitCmpCL, AData[0] + BData[0]) * sqrt(dTemp3) * sqrt(dTemp1 + dTemp2);
		dTemp1 = AData[5] - BData[5] - dTemp;
		dTemp2 = AData[5] - BData[5] + dTemp;
		dTemp3 = AData[8] - BData[8] - dTemp;
		dTemp4 = AData[8] - BData[8] + dTemp;
		sStr1.Format(FC_LINEARHEADER2, dTemp1, dTemp2);
		sStr2.Format(" %f", NANUM); // GRD - Replace missing value with dashes
		sStr1.Replace(sStr2," --");
		sStr1.Replace(sStr2," --");
		wks.SetCell(iRow, 0, " ");
		wks.SetCell(iRow, 1, sStr1);
		for(iCol = 2;iCol <= 5;iCol++)
			wks.SetCell(iRow, iCol, " ");
		iRow++;
		wks.SetCell(iRow, 0, " ");
		if(dTemp1 > dTemp2)
		{
			dTemp = dTemp1;
			dTemp1 = dTemp2;
			dTemp2 = dTemp;
		}
		LT_get_var("OFitCmpCL", &dFitCmpCL);
		dFitCmpCL = 1.0 - dFitCmpCL / 100.0;
		// GRD - Handle missing values
		if((dTemp1 <= 0 && dTemp2 >= 0) || (dTemp1 == NANUM || dTemp2 == NANUM))
		{
			sStr1.Format(FC_LINEARFOOTER1A, dFitCmpCL);
		}
		else
		{
			sStr1.Format(FC_LINEARFOOTER1B, dFitCmpCL);
		}
		wks.SetCell(iRow, 1, sStr1);
		for(iCol = 2;iCol <= 5;iCol++)
			wks.SetCell(iRow, iCol, " ");
		iRow++;

		// Blank row
		for(iCol = 0;iCol <= 5;iCol++)
			wks.SetCell(iRow, iCol, " ");
		iRow++;

		// Slopes
		sStr1.Format(FC_LINEARHEADER3, dTemp3, dTemp4);
		sStr2.Format(" %f", NANUM); // GRD - Replace missing value with dashes
		sStr1.Replace(sStr2," --");
		sStr1.Replace(sStr2," --");
		wks.SetCell(iRow, 0, " ");
		wks.SetCell(iRow, 1, sStr1);
		for(iCol = 2;iCol <= 5;iCol++)
			wks.SetCell(iRow, iCol, " ");
		iRow++;
		wks.SetCell(iRow, 0, " ");
		if(dTemp3 > dTemp4)
		{
			dTemp = dTemp2;
			dTemp3 = dTemp4;
			dTemp4 = dTemp;
		}
		// GRD - Handle missing values
		if((dTemp3 <= 0 && dTemp4 >= 0) || (dTemp3 == NANUM || dTemp4 == NANUM))
		{
			sStr1.Format(FC_LINEARFOOTER2A, dFitCmpCL);
		}
		else
		{
			sStr1.Format(FC_LINEARFOOTER2B, dFitCmpCL);
		}
		wks.SetCell(iRow, 1, sStr1);
		for(iCol = 2;iCol <= 5;iCol++)
			wks.SetCell(iRow, iCol, " ");
		iRow++;

	}

// Blank row
	for(iCol = 0;iCol <= 5;iCol++)
		wks.SetCell(iRow, iCol, " ");
	iRow++;

	LT_execute("type.wks(%M);type;type.EndResults();");
}


