/*------------------------------------------------------------------------------*
 * File Name:Analysis_utils.c 													*
 * Creation: CPY 3/11/03														*
 * Purpose: Origin's basic internal analysis routines							*
 * Copyright (c) Originlab Corp.	2003										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *------------------------------------------------------------------------------*/
 
#include <origin.h> // main Origin C header that is precompiled and already include most headers 

////////////////////////////////////////////////////////////////////////////////////
#define _DBINT(_STR, _INT)	//out_int(_STR, _INT);
#define _DBMGS(_STR)		//out_str(_STR);

#define MAX_MATRIX_POINTS 10000000

bool set_curve_input(const curvebase& cuvInput, TreeNode& trInput, int& i0, int& imax, HWND& hWndRet, bool bNeedInit)// = true);
{
	bool bFullRange = cuvInput.GetSourceRange(i0, imax)? false:true;
	if(bNeedInit) // not supplied, we need to init trInput
	{
		set_active_layer(trInput);
		set_curve(cuvInput, trInput);
	}
	else
	{	
		if(!trInput.Range1.UseRange.IsValid())
			trInput.Range1.UseRange.nVal = 0;
		
		bFullRange = trInput.Range1.UseRange.nVal > 0? false:true;
		// we need to set cuvInput active
		GraphLayer gl = get_graph_layer(trInput);
		if(gl)
		{
			set_active_layer(gl);
			Page pg = gl.GetPage();
			if(pg)
				hWndRet = pg.GetWindow().GetSafeHwnd();
			
			string str;
			str.Format("set %s -a",cuvInput.GetName());
			gl.LT_execute(str);
		}
	}
	return bFullRange;
}

//	add to TreeNode the X and Y dataset names of the given curve
bool set_curve(const curvebase& cuvInput, TreeNode& trNode)
{
	if(!cuvInput)
		return false;
	if(!trNode)
		return false;

	int	i0 = 0, imax = 0;
	trNode.Range1.UseRange.nVal = cuvInput.GetSourceRange(i0, imax) ? 1 : 0;
	trNode.Range1.R1.nVal = cuvInput.GetLowerBound();
	trNode.Range1.R2.nVal = cuvInput.GetUpperBound();

	trNode.Range1.Ydata.strVal = cuvInput.GetName();
	Dataset dsX;
	if(cuvInput.AttachX(dsX))
		trNode.Range1.Xdata.strVal = dsX.GetName();
	else
		trNode.Range1.Xdata.strVal = "";
	return true;
}

// add to TreeNode the name of the active graph page and the active layer number
bool set_active_layer(TreeNode& trNode)
{
	GraphPage gpg = Project.Pages();
	if(gpg)
	{
		GraphLayer glyr = Project.ActiveLayer();
		trNode.Page.strVal = gpg.GetName();
		trNode.Layer.nVal = glyr.GetIndex();
		return true;
	}
	return false;
}

//set the given layer to be active layer, if page not open, will open it to be active page as well
bool set_active_layer(Layer& layr)
{
	if(!layr)
		return false;
	
	Page pg = layr.GetPage();
	if(pg)
	{
		pg.SetShow();
		string strTemp = "page.active=";
		strTemp += layr.GetIndex() + 1;//need to labtalk is 1 offset
		
		pg.LT_execute(strTemp);
		return true;
	}
	return false;
}

//	retrive the graph layer from the tree node
GraphLayer get_graph_layer(const TreeNode& trNode)
{
	GraphLayer gl;
	if(trNode)
	{
		TreeNode trPage = trNode.Page;
		if(trPage)
		{
			GraphPage gpg(trPage.strVal);
			if(gpg)
			{
				TreeNode trTemp = trNode.Layer;
				if(trTemp)
				{
					int nIndex = trTemp.nVal;
					if(nIndex >= 0 && nIndex < gpg.Layers.Count())
					{
						gl = gpg.Layers(nIndex);
						return gl;
					}
				}
			}
		}
	}
	return gl;
}


int curve_in_layer_get_index(const curvebase& cuv, const GraphLayer& gl)
{
	int ii = 0;
	foreach(DataPlot dp in gl.DataPlots)
	{
		string strName = dp.GetDatasetName();
		if(strName.CompareNoCase(cuv.GetName()) == 0)
			return ii;
		
		ii++;
	}
	return -1;
}


BOOL	curve_in_page_get_indices(const curvebase& cuv, const GraphPage &gpg, int *pnLayerIndex, int *pnDataPlotIndex)
{
	int		nLayer = 0;
	foreach(GraphLayer gl in gpg.Layers)
	{
		int		nDPIndex = curve_in_layer_get_index(cuv, gl);
		if (0 <= nDPIndex)
		{
			if (pnLayerIndex)
				*pnLayerIndex = nLayer;
			if (pnDataPlotIndex)
				*pnDataPlotIndex = nDPIndex;
			
			return TRUE;
		}
		
		++nLayer;
	}
	
	if (pnLayerIndex)
		*pnLayerIndex = -1;
	
	if (pnDataPlotIndex)
		*pnDataPlotIndex = -1;
	
	return FALSE;
}

	
//	Search the layer where both two curves are both plotted 
//glayerFound is layer where cuv2 is currently plotted
bool curve_both_in_layer(const curvebase& cuv1, const curvebase& cuv2, GraphLayer& glayerFound)
{
	//----
	// starting layer given, then we should not search if cuv2 already in that layer and 
	// we can assume this is the layer we want
	if(glayerFound && curve_in_layer_get_index(cuv2, glayerFound) >= 0) 
		return curve_in_layer_get_index(cuv1, glayerFound) < 0? false:true;
	//----	
		
	foreach (GraphPage pg in Project.GraphPages)
	{
		foreach(GraphLayer gl in pg.Layers)
		{
			int nCuv2 = curve_in_layer_get_index(cuv2, gl);
			if(nCuv2 < 0)
				continue;
			glayerFound = gl;
			
			int nCuv1 = curve_in_layer_get_index(cuv1, gl);
			if(nCuv1 >= 0)
				return true;
		}
	}
	return false;
}

string curve_get_wks_col_names(curvebase& cc, string& strX, string& strY)
{
	string strTemp;
	string strWks;
	strY.Empty();
	strX.Empty();
	
	if(cc)
	{
		cc.GetName(strTemp);
		int nn = strTemp.Find('_');
		if(nn > 0)
			strWks = strTemp.Left(nn);
		
		Worksheet wks(strWks);
		if(!wks)
		{
			strWks.Empty();
			return strWks;
		}
			
		strY = strTemp.Mid(nn+1);
		
		if(cc.HasX(strTemp))
		{
			nn = strTemp.Find('_');
			strX = strTemp.Mid(nn+1);
		}
	}
	
	return strWks;
}

static bool check_set_col_analysis(Column& cc, bool bSetAsY = true)
{
	int nFormat = cc.GetFormat();

	if(bSetAsY && cc.GetType() != OKDATAOBJ_DESIGNATION_Y)
		cc.SetType(OKDATAOBJ_DESIGNATION_Y);
	
	if(nFormat != OKCOLTYPE_NUMERIC && nFormat != OKCOLTYPE_TEXT_NUMERIC)
	{
		cc.SetFormat(OKCOLTYPE_NUMERIC);
		return true;
	}
	return false;
}

/** >Analysis
		duplicate the active curve in the active graph layer so that analysis routine like derivatives and smoothing can be performed on the copied curve
	Parameters:
		lpcszNewColName = name of the Y column for the copied curve
		bReuseIfInWks = option to always create new copy of can reuse if named column is already in the same worksheet as the original
	Return:
		NULL if no active graph layer with active curve found, or if the operation failed
*/
curvebase& curve_duplicate_active(LPCSTR lpcszNewColName , bool bReuseIfInWks) //= NULL,  = false
{
	GraphPage gpg = Project.Pages();
	string strNewCol = lpcszNewColName;
	if(strNewCol.IsEmpty())
		strNewCol = "A";
	
	if(gpg)
	{
		GraphLayer glyr = Project.ActiveLayer();
		if(glyr)
		{
			string strY, strX;
			string strWks = curve_get_wks_col_names(Project.ActiveCurveBase(), strX, strY);
			if(!strWks.IsEmpty())
			{
				Worksheet wks(strWks);
				Column cx = wks.Columns(strX);
				Column cy = wks.Columns(strY);
				Column cNew = wks.Columns(strNewCol);
				if(!bReuseIfInWks || !cNew.IsValid())
				{
					string strTemp;
					wks.InsertCol(cy.GetIndex()+1, strNewCol, strTemp);
					strNewCol = strTemp;
					cNew = wks.Columns(strNewCol);
				}
				// we need to make sure new col is a Y col
				check_set_col_analysis(cNew);
				cNew.SetUpperBound(cy.GetUpperBound());
				return wks.GetCurve(cx.GetIndex(), cNew.GetIndex());
			}
		}
	}
	return NULL;
}

bool	add_curve_to_graph_layer(curvebase& cuv, GraphLayer &grl, int nColor, bool bRescale, int nPlotType)
{
	int			nPlot = grl.AddPlot(cuv, nPlotType);
	DataPlot	dp = grl.DataPlots(nPlot);
	dp.SetColor(nColor);

	if( bRescale )
		grl.Rescale();

	return true;
}


/**
*/
bool	curve_update_in_page(curvebase& cuv, int nColor, GraphPage &gpg, bool bRescale)
{
	if (!curve_in_page_get_indices(cuv, gpg))	// not present
	{
		GraphLayer		grl = gpg.Layers(0);
		if (!grl)
			return false;
		
		return add_curve_to_graph_layer(cuv, grl, nColor, bRescale);
	}
	else
		gpg.Refresh();

	return true;
}


// find the input curve layer and plot the result curve
bool curve_update(curvebase& cuvFit, const curvebase& cuvInput, const TreeNode& trInput, int nColor, bool bRescale) // bRescale = false
{
	GraphLayer gl = get_graph_layer(trInput);
	if(!curve_both_in_layer(cuvFit, cuvInput, gl) && gl) // there maybe many, we will just use 1st for now
	{
		/*
		// cuvInput plotted in gl but cuvFit not, so we will make the plot
		int nPlot = gl.AddPlot(cuvFit, IDM_PLOT_LINE);
		DataPlot dp = gl.DataPlots(nPlot);
		dp.SetColor(nColor);

		if( bRescale )
			gl.Rescale();
		*/

		return add_curve_to_graph_layer(cuvFit, gl, nColor, bRescale);
	}
	return false;
}


/**
*/
int convert_regular_xyz_to_matrix(vector& vecX, vector& vecY, vector& vecZ, matrix& matData, double& dXmin, double& dXmax, double& dYmin, double& dYmax, bool bCheckData)
{
	bool bRet;
	int ii, jj, id, isize, iXStepLoc, iXNumSteps, iYStepLoc, iYNumSteps, iRet;
	double dev, dXStep, dYStep;

	// Check if X,Y,Z vectors are of equal length - return if not
	int iXsize = vecX.GetSize();
	int iYsize = vecY.GetSize();
	int iZsize = vecZ.GetSize();
	if(iXsize < 4) return 1;		// support only matrix size of greater than 2x2
	if( (iXsize != iYsize) || (iYsize != iZsize) || (iZsize != iXsize) ) return 1;

	// Create a temporary worksheet, and fill in the vectors into first three columns
	Worksheet wksTemp;
	bRet = wksTemp.Create("Origin.otw", CREATE_TEMP);
	string strWksName;
	wksTemp.GetPage().GetName(strWksName);
	wksTemp.DeleteCol(0);
	wksTemp.DeleteCol(0);
	wksTemp.AddCol();
	wksTemp.AddCol();
	wksTemp.AddCol();
	Dataset dsX(wksTemp, 0);
	Dataset dsY(wksTemp, 1);
	Dataset dsZ(wksTemp, 2);
	dsX = vecX;
	dsY = vecY;
	dsZ = vecZ;

	// Create two more columns in temp worksheet for computing x/y step values
	wksTemp.AddCol();
	wksTemp.AddCol();
	Dataset dsXMed(wksTemp, 3);
	Dataset dsYMed(wksTemp, 4);

	// Find step location, number of steps, and step value for x data
	iRet = convert_regular_find_step(0, wksTemp, iXStepLoc, iXNumSteps, dXStep);
	if(iRet != 0) return 2;

	// Find step location, number of steps, and step value for y data
	iRet = convert_regular_find_step(1, wksTemp, iYStepLoc, iYNumSteps, dYStep);
	if(iRet != 0) return 2;

	// Do a check on whether the x and y data groups agree
	// Number of groups in one should match group length in the other
	if ( (iXStepLoc != iYNumSteps) | (iYStepLoc != iXNumSteps) ) return 3;

	// If data should be checked for iregularity, then...
	if(bCheckData)
	{
		// Check for x data for deviations
		iRet = convert_regular_check_data(0, wksTemp, iXStepLoc, iXNumSteps, dXStep);
		if(iRet != 0) return 4;

		// Check for y data for deviations
		iRet = convert_regular_check_data(1, wksTemp, iYStepLoc, iYNumSteps, dYStep);
		if(iRet != 0) return 4;
	}

	// Now sort the worksheet wrt Y ascending as primary and and
	// X ascending as secondary - this format is needed to copy to matrix
	using sort = LabTalk.sort;
	sort.wksname$ = strWksName;
	sort.c1 = 1;					// sort only first 3 cols
	sort.c2 = 3;
	sort.r1 = 1;
	sort.r2 = isize;
	sort.cname1$ = "A: B";
	sort.cname2$ = "A: A";
	sort.wks();

	// Now fill matrix with Z data, and return coordinate values in variables passed
	matData.SetSize(iXStepLoc, iYStepLoc);
	// Get min. max values for use in setting matrix co-ordinates
	dXmin = min(dsX);
	dXmax = max(dsX);
	dYmin = min(dsY);
	dYmax = max(dsY);
	// Fill matrix by just using Z dataset
	matData.SetByVector(dsZ);

	// Success
	return 0;
}

/**
*/
int convert_random_xyz_to_matrix_nag(vector& vecX, vector& vecY, vector& vecZ, matrix& matResult, double& dXmin, double& dXmax, double& dYmin, double& dYmax, int iMethod, double dQIL, double dWFL)
{
	int i, j, m, n, nx, ny, iRet, iErr = 0;
	double xlo, xhi, ylo, yhi;

	// Get number of scatter points, and compute grid size accordingly
	m = vecZ.GetSize();
	nx = 2 * sqrt(m);
	ny = nx;

	// Set up NAG structures for calling gridding routine
	Nag_Scat_Struct comm;
	Nag_2d_Scat_Method method;
	Nag_E01_Opt optional;
	if(iMethod == 0) method = Nag_RC;
	else method = Nag_Shep;

	// Call the appropriate interpolation routine
	if (iMethod == 0)
	{
		// Renka-Cline method
		if( (iRet = nag_2d_scat_interpolant(method, m, vecX, vecY, vecZ, &comm, NULL)) != 0) iErr = iRet;
	}
	else
	{
		// Default values of parameters for Shepard's method
		double nq = 24.0;
		double nw = 12.0;
		double rnq = -1.0;
		// Scale defaults with values passed by user; leave rnq as is
		nq *= dQIL;					// Quadratic Interpolant Locality
		nw *= dWFL;					// Weight Function Locality
		// Set these values for the conversion method
		optional.nq = nq;
		optional.nw = nw;
		optional.rnq = rnq;
		if( (iRet = nag_2d_scat_interpolant(method, m, vecX, vecY, vecZ, &comm, &optional)) != 0) iErr = iRet;
	}

	// Get lo and hi values for x and y data
	xlo = min(vecX);
	ylo = min(vecY);
	xhi = max(vecX);
	yhi = max(vecY);

	// Define vectors for gridded data, and set their size
	vector vecGX, vecGY, vecGZ;
	vecGX.SetSize(nx * ny);
	vecGY.SetSize(nx * ny);
	vecGZ.SetSize(nx * ny);

	// Compute positions of the grid points using the hi and lo scatter values
	n = 0;
	for (j = 0; j < ny; ++j)
	{
		for (i = 0; i < nx; ++i)
		{
			vecGX[i + nx * j] = (1.0 * (nx - i - 1) / (nx - 1)) * xlo + (1.0 * i / (nx - 1)) * xhi;
			vecGY[i + nx * j] = (1.0 * (ny - j - 1) / (ny - 1)) * ylo + (1.0 * j / (ny - 1)) * yhi;
			++n;
		}
	}

	// 	Evaluate two-dimensional interpolant function computed by the interpolant function call
	if( (iRet = nag_2d_scat_eval(&comm, n, vecGX, vecGY, vecGZ)) != 0)
	{
		if(iRet != 179) iErr = iRet;
	}

	// Clean up
	//--- CPY 8/4/03 saw this when doing build 644 of Origin 75
	//iRet = nag_2d_scat_free(&comm);]
	nag_2d_scat_free(&comm);
	//--- 
	// Take gridded data and create a matrix by calling regular conversion
	if((iErr == -251)||(iErr == -249)||(iErr == 0))
	{
		if( ( iRet = convert_regular_xyz_to_matrix(vecGX, vecGY, vecGZ, matResult, dXmin, dXmax, dYmin, dYmax, false)) != 0) iErr = -1;
	}

	// Return iErr
	return iErr;
}


/**
*/
int convert_sparse_xyz_to_matrix(vector& vecX, vector& vecY, vector& vecZ, matrix& matResult, double dXbegin, double dXend, double dXstep, double dYbegin, double dYend, double dYstep, int iPrecision)
{
	
	waitCursor hourGlass; 
	// Check if X,Y,Z vectors are of equal length - return if not
	int iXsize = vecX.GetSize();
	int iYsize = vecY.GetSize();
	int iZsize = vecZ.GetSize();

	if( (iXsize != iYsize) || (iYsize != iZsize) || (iZsize != iXsize) ) return -1;

	// Create a temporary worksheet to hold result data
	Worksheet wksTemp;
	bool bRet = wksTemp.Create("Origin.otw", CREATE_TEMP);
	string strWksName;
	wksTemp.GetPage().GetName(strWksName);
	wksTemp.DeleteCol(0);
	wksTemp.DeleteCol(0);
	wksTemp.AddCol();
	wksTemp.AddCol();
	wksTemp.AddCol();
	Dataset dsXresult(wksTemp, 0);
	Dataset dsYresult(wksTemp, 1);
	Dataset dsZresult(wksTemp, 2);

	// Compute length for temp worksheet columns
	int ii, jj, count;
	int ixlen = nint( 1.0 + (dXend - dXbegin) / dXstep);
	int iylen = nint( 1.0 + (dYend - dYbegin) / dYstep);
	int iSize = ixlen *iylen;
	if(iSize > MAX_MATRIX_POINTS) return -2;			// too may points - step size may be wrong
	
	dsXresult.SetSize(iSize);
	dsYresult.SetSize(iSize);
	dsZresult.SetSize(iSize);

	// Fill Z with missing values to start with
	dsZresult = NANUM;

	// Fill x, y columns of temp worksheet with cyclical values
	for (ii = 0; ii < ixlen; ii++)
	{
		for (jj = 0; jj < iylen; jj++)
		{
			dsYresult[ii * iylen + jj] = dYbegin + jj * dYstep;
		}
	}
	for (ii = 0; ii < ixlen; ii++)
	{
		for (jj = 0; jj < iylen; jj++)
		{
			dsXresult[ii * iylen + jj] = dXbegin + ii * dXstep;
		}
	}

	// Now read thru data in temp worksheet and try filling result vectors
	int ilen = iXsize;
	for (ii=0, count=0; ii < ilen; ii++)
	{
		// Get input x,y, z values
		double dX = vecX[ii];
		double dY = vecY[ii];
		double dZ = vecZ[ii];

		// Find the first occurance of the x value in the temporary worksheet
		int ixRow = Data_list(dX, &dsXresult, iPrecision);

		// If found, continue to look for y value occurance
		if(ixRow != -1)
		{
			// Find the range of rows with this x value
			double dXFound = dsXresult[ixRow];
			for(int ij = ixRow; ij < ixRow + iylen; ij++)
			{
				if(dsXresult[ij] != dXFound) break;
			}

			// Set upper and lower bounds of y dataset to range of x matches found
			dsYresult.SetLowerBound(ixRow);
			dsYresult.SetUpperBound(ij - 1);

			// Find the first y value that matches, in the above range
			int iyRow = Data_list(dY, &dsYresult, iPrecision);

			// If match found, put in the z value into the temporary worksheet at this row
			if(iyRow != -1) dsZresult[iyRow] = dZ;
			// If no match for y, increment count of bad points
			else count++;

			// Reset bounds on y dataset
			dsYresult.SetLowerBound(0);
			dsYresult.SetUpperBound(iSize - 1);
		}
		// If no match for x, increment count of bad points
		else count++;
	}

	// Now fill matrix by using Z dataset
	matResult.SetSize(iylen, ixlen);
	matResult.SetByVector(dsZresult,false);

	// Return number of discarded points
	return count;
}



/**
*/
int convert_sparse_find_min_max_step(vector& vec, double& dMin, double& dMax, double& dStep)
{
	
	// Get min and max values
	Dataset dsTemp(vec);
	dMin = min(dsTemp);
	dMax = max(dsTemp);
	int iSize = dsTemp.GetSize();
	
	// Sort datset, find diff, and sort again
	dsTemp.Sort();
	string str;	
	str = dsTemp.GetName()+"=diff("+dsTemp.GetName() + ");";
	LT_execute(str);
	dsTemp.Sort();
		
	// Use the second-last value in this sorted list as the guess for step size
	dStep = dsTemp[iSize - 2];
		
	return 0;
}



/**
*/
int xyz_remove_duplicates(vector& vecX, vector& vecY, vector& vecZ)
{
	// Check if X,Y,Z vectors are of equal length - return if not
	int iXsize = vecX.GetSize();
	int iYsize = vecY.GetSize();
	int iZsize = vecZ.GetSize();
	if(iXsize < 2) return 0;		// if size less than 2, no need to do anything
	if( (iXsize != iYsize) || (iYsize != iZsize) || (iZsize != iXsize) ) return 1;

	// Create temp worksheet and put xyz data into wks
	Worksheet wksTemp;
	bool bRet = wksTemp.Create("Origin.otw", CREATE_TEMP);
	string strWksName;
	wksTemp.GetPage().GetName(strWksName);
	wksTemp.DeleteCol(0);
	wksTemp.DeleteCol(0);
	wksTemp.AddCol();
	wksTemp.AddCol();
	wksTemp.AddCol();
	Dataset dsX(wksTemp, 0);
	Dataset dsY(wksTemp, 1);
	Dataset dsZ(wksTemp, 2);
	dsX = vecX;
	dsY = vecY;
	dsZ = vecZ;

	// Sort the worksheet wrt X ascending as primary and and Y ascending as secondary
	using sort = LabTalk.sort;
	sort.wksname$ = strWksName;
	sort.c1 = 1;					// sort the first 3 cols
	sort.c2 = 3;
	sort.r1 = 1;
	sort.r2 = iXsize;
	sort.cname1$ = "A: A";
	sort.cname2$ = "A: B";
	sort.wks();

	// Reset vectors to zero size
	vecX.SetSize(0);
	vecY.SetSize(0);
	vecZ.SetSize(0);

	// Loop thru data, look for duplicates, and replace with mean value
	int nDupRow, nTargetRow = 0, nSourceRow = 0;
	while( nSourceRow < dsX.GetSize() )
	{
		vecX.Add(dsX[nSourceRow]);
		vecY.Add(dsY[nSourceRow]);
		vecZ.Add(dsZ[nSourceRow]);

		nDupRow = nSourceRow;
		while( nSourceRow + 1 < dsX.GetSize() && dsX[nSourceRow] == dsX[nSourceRow + 1] && dsY[nSourceRow] == dsY[nSourceRow + 1] )
			nSourceRow++;

		if( nDupRow != nSourceRow )
		{
			// Replace duplicates with mean value
			for( int nRow = nDupRow + 1; nRow <= nSourceRow; nRow++ )
				vecZ[nTargetRow] += dsZ[nRow];
			vecZ[nTargetRow] /= (nSourceRow - nDupRow + 1);
		}
		nSourceRow++;
		nTargetRow++;
	}
	return 0;
}


////////////////////////////////////////////////////////////////////////////////////
// This sub function is used by regular xyz to matrix conversion function
// This function first finds the step location within the dataset, and returns the 
// step location, number of steps, and the median step value.
// To find a step, three successive differences are examined, and if the difference
// in the middle is greater than twice the difference at first and third points, 
// this is defined as a step. 
// Returns 0 for sucess, 1 if no step could be found
//
static int convert_regular_find_step(int iType, Worksheet& wks, int& iStepLoc, int& iStepNum, double& dStep)
{
	Dataset dsData(wks, iType);
	Dataset dsMed(wks, iType + 3);

	int ii, iSize;
	iSize = dsData.GetSize();
		
	// Sort worksheet wrt X/Y primary and Y/X secondary depending on iType value of 0/1
	string strWksName;
	wks.GetPage().GetName(strWksName);
	using sort = LabTalk.sort;
	sort.wksname$ = strWksName;
	sort.c1 = 1;					// sort only first 3 cols
	sort.c2 = 3;
	sort.r1 = 1; 
	sort.r2 = iSize;
	if(iType == 0)
	{	
		sort.cname1$ = "A: A"; 
		sort.cname2$ = "A: B";
	}
	else
	{
		sort.cname1$ = "A: B"; 
		sort.cname2$ = "A: A";
	}
	sort.wks();
	
	// Determine the location of the first step in data
	double step1, step2, step3;
	for(ii=0; ii < (iSize - 3); ii++)
	{
		// compute differences at three points
		step1 = fabs(dsData[ii] - dsData[ii+1]);
		step2 = fabs(dsData[ii+1] - dsData[ii+2]);
		step3 = fabs(dsData[ii+2] - dsData[ii+3]);
		if((step2 > 2 * step1) && (step2 > 2 * step3)) break;
	}
	if (ii == (iSize - 3)) return 1;
	iStepLoc = ii + 2;
	
	// Determine number of groups within X data
	iStepNum = iSize / iStepLoc;

	// Compute the step values and store them in temporary worksheet
	dsMed.SetSize(iStepNum - 1);
	for(ii = 0; ii < (iStepNum - 1); ii++)
	{
		dsMed[ii] = fabs(dsData[ii * iStepLoc + iStepLoc / 2] - dsData[(ii + 1) * iStepLoc + iStepLoc / 2]);
	}	
	
	
	// Sort these step values and pick the median value for the final x/y step size
	sort.wksname$ = strWksName;
	sort.c1 = 4 + iType;					// sort only the 4th/5th col
	sort.c2 = 4 + iType;
	sort.r1 = 1; 
	sort.r2 = iStepNum - 1;
	if(iType == 0) sort.cname1$ = "A: D";
	else sort.cname1$ = "A: E";
	sort.wks();
	dStep = dsMed[(iStepNum - 1) / 2];	
	
	// Success
	return 0;
}



////////////////////////////////////////////////////////////////////////////////////
// This sub function is used by regular xyz to matrix conversion function.
// This function checks all data points to find deviations.
// If the deviation in any point is larger than 1/4th the step size, the data is 
// rejected.
// Returns 0 on success, or 1 if a large deviation is found
// 
static int convert_regular_check_data(int iType, Worksheet& wks, int iStepLoc, int iStepNum, double dStep)
{
	Dataset dsData(wks, iType);
	Dataset dsMed(wks, iType + 3);

	int ii, jj, iSize;
	iSize = dsData.GetSize();

	// Sort worksheet wrt X/Y, ascending
	string strWksName;
	wks.GetPage().GetName(strWksName);
	using sort = LabTalk.sort;
	sort.wksname$ = strWksName;
	sort.c1 = 1;					// sort only first 3 cols
	sort.c2 = 3;
	sort.r1 = 1;
	sort.r2 = iSize;
	if(iType == 0)
	{
		sort.cname1$ = "A: A";
		sort.cname2$ = "A: B";
	}
	else
	{
		sort.cname1$ = "A: B";
		sort.cname2$ = "A: A";
	}
	sort.wks();

	// Build list of median values for X/Y groups by taking medain value of first
	// group and then using the X/Y step value to compute the rest
	dsMed.SetSize(iStepNum);
	for(ii=0; ii < iStepNum; ii++)
	{
		dsMed[ii] = dsData[iStepLoc / 2] + dStep * ii;

	}

	// Now go thru all groups of X/Y data and check for deviations.
	// Replace deviated values with median value for that group.
	// If a deviation is larger than 1/4 th of step size, reject data and return
	for(ii = 0; ii < iStepNum; ii++)
	{
		for(jj = 0; jj < iStepLoc; jj++)
		{
			int id = ii * iStepLoc + jj;
			double dev = fabs(dsData[id] - dsMed[ii]);
			if(dev >= 0.25 * dStep) return 1;
			dsData[id] = dsMed[ii];
		}
	}

	return 0;
}

////////////////////////////////////////////////////////////////////////////////////
// These functions find min and max values of a vector, and are used 
// in random_gridding_nag function.
// Should be replaced by methods in vector class when that is added
//
static double min(vector& vec)
{
	Dataset ds(vec);
	return min(ds);
}

static double max(vector& vec)
{
	Dataset ds(vec);
	return max(ds);
}

