0

I noticed that the following questions were related: Variable Parameters Curve Fit 0, Variable Parameters Curve Fit 1, Hold Parameter Curve Fit 0.

It creates an array of guess parameters and array stipulating which parameters to hold constant. It then generate the correct name to call and then use a wrapper to call it. It works but I want to put my definitions in a separate file called utilites.py and then call the definitions from the initial file. As you see wrapper() needs to find wrapperName(). Is there some way to pass wrapper name to a separate file and then make it accessible by all the definitions in utilities.py? If you try to pass wrapper an extra parameter curve_fit() breaks.

import numpy as np
from numpy import random
import scipy as sp
from scipy.optimize import curve_fit
import matplotlib.pyplot as pl

def exp1(t,a1,tau1):
        val=0.
        val=(a1*np.exp(-t/tau1))*np.heaviside(t,0)
        return val

def exp2(t,a1,tau1,a2,tau2):
        val=0.
        val=(a1*np.exp(-t/tau1)+a2*np.exp(-t/tau2))*np.heaviside(t,0)
        return val

def exp3(t,a1,tau1,a2,tau2,a3,tau3):
        val=0.
        val=(a1*np.exp(-t/tau1)+a2*np.exp(-t/tau2)+a3*np.exp(-t/tau3))*np.heaviside(t,0)
        return val

def exp4(t,a1,tau1,a2,tau2,a3,tau3,a4,tau4):
        val=0.
        val=(a1*np.exp(-t/tau1)+a2*np.exp(-t/tau2)+a3*np.exp(-t/tau3)+a3*np.exp(-t/tau3))*np.heaviside(t,0)
        return val

def multiexp(t,args):
        val=np.zeros(len(t))
        for i in np.arange(0,len(args),2):
            val+=args[i]*np.exp(-t/args[i+1])*np.heaviside(t,0)
        return val

def genWrapperName(hold,p0):
    counter=0
    if len(p0)==2:
        wrapperName='exp1(t,'
    elif len(p0)==4:
        wrapperName='exp2(t,'
    elif len(p0)==6:
        wrapperName='exp3(t,'
    elif len(p0)==8:
        wrapperName='exp4(t,'
    for i in range(0, len(hold)):
        if hold[i]:
            wrapperName+=str(p0[i])
        else:
            if i%2==0:
                wrapperName+='args['+str(counter)+']'
                counter+=1
            else:
                wrapperName+='args['+str(counter)+']'
                counter+=1
        if i<len(hold)-1:
            wrapperName+=','
    wrapperName+=')'
    return wrapperName

def wrapper(t,*args):
    global wrapperName
    return eval(wrapperName)

#inital guess parameters
p0=np.array([1.5,200.,2.,700.,3.,2400.])
#list of parameters to hold in the least squares fit. 1 means hold
hold=[0,1,0,1,0,0]
holdloc=[i for i, x in enumerate(hold) if x]
p1=np.delete(p0,holdloc)
timepoints = np.arange(20.,4000.,20.)
wrapperName=genWrapperName(hold,p0)
y=exp3(timepoints,1.,300.,2.,800,1.,3000.)+np.random.normal(0, .1, size=len(timepoints))

popt, pcov = curve_fit(exp3, timepoints, y, p0=p0)
print 'unheld parameters:', popt
print 'unheld covariance:', pcov

popt, pcov = curve_fit(wrapper, timepoints, y, p0=p1)

for i in range(0, len(hold)):
    if hold[i]:
        popt=np.insert(popt,i,p0[i])
print 'held parameters:', popt
print 'held covariance:', pcov
yfit=multiexp(timepoints,popt)

pl.plot(timepoints,y,timepoints,yfit)
pl.show()
Scott Ellis
  • 1
  • 1
  • That's a lot of code. Do you have something simple that demonstrates the problem? – tdelaney Apr 15 '18 at 20:13
  • You put some stuff in utilities.py... you can use the same trick for shared functions. In fact, utilities.py seems like a good option. – tdelaney Apr 15 '18 at 20:13

1 Answers1

0

In holdCurveFit.py

import numpy as np
import utilities as ut
from numpy import random
import scipy as sp
from scipy.optimize import curve_fit
import matplotlib.pyplot as pl

p0=np.array([1.5,200.,2.,700.,3.,2400.])
hold=[0,1,0,1,0,0]
holdloc=[i for i, x in enumerate(hold) if x]
p1=np.delete(p0,holdloc)
timepoints = np.arange(-200.,4000.,20.)
callFitName= ut.genCallFitName(p0,hold)

y=ut.multiexp(timepoints,p0)+np.random.normal(0, .1, size=len(timepoints))

popt, pcov = curve_fit(ut.multiexp2, timepoints, y, p0=p0)
print 'unheld parameters:', popt

popt, pcov=eval(callFitName)
for i in range(0, len(hold)):
    if hold[i]:
        popt=np.insert(popt,i,p0[i])
print 'held parameters:', popt

yfit=ut.multiexp(timepoints,popt)

pl.plot(timepoints,y,timepoints,yfit)
pl.show()

And in utilities.py

import numpy as np

def multiexp2(t,*args):
    val=np.zeros(len(t))
    for i in np.arange(0,len(args),2):
        val+=args[i]*np.exp(-t/args[i+1])*np.heaviside(t,0)
    return val

def multiexp(t,args):
    val=np.zeros(len(t))
    for i in np.arange(0,len(args),2):
        val+=args[i]*np.exp(-t/args[i+1])*np.heaviside(t,0)
    return val

def genCallFitName(p0,hold):
    counter=0
    callFitName='curve_fit(lambda timepoints, *p: ('
    for i in range(0, len(hold)):
        if hold[i] and i%2==0:
            callFitName+=str(p0[i])
        elif hold[i] and i%2==1:
            callFitName+='*np.exp(-timepoints/'+str(p0[i])+')'
        elif i%2==0 and not hold[i]: 
            callFitName+='p['+str(counter)+']'
            counter+=1
        else:
            callFitName+='*np.exp(-timepoints/p['+str(counter)+'])'
            counter+=1
        if i<len(hold)-1 and i%2==1:
            callFitName+='+'
    callFitName+=')*np.heaviside(timepoints,0), timepoints, y, p0=p1)'
    return callFitName
Scott Ellis
  • 1
  • 1