Minimum Origin Version Required: Origin 8 SR1
// 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; }