Fit on one data two peaks with two functions
Version Info
Minimum Origin Version Required: Origin 8 SR1
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_example4()
{
// Two functions:
// (1) Gauss
string strFDF1 = okutil_get_origin_path(ORIGIN_PATH_SYSTEM, "FitFunc") + "Gauss.FDF";
Tree trFF1;
if( !nlsf_FDF_to_tree( strFDF1, &trFF1 ))
{
out_str("Fail to load FDF function file");
return;
}
// (2) Lorentz
string strFDF2 = okutil_get_origin_path(ORIGIN_PATH_SYSTEM, "FitFunc") + "Lorentz.FDF";
Tree trFF2;
if( !nlsf_FDF_to_tree( strFDF2, &trFF2 ))
{
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 to a total of four peaks, two Gaussian and two Lorentzian peaks,
// with some parameter sharing between peaks.
// Each function will therefore have the multiplicity of 2 (i.e. one replica in each function).
//
// The first function has 4 original paramaters, 3 of which are in replicas unit.
// Since the function will have multiplicity of 2, this means that per function
// (i.e. for its two peaks) there are 7 paramaters:
// 7 = replicas Offset - 1 + nFitFunctionMultiplicity * replicas unit
// since for both Gauss and Lorentz replicas Offset = 2, replicas unit = 3
// and we use nFitFunctionMultiplicity = 2.
int nFitFunctionMultiplicity = 2;
// To set up parameter sharing between peaks, we use an auxiliary vector of integers:
// First we will set the first function (Gauss).
// The first function has 7 parameters
vector<int> vnSharingGroupsGauss(7);
vnSharingGroupsGauss[0] = 0; // y0 (baseline), not shared
vnSharingGroupsGauss[1] = 0; // xc for the first Gauss peak, not shared
vnSharingGroupsGauss[2] = 1; // w for the first Gauss peak, shared with vnSharingGroupsLorentz[4] (see below)
vnSharingGroupsGauss[3] = 2; // A for the first Gauss peak, shared with vnSharingGroupsLorentz[2] (see below)
vnSharingGroupsGauss[4] = 0; // xc for the second Gauss peak, not shared
vnSharingGroupsGauss[5] = 3; // w for the second Gauss peak, shared with vnSharingGroupsLorentz[1] (see below)
vnSharingGroupsGauss[6] = 4; // A for the second Gauss peak, shared with vnSharingGroupsLorentz[5] (see below)
// Set the Gauss fitting function with the multiplicity of 2 (i.e. one replica).
int nn = fit.SetFunction(trFF1, nFitFunctionMultiplicity, -1, vnSharingGroupsGauss, vnSharingGroupsGauss.GetSize());
if (nn <= 0)
{
out_str("Failed setting fitting function!");
return;
}
// The second function does not have the baseline parameter y0 (the baseline parameters are those that
// do not belong to ther replicas unit, both Gauss and Lorentz have exactly one such paramater - y0).
// (only one baseline makes sense and it is always associated with the first function)
// so the total number of params for this function is 6.
// This is the auxiliary vector of integers for specifying the parameter sharing
// for the Lorentz peaks:
// The
vector<int> vnSharingGroupsLorentz(6);
vnSharingGroupsLorentz[0] = 0; // xc for the first Lorentx peak, not shared
vnSharingGroupsLorentz[1] = 3; // w for the first Lorentz peak, shared with vnSharingGroupsGauss[5]
vnSharingGroupsLorentz[2] = 2; // A for the first Lorentz peak, shared with vnSharingGroupsGauss[3]
vnSharingGroupsLorentz[3] = 0; // xc for the second Lorentz peak, not shared
vnSharingGroupsLorentz[4] = 1; // w for the second Lorentz peak, shared with vnSharingGroupsGauss[2]
vnSharingGroupsLorentz[5] = 4; // A for the second Lorentz peak, shared with vnSharingGroupsGauss[6]
// Append the Lorentz fitting function with the multiplicity of 2 (i.e. one replica).
nn = fit.SetFunction(trFF2, nFitFunctionMultiplicity, -1, vnSharingGroupsLorentz, vnSharingGroupsLorentz.GetSize());
if (nn <= 0)
{
out_str("Failed setting fitting function!");
return;
}
///////////////////////////////////////////////////////////////////////////
// 3. 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, 2, "Y"); // Y data from 3th column (Column C) has four 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;
}
///////////////////////////////////////////////////////////////////////////
// 4. Setting paramaters //////////////////////////////////////////////////
// Total number of fitting parameters is dependent on the parameter sharing between peaks.
// It is crucial that the number of parameters as specifed
// in the SetParams() call below agrees with the numbers of parameters in the functions and with the parameter
// sharing as specifed in step 2.
//
// This is how the parameters are lined up:
// First all the parameters for the first function. Their total number is the same as the total
// number of parameters for one Gauss function with one replica, which is, as explained above, 7. Note
// that there is no paramater sharing between paramaters within Gauss function peaks alone, since all
// the values in the vector vnSharingGroupsGauss are different (with the exception of 0, which, by
// denition, means that those parameters are not shared).
// Now we need to add the parameters for Lorentz. The vector vnSharingGroupsLorentz has 6 elements,
// but 4 of the 6 have nonzero values which have already appeared in vnSharingGroupsGauss (which means they
// are shared with some paramaters in Gauss), so only 2 more parameters (vnSharingGroupsLorentz[0]
// and vnSharingGroupsLorentz[3]) are contributed by Lorentz peaks.
// That amounts to a total of 9.
int nNumParams = 9;
vector vParams(nNumParams); // this vector should be initialized to the total number of paramaters
// The vector vParams will be used to supply the initial values of parameters, and it will also receive
// the parameter values after fitting.
// The paramaters are lined up in this vector such that the first come all the parameters for the Gauss, and
// then the additional ones for the Lorentz:
vParams[0] = 12; // the baseline y0 (vnSharingGroupsGauss[0])
vParams[1] = 6.3; // xc for the first Gauss peak (vnSharingGroupsGauss[1])
vParams[2] = 0.1; // w for the first Gauss peak (vnSharingGroupsGauss[2]) and for the second Lorentz peak (vnSharingGroupsLorentz[4]).
vParams[3] = 85; // A for the first Gauss peak (vnSharingGroupsGauss[3]) and for the first Lorentz peak (vnSharingGroupsLorentz[2]).
vParams[4] = 8; // xc for the second Gauss peak (vnSharingGroupsGauss[4])
vParams[5] = 0.1; // w for the second Gauss peak (vnSharingGroupsGauss[5]) and for the first Lorentz peak (vnSharingGroupsLorentz[1]).
vParams[6] = 40; // A for the second Gauss peak (vnSharingGroupsGauss[6]) and for the second Lorentz peak (vnSharingGroupsLorentz[5]).
vParams[7] = 3.8; // xc for the first Lorentz peak (vnSharingGroupsLorentz[0])
vParams[8] = 1.2; // xc for the second Lorentz peak (vnSharingGroupsLorentz[3])
nn = fit.SetParams(nNumParams, vParams);
if ( 0 != nn )
{
out_str("Invalid paramater setting!");
return;
}
///////////////////////////////////////////////////////////////////////////
// 5. Fitting /////////////////////////////////////////////////////////////
int nMaxNumIterations = 100;
nn = fit.Fit(nMaxNumIterations);
printf("Fit(%ld) returned: %ld\n", nMaxNumIterations, nn);
if( 0 == nn )
printf("Fit Done!\n");
///////////////////////////////////////////////////////////////////////////
// 6. 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]);
}
_calc_fitted_data_for_multi_funcs_fit(trFF1, trFF2, vParams, vX, wks);
return;
}
static void _calc_fitted_data_for_multi_funcs_fit(const TreeNode& trFF1, const TreeNode& trFF2, const vector& vParams, const vector& vX, Worksheet& wks)
{
vector vOnePeakParams;
vector vY1, vY2, vY3, vY4; // Y data for each peak
int nNumPts = vX.GetSize();
// calculate for first peak
vOnePeakParams.SetSize(4); // Gauss and Lorentz original both have 4 parameters
vOnePeakParams[0] = vParams[0];
vOnePeakParams[1] = vParams[1];
vOnePeakParams[2] = vParams[2];
vOnePeakParams[3] = vParams[3];
NumericFunction NF;
if(!NF.SetTree(trFF1))
{
out_str("Invalid function");
return;
}
vY1.SetSize(nNumPts);
BOOL bb = NF.Evaluate(vOnePeakParams, vY1, vX, NULL, nNumPts);
if (!bb)
{
out_str("Failed!");
return;
}
// calculate for second peak
vOnePeakParams[0] = vParams[0];
vOnePeakParams[1] = vParams[4];
vOnePeakParams[2] = vParams[5];
vOnePeakParams[3] = vParams[6];
vY2.SetSize(nNumPts);
bb = NF.Evaluate(vOnePeakParams, vY2, vX, NULL, nNumPts);
if (!bb)
{
out_str("Failed!");
return;
}
// calculate for third peak
vOnePeakParams[0] = vParams[0];
vOnePeakParams[1] = vParams[7];
vOnePeakParams[2] = vParams[5];
vOnePeakParams[3] = vParams[3];
if(!NF.SetTree(trFF2))
{
out_str("Invalid function");
return;
}
vY3.SetSize(nNumPts);
bb = NF.Evaluate(vOnePeakParams, vY3, vX, NULL, nNumPts);
if (!bb)
{
out_str("Failed!");
return;
}
// calculate for fourth peak
vOnePeakParams[0] = vParams[0];
vOnePeakParams[1] = vParams[8];
vOnePeakParams[2] = vParams[2];
vOnePeakParams[3] = vParams[6];
vY4.SetSize(nNumPts);
bb = NF.Evaluate(vOnePeakParams, vY4, vX, NULL, nNumPts);
if (!bb)
{
out_str("Failed!");
return;
}
vector vYOut;
vYOut = vY1 + vY2 + vY3 + vY4 - vParams[0] * 3; // added fitted value for 4 peaks and reduce duplicate y0
// put fitted data into new column
int nCol = wks.AddCol();
Dataset dsOut(wks, nCol);
dsOut = vYOut;
printf("Fitted data putted into Column %s\n", wks.Columns(nCol).GetName());
return;
}
|