/*------------------------------------------------------------------------------*
 * File Name: FitToolsGetXfromY.c 												*
 * Creation: ER, 01/09/2002														*
 * Purpose: Function to find X value corresponding to Y value for Fit Tools		*
 * Copyright (c) OriginLab Corp.2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007	*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *------------------------------------------------------------------------------*/
 
////////////////////////////////////////////////////////////////////////////////////
//
#include <origin.h>
//
////////////////////////////////////////////////////////////////////////////////////
//
// Prototype of functions in this file:
double PolyGetXfromY(double dYval, string strFitLine, double dA, double dB1, double dB2 = 0, 
			double dB3 = 0, double dB4 = 0, double dB5 = 0, double dB6 = 0, double dB7 = 0, 
			double dB8 = 0, double dB9 = 0);
double SigmoidGetXfromY(double dYval, string strFitLine, double dP1, double dP2, 
			double dP3, double dP4, int iSigType);
double GetXfromY(string strFitLine, double *dCoeff, int iFlag);
double ComputeYfromX(double dX, double *dCoeff, int iFlag);
//
////////////////////////////////////////////////////////////////////////////////////
//
//
//
//
//
////////////////////////////////////////////////////////////////////////////////////
// 
// PolyGetXfromY: Gets X value corresponding to Y value for polynomial fit tool
// Parameters:
//					dYval : Y value for which X value is desired
//					strFitLine: name of fit curve created by the tool
//					dA, dB1,...dB9: fit coefficients
// Return:
//					Returns X value
//					Returns missing value if:
//								Y value passed is outside of fit curve limits, or
//								bisectional search algorithm fails
//
////////////////////////////////////////////////////////////////////////////////////
double PolyGetXfromY(double dYval, string strFitLine, double dA, double dB1, double dB2, 
			double dB3, double dB4, double dB5, double dB6, double dB7, 
			double dB8, double dB9)
{
	// Stuff all the coefficients of the polynomial into an array
	double dCoeff[10];
	dCoeff[0] = dA;
	dCoeff[1] = dB1;	
	dCoeff[2] = dB2;
	dCoeff[3] = dB3;
	dCoeff[4] = dB4;
	dCoeff[5] = dB5;
	dCoeff[6] = dB6;
	dCoeff[7] = dB7;
	dCoeff[8] = dB8;
	dCoeff[9] = dB9;
	int iFlag = 0;
	return GetXfromY(dYval, strFitLine, dCoeff, iFlag);
}



////////////////////////////////////////////////////////////////////////////////////
// 
// SigmoidGetXfromY: Gets X value corresponding to Y value for sigmoidal fit tool
// Parameters:
//					dYval : Y value for which X value is desired
//					strFitLine: name of fit curve created by the tool
//					dP1...dP4: fit coefficients
//					iSigType: Type of sigmoidal function:
//							1: Boltzmann  2: Dose Resp   3: Logistic
// Return:
//					Returns X value
//					Returns missing value if:
//								Y value passed is outside of fit curve limits, or
//								bisectional search algorithm fails
//
////////////////////////////////////////////////////////////////////////////////////
double SigmoidGetXfromY(double dYval, string strFitLine, double dP1, double dP2, 
			double dP3, double dP4, int iSigType)
{
	// Stuff all the coefficients into an array
	double dCoeff[4];
	dCoeff[0] = dP1;
	dCoeff[1] = dP2;	
	dCoeff[2] = dP3;
	dCoeff[3] = dP4;
	return GetXfromY(dYval, strFitLine, dCoeff, iSigType);
}



////////////////////////////////////////////////////////////////////////////////////
// 
// GetXfromY: 		Gets X value corresponding to Y value for both tools
//					This function first locates the (first) pair of fit data points 
//					within which dYval lies, and then performs a bisectional search
//					within this interval to determine the corresponding x value
//					within a precision of 0.00001% or within 200 iterations 
// Parameters:
//					dYval : Y value for which X value is desired
//					strFitLine: name of fit curve created by the tool
//					dCoeff: Array holding fit coefficients
//					iFlag: determines type of fit tool and function
//							0: Polynomial
//							1: Sigmoidal (Boltzmann)
//							2: Sigmoidal (DoseREsp)
//							3: Sigmoidal (Logistic)
// Return:
//					Returns X value
//					Returns missing value if:
//								Y value passed is outside of fit curve limits, or
//								bisectional search algorithm fails
//
////////////////////////////////////////////////////////////////////////////////////
double GetXfromY(double dYval, string strFitLine, double *dCoeff, int iFlag)

{
	// Create a curve object of the fit line dataset
	Curve crvFitLine(strFitLine);
	// If curve does not exist, return missing value
	if(!crvFitLine.IsValid()) return NANUM;
	// Create x dataset from fit line
	Dataset dsXFitLine;
	crvFitLine.AttachX(dsXFitLine);
	// Also create a y dataset of the fit line
	Dataset dsYFitLine(strFitLine);	

	// If yval is outside of [ymin,ymax] of yfit dataset, return missing value
	if( (dYval < min(dsYFitLine)) || (dYval > max(dsYFitLine)) ) return NANUM;
	
	// Now check if the yval for which user seeks x, is an exact point in the fit dataset
	int iIsPoint = Data_list(dYval, &dsYFitLine);
	// If it is a point, just return the x corresponding to the index
	if(iIsPoint != -1) return dsXFitLine[iIsPoint];

		
	// Now find (first) x interval in which yval belongs
	int iSize = dsYFitLine.GetSize();
	for(int ik = 0; ik < (iSize - 1); ik++)
	{
		double dSign = (dsYFitLine[ik] - dYval) * (dsYFitLine[ik+1] - dYval);
		if (dSign < 0) break;
	}
	if (ik == (iSize - 1)) return NANUM;
	double dXLo = dsXFitLine[ik];
	double dXHi = dsXFitLine[ik+1];	
	
	// Now start iterative process in the interval [dXLo, dXHi]
	double dXMid, dYatXLo, dYatXHi, dYOld, dYNew, dYChange;
	dXMid = (dXLo + dXHi) / 2.0;	
	dYatXLo = ComputeYfromX(dXLo, dCoeff, iFlag);
	dYatXHi = ComputeYfromX(dXHi, dCoeff, iFlag);	 
	dYOld = ComputeYfromX(dXMid, dCoeff, iFlag);
		
	int ii = 0;
	do
	{
		ii++;
		if(dYatXLo > dYatXHi)
		{
			if(dYOld > dYval) dXLo = dXMid;
			else dXHi = dXMid;
			
		}
		else
		{
			if(dYOld < dYval) dXLo = dXMid;
			else dXHi = dXMid;
		}
		dYatXLo = ComputeYfromX(dXLo, dCoeff, iFlag);
		dYatXHi = ComputeYfromX(dXHi, dCoeff, iFlag);	 
		dXMid = (dXLo + dXHi) / 2.0;
		dYNew = ComputeYfromX(dXMid, dCoeff, iFlag);
		// Compute percentage change in y
		dYChange = 100 * fabs( (dYNew - dYval) / dYval);
		dYOld = dYNew;
		// Loop until % change in y is less than 0.0001 and loop has not been executed 100 times
	} while ( (dYChange > 0.00001) && (ii < 200) ); 		
	
	// If number of iterations equals 200, did not converge and so return missing value
	if(ii == 200) return NANUM;
	
	return dXMid;	// Return the current dXMid value - this is the x corresponding to dYval
}



////////////////////////////////////////////////////////////////////////////////////
// 
// ComputeYfromX: 	Computes Y value given an X value, using appropriate equation
// Parameters:
//					dX: X value for which Y value is to be computed
//					dCoeff: Array holding fit coefficients
//					iFlag: determines equation to be used
//							0: Polynomial
//							1: Sigmoidal (Boltzmann)
//							2: Sigmoidal (DoseREsp)
//							3: Sigmoidal (Logistic)
// Return:
//					Returns computed Y value
//
////////////////////////////////////////////////////////////////////////////////////
double ComputeYfromX(double dX, double *dCoeff, int iFlag)
{
	if(iFlag == 0)
	{
		double dYval = dCoeff[0];
		for(int ii = 1; ii <= 9; ii ++)
		{
			dYval += dCoeff[ii] * dX^ii;
		}
		return dYval;
	}
	else if(iFlag == 1) return nlfxBoltzman(dX,dCoeff[0],dCoeff[1],dCoeff[2],dCoeff[3]);
	else if(iFlag == 2) return nlfxResponse1(dX,dCoeff[0],dCoeff[1],dCoeff[2],dCoeff[3]);
	else return nlfxLogistic(dX,dCoeff[0],dCoeff[1],dCoeff[2],dCoeff[3]);
}


// end of all function in file



