Curve fitting with a system of equations + convolution
I am fitting data that is a convolution of a gaussian instrument response with exponential decay / growth associated with a kinetic model. Because we often change models, I use integrateODE to solve the kinetic model for populations at each fit point (the signal will be the sum of amplitudes*populations) , and the convolution is done numerically from zero up to each fit point, within the instrument response region. The fit parameters end up being rate constants within the kinetic model and the populations' amplitudes.
This is pretty slow because we're repeating the convolution / model integration a lot, and I further use it for global fitting.
Generally, the decay times we see span ~4 orders of magnitude, so the fit needs to have high resolution at shorter times. I can pretty easily write something to do all-at-once fits if I know what the fit points will be, but all-at-once seems to choose evenly spaced points throughout the data range. Is there a way to control this or a good way to account for uneven data spacing within the fit? I think the way I'm doing this is clunky and I am looking for an efficient way to do the fits.
Here is an example fit function. I have a separate mode that has the population connectivities. I can provide an example with data if that would be helpful.
Thanks in advance
Function NRF_SHORT_BR1(w,x) : FitFunc //A-(.95)>B->C->D
Wave w //A(.05)->E E is S0
Variable x
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(x) =
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ x
//CurveFitDialog/ Coefficients 13
//CurveFitDialog/ w[0] = Y_offset
//CurveFitDialog/ w[1] = tzero
//CurveFitDialog/ w[2] = width
//CurveFitDialog/ w[3] = amp1
//CurveFitDialog/ w[4] = tau1
//CurveFitDialog/ w[5] = amp2
//CurveFitDialog/ w[6] = tau2
//CurveFitDialog/ w[7] = amp3
//CurveFitDialog/ w[8] = tau3
//CurveFitDialog/ w[9] = amp4
//CurveFitDialog/ w[10] = amp5
variable conv, step, t, points
Variable V_fitOptions=4
variable V_Fiterror = 0
wave pw, y_w
make /o/n=3 pw
make/o/n=(2,5) y_w
pw = {W[4],W[6],W[8]} //these are the time constants
y_w[0][0] = .66 //populations for the model
y_w[1][0] = 0
y_w[][1] = 0
y_w[][2]=0
y_w[0][3]=.34
y_w[1][3] = 0
y_w[0][4] = 0
y_w[1][4] = 0
//y_w[][5] = 0
make/o/n=(2,6) xwave
xwave[0][] = w[1]
xwave[1][] = x
make/o/n=(2,6) stepwave
//y_w*=y_w*w[p]
points = ceil(w[2]/3E-4) //sets the numper of points in wave according to the width of pulse
make/o/n = (points) gwave //gaussian wave
make/o/n = (points) tgwave //x-wave for gaussian wave
make/o/n = (points) fwave //fitting wave
make/o/n = (points) convwave //convolved wave
if (x<=(w[1]-w[2]*6)) //check to see if time point is less than 6 sigma away from time zero
conv = 0
elseif (x>=(w[1]+w[2]*6)) // if greater than 6 sigma for the numerical convolution range solve only the exponential only
integrateODE /x=xwave model_br1, pw, Y_W
conv = w[3]*y_w[1][0] + w[5]*y_w[1][1] + w[7]*y_w[1][2]+w[10]*y_w[1][4] + (1-y_w[1][3])*w[9] //integrate using the model up to build the signal
else // solve the numerical convolution
step = -w[2]*6
t = 0
do
if ((step+x)>=w[1])
xwave[1][] = x + step
integrateODE/x=xwave model_br1,pw, Y_W
fwave[t] =w[3]*y_w[1][0] + w[5]*y_w[1][1] + w[7]*y_w[1][2]+w[10]*y_w[1][4] + (1-y_w[1][3])*w[9] //integrate to build the signal
else
fwave[t] = 0 //before time zero is 0
endif
gwave[t] = (0.5*exp(-((step)^2/(2*w[2]^2)))) //build gaussian wave
tgwave[t] = step
step+= ((w[2]*6)*2)/(points)
t+=1
while (t<points)
tgwave=tgwave+x
gwave = gwave*(2/(sqrt(2*pi)*w[2])) //converting gaussian wave into a normalized gaussian wave
convwave = fwave
convolve/a gwave, convwave //convolve the fit wave and the guassian wave witht the convwave as the desitation wave
convwave = convwave*(((w[2]*6)*2)/points) //need to account for step size between points
conv = (convwave[(numpnts(fwave)/2)-1])
endif
return conv
End
Assuming that your fitting process is optimized (I'm not commenting on that part), you would want your fitting function to execute as fast as possible. Looking at your code I see many possible optimizations.
1. At the top of your function you have a bunch of Make commands. Ideally, you should have one set of Make(s) like this for each running thread and have it called outside the fitting function.
2. Inside the function compute w[2]*6 only once and assign to a local variable.
3. Inside your do loop you the quantity 2*w[2]^2 does not change so you should factor it out of the loop and assign to a local variable.
4. Similary, 2/(sqrt(2*pi)*w[2] does not change. In fact, you should also replace the sqrt(2*pi) with something like
5. Replace gwave=... with a MatrixOP assignment.
6. Delete the convwave=fwave line.
7. Replace the Convolve and the following line with a single MatrixOP assignment.
I hope this helps.
A.G.
July 29, 2021 at 12:00 pm - Permalink
You can use an x-wave with an all-at-once fit function, see DisplayHelpTopic "User-Defined Fitting Function: Detailed Description" Pay attention to the last line.
All-At-Once Fitting Functions
The scheme of calculating one Y value at a time doesn't work well for some fitting functions. This is true of functions that involve a convolution such as might arise if you are trying to fit a theoretical signal convolved with an instrument response. Fitting to a solution to a differential equation might be another example.
For this case, you can create an "all at once" fit function. Such a function provides you with an X and Y wave. The X wave is input to the function; it contains all the X values for which your function must calculate Y values. The Y wave is for output -- you put all your Y values into the Y wave.
Because an all-at-once function is called only once for a given set of fit coefficients, it will be called many fewer times than a basic fit function. Because of the saving in function-call overhead, all-at-once functions can be faster even for problems that don't require an all-at-once function.
Here is the format of an all-at-once fit function:
Function myFitFunc(pw, yw, xw) : FitFunc
WAVE pw, yw, xw
yw = <expression involving pw and xw>
End
Note that there is no return statement because the function result is put into the wave yw. Even if you include a return statement the return value is ignored during a curve fit.
The X wave contains all the X values for your fit, whether you provided an X wave to the curve fit or not. If you did not provide an X wave, xw simply contains equally-spaced values derived from your Y wave's X scaling.
July 29, 2021 at 12:20 pm - Permalink
A quick comment about gwave = gwave*(2/(sqrt(2*pi)*w[2]))
In addition to calculating the constant beforehand as suggested by Igor you could also use FastOP. It's very fast when you only need to add or multiply with a constant.
July 29, 2021 at 12:30 pm - Permalink
First- YES you need an all-at-once fit function. You are integrating over an X range over and over and over...
Since you are using Convolve, your model must have evenly-space X values. But since your data is not evenly spaced, you will need to make an evenly-spaced model as an intermediate step, and then pick out the appropriate Y values based on your xw wave.
If your system has decay constants that vary of orders of magnitude, you may have a "stiff" system. In that case, use IntegrateODE/M=3 to use the BDF method of integration.
Have you run IntegrateODE outside of the fitting function? That is pretty essential to make sure that part of the computation is working correctly. It can be quite difficult to get it right.
y_w[1][0] = 0
y_w[][1] = 0
y_w[][2]=0
y_w[0][3]=.34
y_w[1][3] = 0
y_w[0][4] = 0
y_w[1][4] = 0
Initial conditions go in the first (zero) row of the output wave. The second line here sets the second row (row 1, of course :). If your goal is to initialize the wave to all zeroes, you probably want something like this:
// Now set the initial conditions in the zero row
y_w[0][0] = .66 //populations for the model
y_w[0][1] = 0
y_w[0][2] = 0
y_w[0][3] = .34
y_w[0][4] = 0
olelytken has quoted a part of the documentation for all-at-once fitting functions. To read all about it, execute this command in Igor: DisplayHelpTopic "All-At-Once Fitting Functions"
It sounds like you have a pretty difficult project. Try to apply what we have suggested here, then come back. I'm sure you will run into further blocks.
The advice Igor is giving you is very good and sound. But I would get the basic implementation doing the right thing first, and then work on the types of optimizations he is suggesting.
July 29, 2021 at 12:47 pm - Permalink