4.2.2.10 User Defined Fitting Function using GNU Scientific LibraryUDF-GNU
This article demonstrates how to use GSL function as fit function.
Minimum Origin Version Required: Origin 8.0 SR6
1. We will fit the sample Data below by the following model:
0.1 0.10517
0.2 0.2214
0.3 0.34986
0.4 0.49182
0.5 0.64872
0.6 0.82212
0.7 1.01375
0.8 1.22554
0.9 1.4596
1 1.71828
1.1 2.00417
1.2 2.32012
1.3 2.6693
1.4 3.0552
1.5 3.48169
1.6 3.95303
1.7 4.47395
1.8 5.04965
1.9 5.68589
2 6.38906
2.1 7.16617
2.2 8.02501
2.3 8.97418
2.4 10.02318
2.5 11.18249
2.6 12.46374
2.7 13.87973
2.8 15.44465
2.9 17.17415
3 19.08554
3.1 21.19795
3.2 23.53253
2. Add the file ocgsl.h in (User Files folder), before next step, make sure the gsl dlls are copied to this same location, see Calling GNU Scientific Library.
ocgsl.h
#pragma dll(libgsl, header)
// this is OC special pragma,
// header keyword is to indicate libgsl.dll is in same location as this file
#define GSL_EXPORT // for OC, this is not needed, so make it empty
// you can directly search and copy gsl function prototypes here
typedef double (* FUNC)(double x, void * params);
struct gsl_function_struct
{
FUNC function;
void * params;
};
typedef struct gsl_function_struct gsl_function ;
typedef struct
{
size_t limit;
size_t size;
size_t nrmax;
size_t i;
size_t maximum_level;
double *alist;
double *blist;
double *rlist;
double *elist;
size_t *order;
size_t *level;
}
gsl_integration_workspace;
GSL_EXPORT gsl_integration_workspace *gsl_integration_workspace_alloc (const size_t n);
GSL_EXPORT void gsl_integration_workspace_free (gsl_integration_workspace * w);
GSL_EXPORT int gsl_integration_qag (const gsl_function * f,
double a, double b,
double epsabs, double epsrel, size_t limit,
int key,
gsl_integration_workspace * workspace,
double *result, double *abserr);
3. Press F9 to open the Fitting Function Organizer and then add a new function as follows:
4. Press the button on the right hand side of the Function Field to open the code builder and add the following codes and compile:
_nlfgsl_integration_qag.fit
#include "..\ocgsl.h"
static double f_callback(double x, void * params)
{
double alpha = *(double *)params;
return exp(alpha*x);
}
void _nlsfgsl_integration_qag(
// Fit Parameter(s):
double y0, double a, double beta,
// Independent Variable(s):
double x,
// Dependent Variable(s):
double& y)
{
// Beginning of editable part
double result, err, expected = -4.0;
// Allocates a workspace suffcient to hold 1000 double precision intervals,
// their integration results and error estimates
gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000);
gsl_function F;
F.function = f_callback;
F.params = &beta ;
// integral interval (0, x), within the desired absolute
// error 0 and relative error 1e-7
gsl_integration_qag(&F, 0, x, 0, 1e-7, 1000, 0, ww, &result, &err);
// frees the memory associated with the workspace w
gsl_integration_workspace_free (ww);
y = y0 + a*result;
// End of editable part
}
Furthermore, a more elaborate but efficient version of the fitting function is given as follows
//----------------------------------------------------------
//
#include <ONLSF.h>
#include "..\ocgsl.h"
static double f_callback(double x, void * params)
{
double alpha = *(double *)params;
return exp(alpha*x);
}
void _nlsfgsl_integration_qag(
// Fit Parameter(s):
double y0, double a, double beta,
// Independent Variable(s):
double x,
// Dependent Variable(s):
double& y)
{
// Beginning of editable part
NLFitContext *pCtxt = Project.GetNLFitContext();
if ( pCtxt )
{
static vector vInteg;
NLSFCURRINFO stCurrInfo;
pCtxt->GetFitCurrInfo(&stCurrInfo);
int nCurrentIndex = stCurrInfo.nCurrDataIndex;
BOOL bIsNewParamValues = pCtxt->IsNewParamValues();
if ( bIsNewParamValues )
{
vector vx;
pCtxt->GetIndepData(&vx);
int nSize = vx.GetSize();
vInteg.SetSize(nSize);
// Allocates a workspace suffcient to hold 1000 double precision intervals,
// their integration results and error estimates
gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000);
gsl_function F;
F.function = f_callback;
F.params = &beta ;
double result, err, expected = -4.0;
for(int ii=0; ii<nSize; ++ii)
{
// integral interval (0, vx[ii]), within the desired absolute
// error 0 and relative error 1e-7
gsl_integration_qag(&F, 0, vx[ii], 0, 1e-7, 1000, 0, ww, &result, &err);
vInteg[ii] = result;
}
// frees the memory associated with the workspace w
gsl_integration_workspace_free (ww);
}
y = y0 + a*vInteg[nCurrentIndex];
x;
}
// End of editable part
}
5. Add the following initilization codes:
Parameter Init
//Code to be executed to initialize parameters
sort( x_y_curve );
double coeff[2];
fitpoly( x_y_curve, 1, coeff);
a = coeff[0];
y0 = coeff[1];
beta=1.0
6. Fit using the user-defined function gsl_integration_qag, here are the results:
y0 = -1.06363E-6
a = 1
beta =1
|