Help with optimizing minimization/objective function
32 views
Skip to first unread message
Sarah Kathryn Loughran
Oct 23, 2023, 3:18:34 PM10/23/23
to lmfit-py
Hello all,
I am working on a complex minimization problem. I have an observational dataset of Doppler velocities across a planet, and I need to retrieve a latitudinal wind profile from it. to do this, I am generating a synthetic dataset from a first guess latitudinal profile. The residuals between the observations and the synthetic dataset are minimized y changing the latitudinal profile. ie. the latitudinal profile is the dependant parameter.
My current script is functional, but extremely slow. It was working for over 24 hours and did not get close to converging on a solution. Currently, the first guess profile is split into invidiual parameters for Minimizer(). I believe this is the source of the slow down, but because LMFIT requires parameters to be scalar, I am unsure of how to proceed. Reducing the number of parameters is not a final solution because a certian resolution is needed for the results to be meaningful.
How can I best optimize my code so that it produces a result in a reasonable time (a few hours runtime is expected)? The desired outputs are:
- final result array, 2D array - synthetic dataset that has minimized residuals compared to the observations
- fonal profile/fit parameters, 1D array - the optimized latitudinal wind profile that rsults in minimized residuals.
If LMFIT is not a good option for this problem, what other programs/packages/functions would you recommend instead?
Many thanks,
Sarah Loughran
Please note that the code below contains large sections that are required for the functions to run properly, but are not relevant to LMFIT. I did my best to comment extensively and make it clear what sections of the code are relevant, and which are not. Please let me know if you have any clarifying questions.
###############################################################################
import numpy as np
import lmfit
from scipy.interpolate import interp1d
from scipy.ndimage import convolve1d
from astropy.convolution import Gaussian2DKernel, convolve
from scipy.interpolate import RegularGridInterpolator
gauss = lambda x,x0,sigma: np.exp(-(x-x0)**2./(2.*sigma**2.)) # gaussian function
J = lambda nu,T: 6.6260755e-34*nu/(1.380658e-23*(np.exp(6.6260755e-34*nu/(1.380658e-23*T))-1)) # emission function
# Calculate latitudinal profile over all latitudes
def vfun(v,l):
global latitudes,planet_latitudes
vi=interp1d(latitudes,v,kind='linear',bounds_error=True)
sigma=15/2.35482
ggrid = np.arange(0,len(planet_latitudes))
psf = (1.0/sigma)*(2.0*np.pi)**-0.5*gauss(ggrid,len(planet_latitudes)/2,sigma)
vsmth=convolve1d(vi(planet_latitudes),psf,mode='constant',cval=0.0)
#Generate exponential taper
taper=1/(1.+np.exp(-(planet_latitudes-(5-90.))/1)) - 1/(1.+np.exp(-(planet_latitudes-(90.-5))/1))
vsmthi=interp1d(planet_latitudes,vsmth*taper,kind='linear',bounds_error=True)
return vsmthi(l)
# Important for LMFIT
# Calcualte synthetic dataset from a given latitudinal profile (the main function of the program)
# potential candidate for using lmfit.Model?
def getWind(v):
# Doppler shift of each cell.
dop = (d/rho) * vfun(v,lats) * np.cos(np.deg2rad(26.5))
# Set up model data cube
modelCube = np.zeros([len(vgrid),len(y_c),len(x_c)])
for i in range(len(x_c)):
for j in range(len(y_c)):
if ((x_c[i]**2+y_c[j]**2)**0.5 <= 2575 + 1200): # only consider cells within the atmosphere
# Set the background temperature
if ((x_c[i]**2+y_c[j]**2)**0.5 <= 2575):
Tbg = 94
else:
Tbg = 2.73
k=0
# apply mask to atmosperhic property arrays
mask = atmosShadowMask[i,j,:]
TBits=Tgrid[i,j,:][mask]
tauBits=tauCube[i,j,:][mask]
dopsigmaBits=dopsigmaCube[i,j,:][mask]
dopBits=dop[i,j,:][mask]
for v in vgrid:
# Get the Doppler-shifted atmoapsheric properties
tauBitsv=tauBits*gauss(v,dopBits,dopsigmaBits)
emission = 0.
continuum = Tbg
# Integrate the atmospheric properties along the line of sight
for tauBitv,TBit in zip(tauBitsv[::-1],TBits[::-1]):
rad = J(349446986900.0,TBit)-J(349446986900.0,continuum)
emission += rad * (1.-np.exp(-tauBitv))
continuum = emission + Tbg
modelCube[k,j,i] = emission
k += 1
# Convolve the model cube with beam size
modelCubec = np.zeros([len(vgrid),len(y_c),len(x_c)])
k=0
for plane in modelCube:
modelCubec[k] = convolve(plane, kernel)
k += 1
# calculate Weighted mean velocities
modelCubecv = modelCubec * vgrid[:,np.newaxis,np.newaxis]
vmodel = (np.ma.sum(modelCubecv,axis=0) / np.sum(modelCubec,axis=0)).transpose()
vmodelI=RegularGridInterpolator((x_c,y_c),vmodel,bounds_error=False, fill_value = None)
print('Run of getWind completed')
return vmodelI((akm,bkm))
# Important for LMFIT
# Objective function for preforming the fit
def getChisq(params,v_in,obsWind):
global xkm,ykm,akm,bkm,r #observed distances and winds
model=[]
# Parameters for this run in array format - input for getWind
v_alt = np.fromiter(params.values(), dtype=float)
print(v_alt)
# Get the velocity field model
model = getWind(v_alt)
# Get the normalized residuals
residuals = (v_obs - model)/((1.29*1000./2.)*(0.0114686+0.00858804)) # residuals multipled by error
return residuals
###############################################################################
# This section is necessary for the proper functioning of the example, but not LMFIT directly
latitudes = np.asarray([-90., -75., -60., -45, -30., -15., 0, 15., 30., 45, 60., 75., 90.])
velocities = np.asarray([0,50,120,150,160,170,180,150,120,100,70,50,0])
planet_latitudes=np.arange(-90,91)
# Define cartesian cube
cell_size = 100 # size of pixels in km
grid_size = 12000./100
# Orbital info
beam_maj = 0.274* 725.27 * 9.21 # Major axis of beam
beam_min = 0.253 * 725.27 * 9.21 # Minor axis of beam
bpa = -7.0 # Beam position/rotation angle
platescale = 0.025 #Pixel size of model in arcsec
kernel = Gaussian2DKernel(x_stddev=beam_maj/cell_size/2.3548,y_stddev=beam_min/cell_size/2.3548,theta=np.deg2rad(-bpa))
x_c = ((np.arange(grid_size)-grid_size/2.)+0.5)*cell_size
y_c = ((np.arange(grid_size)-grid_size/2.)+0.5)*cell_size
z_c = ((np.arange(grid_size)-grid_size/2.)+0.5)*cell_size
a_c,b_c,c_c = np.meshgrid(x_c,y_c,z_c)
c_radii = np.sqrt((a_c)**2+(b_c)**2+(c_c)**2)
#Rotate the coordinates
# Sky-projected distance of each point from the polar vector
d = (a_c - (-b_c * np.tan(np.deg2rad(5.3)))) * np.cos(np.deg2rad(5.3))
# Cylindrical polar radii
rho = (d**2 + c_c**2)**0.5
# Generate the latitude for each cell
# N pole vector
northx = -np.sin(np.deg2rad(5.3))*np.cos(np.deg2rad(26.5))
northy = np.cos(np.deg2rad(5.3))*np.cos(np.deg2rad(26.5))
northz = np.sin(np.deg2rad(26.5))
#dot product of north pole vector and cell vectors
dprod = (northx*a_c + northy*b_c + northz*-c_c)/c_radii
lats = 90. - np.degrees(np.arccos(dprod))
vgrid=np.arange(-1500,1600,100)
Tgrid = np.full([120,120,120],94)
dopsigmaCube = np.full([120,120,120],33.2754)
tauCube = np.full([120,120,120],1)
with np.errstate(invalid='ignore'):
atmosShadowMask=np.invert((c_radii>(2575+1200)) | (((a_c**2+b_c**2)**0.5<2575) & (c_c>-(2575**2-a_c**2-b_c**2)**0.5)))
# Final array
pix = np.arange(0,55)
obs_x = np.tile(pix,(55,1))
obs_y = np.tile(pix,(55,1))
obs_x = obs_x.transpose()
xkm = (obs_x[:,0]-len(obs_x[:,0])/2.+0.5)*platescale*9.21*725.27
ykm = (obs_y[0]-len(obs_y[0])/2.+0.5)*platescale*9.21*725.27
akm,bkm=np.meshgrid(xkm,ykm)
#plotting_mask = (((akm**2+bkm**2)**0.5)>4000)
###############################################################################
# Example observations
v_lat = np.arange(-90,90,3.3)
obs = [0.0757,61.2186,81.7521,104.7237,123.5139,154.2670,120.6526,125.6283,71.0417,63.7837,46.2720,48.1143,0.0772]
v_obs = getWind(obs)
#print(v_obs)
# Using LMFIT
# These parameters are the latitude and velocity pairs from the variables "latitudes" and "velocities" above
params = lmfit.Parameters()
params.add('v90n',value=0,vary=False)
params.add('v75n',value=50)
params.add('v60n',value=120)
params.add('v45n',value = 150)
params.add('v30n',value = 160)
params.add('v15n',value = 170)
params.add('v0',value = 180)
params.add('v15',value=150)
params.add('v30',value = 120)
params.add('v45',value = 100)
params.add('v60',value=70)
params.add('v75',value = 50)
params.add('v90',value = 0, vary=False)
fit = lmfit.Minimizer(getChisq, params,fcn_args=(velocities,v_obs))
print('Starting LMFIT...')
result = fit.minimize()
final_profile = np.fromiter(params.values(), dtype=float)
print(result)
print(final_profile)
import lmfit
from scipy.interpolate import interp1d
from scipy.ndimage import convolve1d
from astropy.convolution import Gaussian2DKernel, convolve
from scipy.interpolate import RegularGridInterpolator
gauss = lambda x,x0,sigma: np.exp(-(x-x0)**2./(2.*sigma**2.)) # gaussian function
J = lambda nu,T: 6.6260755e-34*nu/(1.380658e-23*(np.exp(6.6260755e-34*nu/(1.380658e-23*T))-1)) # emission function
# Calculate latitudinal profile over all latitudes
def vfun(v,l):
global latitudes,planet_latitudes
vi=interp1d(latitudes,v,kind='linear',bounds_error=True)
sigma=15/2.35482
ggrid = np.arange(0,len(planet_latitudes))
psf = (1.0/sigma)*(2.0*np.pi)**-0.5*gauss(ggrid,len(planet_latitudes)/2,sigma)
vsmth=convolve1d(vi(planet_latitudes),psf,mode='constant',cval=0.0)
#Generate exponential taper
taper=1/(1.+np.exp(-(planet_latitudes-(5-90.))/1)) - 1/(1.+np.exp(-(planet_latitudes-(90.-5))/1))
vsmthi=interp1d(planet_latitudes,vsmth*taper,kind='linear',bounds_error=True)
return vsmthi(l)
# Important for LMFIT
# Calcualte synthetic dataset from a given latitudinal profile (the main function of the program)
# potential candidate for using lmfit.Model?
def getWind(v):
# Doppler shift of each cell.
dop = (d/rho) * vfun(v,lats) * np.cos(np.deg2rad(26.5))
# Set up model data cube
modelCube = np.zeros([len(vgrid),len(y_c),len(x_c)])
for i in range(len(x_c)):
for j in range(len(y_c)):
if ((x_c[i]**2+y_c[j]**2)**0.5 <= 2575 + 1200): # only consider cells within the atmosphere
# Set the background temperature
if ((x_c[i]**2+y_c[j]**2)**0.5 <= 2575):
Tbg = 94
else:
Tbg = 2.73
k=0
# apply mask to atmosperhic property arrays
mask = atmosShadowMask[i,j,:]
TBits=Tgrid[i,j,:][mask]
tauBits=tauCube[i,j,:][mask]
dopsigmaBits=dopsigmaCube[i,j,:][mask]
dopBits=dop[i,j,:][mask]
for v in vgrid:
# Get the Doppler-shifted atmoapsheric properties
tauBitsv=tauBits*gauss(v,dopBits,dopsigmaBits)
emission = 0.
continuum = Tbg
# Integrate the atmospheric properties along the line of sight
for tauBitv,TBit in zip(tauBitsv[::-1],TBits[::-1]):
rad = J(349446986900.0,TBit)-J(349446986900.0,continuum)
emission += rad * (1.-np.exp(-tauBitv))
continuum = emission + Tbg
modelCube[k,j,i] = emission
k += 1
# Convolve the model cube with beam size
modelCubec = np.zeros([len(vgrid),len(y_c),len(x_c)])
k=0
for plane in modelCube:
modelCubec[k] = convolve(plane, kernel)
k += 1
# calculate Weighted mean velocities
modelCubecv = modelCubec * vgrid[:,np.newaxis,np.newaxis]
vmodel = (np.ma.sum(modelCubecv,axis=0) / np.sum(modelCubec,axis=0)).transpose()
vmodelI=RegularGridInterpolator((x_c,y_c),vmodel,bounds_error=False, fill_value = None)
print('Run of getWind completed')
return vmodelI((akm,bkm))
# Important for LMFIT
# Objective function for preforming the fit
def getChisq(params,v_in,obsWind):
global xkm,ykm,akm,bkm,r #observed distances and winds
model=[]
# Parameters for this run in array format - input for getWind
v_alt = np.fromiter(params.values(), dtype=float)
print(v_alt)
# Get the velocity field model
model = getWind(v_alt)
# Get the normalized residuals
residuals = (v_obs - model)/((1.29*1000./2.)*(0.0114686+0.00858804)) # residuals multipled by error
return residuals
###############################################################################
# This section is necessary for the proper functioning of the example, but not LMFIT directly
latitudes = np.asarray([-90., -75., -60., -45, -30., -15., 0, 15., 30., 45, 60., 75., 90.])
velocities = np.asarray([0,50,120,150,160,170,180,150,120,100,70,50,0])
planet_latitudes=np.arange(-90,91)
# Define cartesian cube
cell_size = 100 # size of pixels in km
grid_size = 12000./100
# Orbital info
beam_maj = 0.274* 725.27 * 9.21 # Major axis of beam
beam_min = 0.253 * 725.27 * 9.21 # Minor axis of beam
bpa = -7.0 # Beam position/rotation angle
platescale = 0.025 #Pixel size of model in arcsec
kernel = Gaussian2DKernel(x_stddev=beam_maj/cell_size/2.3548,y_stddev=beam_min/cell_size/2.3548,theta=np.deg2rad(-bpa))
x_c = ((np.arange(grid_size)-grid_size/2.)+0.5)*cell_size
y_c = ((np.arange(grid_size)-grid_size/2.)+0.5)*cell_size
z_c = ((np.arange(grid_size)-grid_size/2.)+0.5)*cell_size
a_c,b_c,c_c = np.meshgrid(x_c,y_c,z_c)
c_radii = np.sqrt((a_c)**2+(b_c)**2+(c_c)**2)
#Rotate the coordinates
# Sky-projected distance of each point from the polar vector
d = (a_c - (-b_c * np.tan(np.deg2rad(5.3)))) * np.cos(np.deg2rad(5.3))
# Cylindrical polar radii
rho = (d**2 + c_c**2)**0.5
# Generate the latitude for each cell
# N pole vector
northx = -np.sin(np.deg2rad(5.3))*np.cos(np.deg2rad(26.5))
northy = np.cos(np.deg2rad(5.3))*np.cos(np.deg2rad(26.5))
northz = np.sin(np.deg2rad(26.5))
#dot product of north pole vector and cell vectors
dprod = (northx*a_c + northy*b_c + northz*-c_c)/c_radii
lats = 90. - np.degrees(np.arccos(dprod))
vgrid=np.arange(-1500,1600,100)
Tgrid = np.full([120,120,120],94)
dopsigmaCube = np.full([120,120,120],33.2754)
tauCube = np.full([120,120,120],1)
with np.errstate(invalid='ignore'):
atmosShadowMask=np.invert((c_radii>(2575+1200)) | (((a_c**2+b_c**2)**0.5<2575) & (c_c>-(2575**2-a_c**2-b_c**2)**0.5)))
# Final array
pix = np.arange(0,55)
obs_x = np.tile(pix,(55,1))
obs_y = np.tile(pix,(55,1))
obs_x = obs_x.transpose()
xkm = (obs_x[:,0]-len(obs_x[:,0])/2.+0.5)*platescale*9.21*725.27
ykm = (obs_y[0]-len(obs_y[0])/2.+0.5)*platescale*9.21*725.27
akm,bkm=np.meshgrid(xkm,ykm)
#plotting_mask = (((akm**2+bkm**2)**0.5)>4000)
###############################################################################
# Example observations
v_lat = np.arange(-90,90,3.3)
obs = [0.0757,61.2186,81.7521,104.7237,123.5139,154.2670,120.6526,125.6283,71.0417,63.7837,46.2720,48.1143,0.0772]
v_obs = getWind(obs)
#print(v_obs)
# Using LMFIT
# These parameters are the latitude and velocity pairs from the variables "latitudes" and "velocities" above
params = lmfit.Parameters()
params.add('v90n',value=0,vary=False)
params.add('v75n',value=50)
params.add('v60n',value=120)
params.add('v45n',value = 150)
params.add('v30n',value = 160)
params.add('v15n',value = 170)
params.add('v0',value = 180)
params.add('v15',value=150)
params.add('v30',value = 120)
params.add('v45',value = 100)
params.add('v60',value=70)
params.add('v75',value = 50)
params.add('v90',value = 0, vary=False)
fit = lmfit.Minimizer(getChisq, params,fcn_args=(velocities,v_obs))
print('Starting LMFIT...')
result = fit.minimize()
final_profile = np.fromiter(params.values(), dtype=float)
print(result)
print(final_profile)
Matt Newville
Oct 25, 2023, 5:03:28 PM10/25/23
to lmfi...@googlegroups.com
Hi Sarah,
Sorry for the delay on this. I'm not sure I have a definitive answer. It doesn't seem like your array sizes or number of parameters is *that* large for a runtime to be that slow, but it does seem that your `getWind` function - which is called at each iteration by the objective function - has a triple loop in it. I think that might be the most obvious place to look for "slow code". If that could be "vectorized" in any way, I would bet that could give a speedup, but I have to admit that nothing is really jumping out at me as "oh, that's the slow part that can be trivially made into a numpy ufunc call". Perhaps the masks could be precalculated or something?
--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/b51c0e34-014d-4288-b7d8-b7a805602486n%40googlegroups.com.
--Matt
Sarah Loughran
Oct 26, 2023, 10:34:11 AM10/26/23
to lmfi...@googlegroups.com
While I plan on editing the getWind function in future, I am confident that is not the main problem. It takes less than 30 seconds to run in its current state.
I believe how I have the minimization set up (11 free/variable parameters) is the problem. I just tried having all but one parameter held fixed (v0 was variable), and the script finished in just under 8 minutes. The final result was the same as the input, which is a different problem.
Is there a better way to set up the parameters? When I have the parameters of each run of getWind print to the screen, it shows that only one parameter is altered at a time and by a small amount (<1). With 11 variable parameters, coming to a solution in a reasonable time frame is impossible.
Would it be better to set up the minimization as a Model()? I initially attempted this with getWind, but I could not get it to work as I wanted; the final wind array would be altered/minimized, but the parameters would remain the same.
Sarah
You received this message because you are subscribed to a topic in the Google Groups "lmfit-py" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/lmfit-py/IYoj4jGcF4c/unsubscribe.
To unsubscribe from this group and all its topics, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESboC5g3wa-z%3D%3DL8VDFPv4JhE-OH797%3DJiXQP7f56Uwqefg%40mail.gmail.com.
Reply all
Reply to author
Forward
0 new messages