3.9.2.2.3 Fit on one data two peaks with Gauss function with parameters W and A are shared

Version Info

Minimum Origin Version Required: Origin 8 SR0

Examples

// before running this function, import \Samples\Curve Fitting\Multiple Peaks.dat to worksheet.

#include <Origin.h>
#include <FDFTree.h>
#include <ONLSF.h>

void	NLFit_example3()
{
	string			strFDF = okutil_get_origin_path(ORIGIN_PATH_SYSTEM, "FitFunc") + "Gauss.FDF";
	Tree			trFF;
	if( !nlsf_FDF_to_tree( strFDF, &trFF ))
	{
		out_str("Fail to load FDF function file");
		return;
	}
	
	//////////////////////////////////////////////////////////////////////////////
	// 1. Setting the fitting function and parameter sharing for multipeak fit  //
	// The fit object:
	NLFit			fit;
				
	// We will fit two peaks using the same fitting function. We call this "one-replica" case.
	// This means that in addition to one full set of paramaters for the first peak, there are
	// also additional paramaters, but only those within the "replicas unit", for the second peak.
	// Only the fitting functions that allow replicas in its .FDF file can be used for multipeak fitting.
	// The [CONTROLS] section of the .FDF file must contain the following two entries to enable replicas:
	//		Duplicate Offset=2
	//		Duplicate Unit=3
	// (the actual values may be different for different functions). The value Duplicate Offset is the 1-offset index
	// of the first function paramater which belongs to the replicas unit. Duplicate Unit is the size of the replicas
	// unit. The following must hold for all such functions:
	//		total number of parameters for the function = Duplicate Offset - 1 + Duplicate Unit
	// If the above two entries are not present, then the following two entries must be present to allow
	// multipeak fitting:
	//		Replicas Offset=2
	//		Replicas Unit=3
	// The meanings are the same as described above. If you create your own fitting function that you wish
	// to use for multipeak fitting, you MUST use the latter entries.

	//
	// The multiplicity of the function is equal to the number of peaks.
	int				nFitFunctionMultiplicity = 2;

	
	// We want also to specify parameter sharing between peaks. This is accomplished via the following array:
	vector<int>		vnSharingGroups(7);
	// The size of this vector for specifying sharing between peaks must be equal to:
	// 		replicas Offset - 1 + nFitFunctionMultiplicity * replicas unit
	// In the case of the Gauss function with one replica (i.e. with two peaks) this is:
	//		2 - 1 + 2 * 3 = 7
	// This means that in the vector vnSharingGroups there is one element for each function parameter, including
	// replicas, if any: first all the function parameters for the first peak, followed by only the replicated
	// parameters (i.e. those in the replicas unit) for the additional peaks.
	// The meaning of the array vnSharingGroups: for each function paramater including replicas the value
	// says whether the argument is shared or not, and if shared, how it is shared.
	//		If the value for a particular argument is 0 or less, the parameter is not shared between replicas.
	//		If the value is greater than 0, it is shared with all the other paramaters whose value in this vector
	//		is the same. This value is the so-called group index. The concrete value is not important, as long as
	//		it is greater than 0.
	
	// We will share paramaters w and A (they are offset 2 and offset 3 in the parameter list for the Gauss function)
	// respectively, that is: w from the first peak will be the same as w for the second peak, whereas
	// A for the first peak will be the same as A for the second peak:
	vnSharingGroups[0] = 0;			// y0 ("baseline"), not shared
	vnSharingGroups[1] = 0;			// xc for the first peak, not shared
	vnSharingGroups[2] = 1;			// w, shared (group index 1) with vnSharingGroups[5] 
	vnSharingGroups[3] = 2;			// A, shared (group index 2) with vnSharingGroups[6]
	vnSharingGroups[4] = 0;			// xc for the second peak (i.e. for the first replica), not shared
	vnSharingGroups[5] = 1;			// w, shared (group index 1) with vnSharingGroups[2]
	vnSharingGroups[6] = 2;			// A, shared (group index 2) with vnSharingGroups[3]
	// Set the fitting function with the multiplicity of 2 (i.e. one replica).
	int				nn = fit.SetFunction(trFF, nFitFunctionMultiplicity, -1, vnSharingGroups, vnSharingGroups.GetSize());
	if (nn <= 0)
	{
		out_str("Failed setting fitting function!");
		return;
	}	

	
	
	///////////////////////////////////////////////////////////////////////////
	// 2. Setting data for fitting ////////////////////////////////////////////
	Worksheet 	wks = Project.ActiveLayer();
	if( !wks )
	{
		out_str("No active worksheet to get input data");
		return;
	}
	
	DataRange 	dr;
	dr.Add(wks, 0, "X");
	dr.Add(wks, 4, "Y"); // Y data from 5th column (Column E) has two peaks
	
	DWORD 		dwRules = DRR_GET_DEPENDENT | DRR_NO_FACTORS;
	int 		nNumData = dr.GetNumData(dwRules);
	ASSERT( 1 == nNumData );
	if( nNumData <= 0 || nNumData > 1)
	{
		out_str("Not proper input data to do fit");
		return;	
	}
	
	DWORD			dwPlotID;
	vector			vX, vY;
	dr.GetData( dwRules, 0, &dwPlotID, NULL, &vY, &vX);
	
	NLSFONEDATA		stDep, stIndep;	
	stIndep.pdData = vX;
	stIndep.nSize = vX.GetSize();
	stDep.pdData = vY;
	stDep.nSize = vY.GetSize();	
	
	nn = fit.SetData(1, &stDep, &stIndep);
	if (nn != 0)
	{
		out_str("SetData() failed!");
		return;
	}
	
	
	///////////////////////////////////////////////////////////////////////////
	// 3. Setting paramaters //////////////////////////////////////////////////

	// Total number of fitting parameters is 5, which is 7 (see above) reduced by two due to two
	// shared paramaters.
	int				nNumParams = 5;
	
	vector			vParams(nNumParams);			// this vector should be initialized to the total number of paramaters
	// This vector will be used both to set initial paramater values, and
	// to receive the parameter values after fitting:
	vParams[0] = 1.66;				// y0, not shared 
	vParams[1] = 2.5;				// xc, not shared
	vParams[2] = 0.3;				// w, shared (group 1)
	vParams[3] = 50;				// A, shared (group 2)
	vParams[4] = 7.5;				// xc for the second peak, not shared
	nn = fit.SetParams(nNumParams, vParams);
	if ( 0 != nn )
	{
		out_str("Invalid paramater setting!");
		return;
	}	
	
	
	///////////////////////////////////////////////////////////////////////////
	// 4. Fitting /////////////////////////////////////////////////////////////
	int				nMaxNumIterations = 100;
	nn = fit.Fit(nMaxNumIterations);
	printf("Fit(%ld) returned: %ld\n", nMaxNumIterations, nn);
	if( 0 == nn )
		printf("Fit Done!\n");
	
	
	///////////////////////////////////////////////////////////////////////////
	// 5. Dump all the results and compare with what was expected /////////////
	for (int iparam = 0; iparam < nNumParams; iparam++)
	{
		printf("param index = %d   value obtained = %lf\n", iparam + 1, vParams[iparam]);
	}
	
	return;
	
}