optimizing user fit function
sleepingawake86
I am trying to fit data which is modeled by a pair of coupled partial differential equations(CPDEs). I wrote code to solve the CPDEs and it works well, and I implemented it into a user fit function and it fits fine, but it is extremely slow. The issue is that I end up having to do a bunch of large matrix(typically 500x500) calculations for the solution and its very slow.
I have attached the code below, any help optimizing the code would be great! thanks!
function intpopmodel()
NVAR a,b,su,sd,jmax,kmax, dz,dt,I0,gt
variable jj,kk
make /o /d /n=(jmax,kmax) nwave, dn, intwave, dI
nwave[0][]=1
intwave[][0]=I0
For(jj=0;jj<=jmax-1;jj+=1)
For(kk=0;kk<=kmax-1;kk+=1)
nwave[jj][0]=b/(b+a*I0)+a*I0/(b+a*I0)*exp(-(b+a*I0)*jj*dt)
intwave[0][kk]=I0*exp(-su*kk*dz)
endfor
endfor
For(jj=0;jj<=jmax-1;jj+=1)
For(kk=0;kk<=kmax-1;kk+=1)
dn[jj][kk]=(-a*intwave[jj][kk]*nwave[jj][kk]+b*(1-nwave[jj][kk]))*dt
dI[jj][kk]=-su*intwave[jj][kk]*nwave[jj][kk]*dz-sd*intwave[jj][kk]*(1-nwave[jj][kk])*dz
nwave[jj+1][kk]=nwave[jj][kk]+dn[jj][kk]
intwave[jj][kk+1]=intwave[jj][kk]+dI[jj][kk]
endfor
endfor
variable Nt=gt/dt
variable intend=intwave[Nt][kmax]/I0
return intend
end
NVAR a,b,su,sd,jmax,kmax, dz,dt,I0,gt
variable jj,kk
make /o /d /n=(jmax,kmax) nwave, dn, intwave, dI
nwave[0][]=1
intwave[][0]=I0
For(jj=0;jj<=jmax-1;jj+=1)
For(kk=0;kk<=kmax-1;kk+=1)
nwave[jj][0]=b/(b+a*I0)+a*I0/(b+a*I0)*exp(-(b+a*I0)*jj*dt)
intwave[0][kk]=I0*exp(-su*kk*dz)
endfor
endfor
For(jj=0;jj<=jmax-1;jj+=1)
For(kk=0;kk<=kmax-1;kk+=1)
dn[jj][kk]=(-a*intwave[jj][kk]*nwave[jj][kk]+b*(1-nwave[jj][kk]))*dt
dI[jj][kk]=-su*intwave[jj][kk]*nwave[jj][kk]*dz-sd*intwave[jj][kk]*(1-nwave[jj][kk])*dz
nwave[jj+1][kk]=nwave[jj][kk]+dn[jj][kk]
intwave[jj][kk+1]=intwave[jj][kk]+dI[jj][kk]
endfor
endfor
variable Nt=gt/dt
variable intend=intwave[Nt][kmax]/I0
return intend
end
This double loop has two statements. However, both statements are independent in the sense that one depends only on jj and the other depends only on kk. So rewrite it like this:
nwave[jj][0]=b/(b+a*I0)+a*I0/(b+a*I0)*exp(-(b+a*I0)*jj*dt)
endfor
for(kk=0;kk<=kmax-1;kk+=1)
intwave[0][kk]=I0*exp(-su*kk*dz)
endfor
Which will be much faster because you're now doing kmax + jmax assignments instead of kmax * jmax.
Rewrite this using MatrixOP. For example:
MatrixOP /O dn = (-a * intwave * nwave +b* (M_Ones - nwave)) * dt
MatrixOP /O dI = -su*intwave * nwave *dz -sd * intwave * (M_Ones - nwave) * dz
Note that you can't use MatrixOP directly for the other two assignments due to the '+1' index in
nwave[jj+1][kk]=nwave[jj][kk]+dn[jj][kk]. So either keep the full loop or do some trickery.That said: you're only returning a single value from your function:
So how about you avoid calculating the full matrices, and you calculate only the one point that will be used? That will be the biggest speedup by far.
By the way, you may want to convert the globals (NVAR) in variables that are passed into the function, or variables that are passed in using a parameter wave.
September 21, 2012 at 12:33 pm - Permalink