| Author |
 Topic  |
|
|
nmr2013
Germany
3 Posts |
 Posted - 06/18/2013 : 06:49:10 AM
|
Origin Ver. and Service Release (Select Help-->About Origin): 8.5.0G SR1
Operating System: Windows XP
Hi everybody,
I'm trying to use a convolution of a rectangular and a gaussian as fitting function for some experimental data. I used the codes given in the online tutorials (link1 and link2) as a starting point, but unfortunately I can't get the code to work properly, even the one given in the tutorials doesn't produce a proper fit when I follow it. I printed the functions and their convolution results to a datasheet and this part seems to work properly. Unfortunately fitting only calculates the first x/y pair and then stops saying that the fit doesn't converge. The online tutorial files produce the same error (only the first x/y pair is calculated) and I can't identify the reason. Any ideas?
Thanks a lot
Sven
________________________
#pragma warning(error : 15618)
#include <origin.h>
#include <ONLSF.h>
#include <fft_utils.h>
NLFitContext *pCtxt = Project.GetNLFitContext();
if (pCtxt )
{
Worksheet wks = Project.ActiveLayer();
vector<double> vNu;
Dataset Nu(wks, 0);
Dataset test1(wks, 2);
Dataset test2(wks, 3);
Dataset test3(wks, 4);
Dataset test4(wks, 5);
Dataset test5(wks, 6);
Dataset test6(wks, 7);
vNu=Nu;
int iSize=vNu.GetSize();
vector<double> vGauss(iSize), vA(iSize), vFalt(2*iSize-1), vFalt2(2*iSize-1), vRet(iSize);
vGauss=A*exp(-1.0*(vNu-mu)^2/(2.0*(M2*10^3)^2));
double lower=mu-d;
double upper=mu+d;
for(int ii=0; ii<iSize; ii++) { vA[ii]=A;}
for(ii=0; ii<iSize; ii++) {
if(vNu[ii]<lower) vA[ii]=0;
if(vNu[ii]>upper) vA[ii]=0;
continue;
}
BOOL bIsNewParamValues = pCtxt->IsNewParamValues();
if ( bIsNewParamValues )
{
vFalt2.SetSubVector(vA);
vFalt.SetSubVector(vGauss);
test1=vNu; test2=vGauss; test3=vA; test4=vFalt; test5=vFalt2;
int iRet = fft_fft_convolution(2*iSize-1, vFalt, vFalt2);
vFalt.GetSubVector(vRet, floor(iSize/2), iSize + floor(iSize/2)-1 );
test6=vRet;
}
NLSFCURRINFO stCurrInfo;
pCtxt->GetFitCurrInfo(&stCurrInfo);
int nCurrentIndex = stCurrInfo.nCurrDataIndex;
y=vRet[nCurrentIndex];
x;
}
// End of editable part
} |
|
|
Penn
China
644 Posts |
Posted - 06/19/2013 : 02:26:55 AM
|
Hi Sven,
I have tried both tutorials you mentioned, and they work. Now, you are going to merge these two tutorials. Though the fitting function is able to define and compile successfully, without your data, I am not sure whether it works or not. If you can provide your data, definition of your fitting function, and the proper initial values for each parameters, and how you want to fit your data, then we can have a try.
Penn |
 |
|
|
nmr2013
Germany
3 Posts |
Posted - 06/19/2013 : 08:00:12 AM
|
Hi,
thanks for your help. Do you have the same Origin Version?
The desired fitting function is simply the convolution of a gaussian and a rectangular, as given in the code.
The parameters I used are to be interpreted as follows:
- A: Intensity of the functions (good starting value depends on dataset, the maximum intensity value should be ok)
- mu: Position of the gaussian/Rectangular (good starting value is around the position of the maximum of the experimental data)
- d: half width of the rectangular (s. lower and upper in the code, should be comparable to the full width at half height of the peak)
- M2: second moment of the gaussian (-> width parameter, 10 should be a good value)
The ASCII of the dataset can be downloaded here: http://rapidshare.com/files/61950759/spect.txt
Thanks again, your help is very much appreciated!
Sven |
Edited by - nmr2013 on 06/19/2013 08:16:05 AM |
|
|
|
Penn
China
644 Posts |
Posted - 06/20/2013 : 03:52:29 AM
|
Hi Sven,
I tried both tutorials in the same version as you are using, and both work.
If just using your code directly, there is a speed problem, so I change the code as below to improve the speed.
NLFitContext *pCtxt = Project.GetNLFitContext();
if (pCtxt )
{
static vector vRet;
BOOL bIsNewParamValues = pCtxt->IsNewParamValues();
if ( bIsNewParamValues )
{
Worksheet wks = Project.ActiveLayer();
vector<double> vNu;
Dataset Nu(wks, 0);
Dataset test1(wks, 2);
Dataset test2(wks, 3);
Dataset test3(wks, 4);
Dataset test4(wks, 5);
Dataset test5(wks, 6);
Dataset test6(wks, 7);
vNu=Nu;
int iSize=vNu.GetSize();
//vector<double> vGauss(iSize), vA(iSize), vFalt(2*iSize-1), vFalt2(2*iSize-1), vRet(iSize);
vector<double> vGauss(iSize), vA(iSize), vFalt(2*iSize-1), vFalt2(2*iSize-1);
vRet.SetSize( iSize );
vGauss=A*exp(-1.0*(vNu-mu)^2/(2.0*(M2*10^3)^2));
double lower=mu-d;
double upper=mu+d;
for(int ii=0; ii<iSize; ii++) { vA[ii]=A;}
for(ii=0; ii<iSize; ii++) {
if(vNu[ii]<lower) vA[ii]=0;
if(vNu[ii]>upper) vA[ii]=0;
continue;
}
vFalt2.SetSubVector(vA);
vFalt.SetSubVector(vGauss);
test1=vNu; test2=vGauss; test3=vA; test4=vFalt; test5=vFalt2;
int iRet = fft_fft_convolution(2*iSize-1, vFalt, vFalt2);
vFalt.GetSubVector(vRet, floor(iSize/2), iSize + floor(iSize/2)-1 );
test6=vRet;
}
NLSFCURRINFO stCurrInfo;
pCtxt->GetFitCurrInfo(&stCurrInfo);
int nCurrentIndex = stCurrInfo.nCurrDataIndex;
y=vRet[nCurrentIndex];
x;
}
There are two main changes (do like the two tutorials) in the code to improve the speed.
1. Move most of the code into the if(bIsNewParamValues) branch.
2. vRet is defined as static.
With the following initial values for parameters, It is able to get the final results in 8.5.0 SR1, so you can have a try.
A = 4
mu = 777
d = 5249
M2 = 10
Also, there are some issues to be payed attention to.
1. To do such convolution as fitting function, the X data is better with the ascending order and even-space.
2. Please check whether you need to consider the dx, as the second tutorial does.
Penn |
Edited by - Penn on 06/20/2013 03:56:36 AM |
|
|
|
nmr2013
Germany
3 Posts |
Posted - 06/24/2013 : 03:39:27 AM
|
| Thank you very much, the code you posted worked just fine! |
|
|
| |
Topic |
|
|
|
Using a convolution for fitting