Curve fitting - Constraints

Sometimes we need to put some constraints on the problem, meaning that we would like to have a set of best-fit parameters which correspond to some wanted conditions, usually more sophisticated than just introducing bounds.

Let's start with defining the objective function and Jacobian matrix. The fitting function we will use is scipy.optimize.minimize. It is the only function that I found which takes complicated constraints.

from scipy.optimize import minimize
import scipy
import math
from scipy import special
import pyfits

## the model

def Erfc(x,sigma):

 y = special.erfc(x/(sigma*np.sqrt(2)))
 return y


def SL_fit(x,w1,w2,w3,s1,s2,s3):

  f = 0.5*(w1*Erfc(x,s1)+w2*Erfc(x,s2)+w3*Erfc(x,s3)+ (1-w1-w2-w3))
  return f

## the customised objective function for the minimize function.
def Minimize(params,x,y,w):
 
 w1 = params[0]
 w2 = params[1]
 w3 = params[2]
 s1 = params[3]
 s2 = params[4]
 s3 = params[5]
 model = 0.5*(w1*Erfc(x,s1)+w2*Erfc(x,s2)+w3*Erfc(x,s3)+(1-w1-w2-w3))
 return np.sum(((1/w)*(model-y))**2)

## The Jacobian matrix
def Jac(popt,x,y,w):
 w1 = popt[0]
 w2 = popt[1]
 w3 = popt[2]
 s1 = popt[3]
 s2 = popt[4]
 s3 = popt[5]

 dfw1 = np.sum(2*(1/w**2)*(SL_fit(x,w1,w2,w3,s1,s2,s3)-y)*(0.5*(Erfc(x,s1)-1)))
 dfw2 = np.sum(2*(1/w**2)*(SL_fit(x,w1,w2,w3,s1,s2,s3)-y)*(0.5*(Erfc(x,s2)-1)))
 dfw3 =  np.sum(2*(1/w**2)*(SL_fit(x,w1,w2,w3,s1,s2,s3)-y)*(0.5*(Erfc(x,s3)-1)))
 dfs1 = np.sum(2*(1/w**2)*(SL_fit(x,w1,w2,w3,s1,s2,s3)-y)*(w1*x*np.exp(-x**2/(2*s1**2))/ (np.sqrt(2*np.pi)*s1**2)))
 dfs2 =np.sum( 2*(1/w**2)*(SL_fit(x,w1,w2,w3,s1,s2,s3)-y)*(w2*x*np.exp(-x**2/(2*s2**2))/ (np.sqrt(2*np.pi)*s2**2)))
 dfs3 = np.sum(2*(1/w**2)*(SL_fit(x,w1,w2,w3,s1,s2,s3)-y)*(w3*x*np.exp(-x**2/(2*s3**2))/ (np.sqrt(2*np.pi)*s3**2)))
 return np.transpose(np.array([dfw1,dfw2,dfw3,dfs1,dfs2,dfs3]))

import numpy as np
import matplotlib.pyplot as plt

path = "/home/fatima/Desktop/project_3/"
file = np.loadtxt(path+'limb_profile_av_norm_shifted')
x = file[:,0]
y = file[:,1]
ind = np.where(x>=0)
x = x[ind]
y = y[ind]
weights = np.sqrt(np.abs(y))
#p0=[0.3, 0.3, 0.2, 1, 2, 3]
p0  = [0.544,0.291,0.068,0.3,0.9,4.7]

For unconstrained fit, the minimize function is called. Note that the args parameter of the function should have the same format as the parameters passed to the objective function.

mini = scipy.optimize.minimize(Minimize,p0,args=(x,y,weights),jac=Jac,options={'xtol':1e-8,'gtol': 1e-8,'maxit':5000, 'disp': True})
print mini.x

Optimization terminated successfully.
         Current function value: 0.001221
         Iterations: 40
         Function evaluations: 51
         Gradient evaluations: 51
[ 0.52750273  0.28882469  0.10191702  0.25905394  0.76540951  2.83344876]

/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:1: OptimizeWarning: Unknown solver options: xtol, maxit
  """Entry point for launching an IPython kernel.

Now let's introduce some easy-defined constraints. Usually the defined constraint function which will be passed later to the fitting function minimize, returns a zero or a positive value. For example if we want to have the parameter s2 to be larger than s1 in the resultant fitting process, we can type the following:

def Con_1(params):
 w1 = params[0]
 w2 = params[1]
 w3 = params[2]
 w4 = 1-(w1+w2+w3)
 s1 = params[3]
 s2 = params[4]
 s3 = params[5] 
 return s2-s1    

And then we add this constraint function to the object cons and specify whether it is an equality or inequality constraint. The reason for defining this object is to include as many constraint functions as possible. This object is passed to the function minimize under the parameter constraints.

cons = ({'type':'ineq','fun':Con_1}) ##this tells minimize that we want the returned value of the constraint function should be positive

Now we pass this constraint to minimize but we have to also change the method of optimization which is defined for constrained problems. Check scipy.optimize.minimize for more details on the optimization methods available.

mini = scipy.optimize.minimize(Minimize,p0,args=(x,y,weights),jac=Jac,constraints=cons,method='SLSQP',options={'xtol':1e-8,'gtol': 1e-8,'maxit':5000, 'disp': True})
popt = mini.x ## the best-fit parameters
print popt

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.00233234740213
            Iterations: 8
            Function evaluations: 17
            Gradient evaluations: 8
[ 0.56554532  0.27822938  0.09418257  0.27118045  0.8841216   4.69687872]

/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:2: OptimizeWarning: Unknown solver options: xtol, maxit, gtol
  

Now I will introduce a more complicated constraint function. For a better understanding of this function, I recommend going back to data_cleaning which explaines how to use the model introduced here to cleaning my data for stray light.

In my project, I have used this constraint function to force minimize to fit the model to the data in a way that the RMS contrast (check RMS_contrast) of the corrected image should be equal to a certain value, usually that returned by theoretical expectations (Magnetohydrodynamical simulations of the solar convection).

## Here I am loading the image to be corrected and embedding it to preparing it for use by the constraint function

image = 'sufi300_lev3_2009.fits'
imscale=0.0198
Input = pyfits.getdata(path+image, ignore_missing_header=True)
w_in = Input.shape[1]
h_in = Input.shape[0]
Input_emb = np.zeros((h_in,h_in))
Input_emb[:,(h_in-w_in)/2:(h_in+w_in)/2] = Input  
mtf_pd = pyfits.getdata(path+'MTF_PD_large.fits')[0]
freqscale = 1./(h_in*imscale)
otf = np.zeros((h_in,h_in))
i,j = np.indices([h_in,h_in])*freqscale
center = np.array([(i.max()-i.min())/2.0, (j.max()-j.min())/2.0])
xc = int(center[1])
yc = int(center[0])
r = np.hypot(i - yc, j - xc)

## The constraint function
def Con(params):
 w1 = params[0]
 w2 = params[1]
 w3 = params[2]
 w4 = 1-(w1+w2+w3)
 s1 = params[3]
 s2 = params[4]
 s3 = params[5]
 s_t = s1
 ss2 = np.sqrt(s2**2 - s_t**2)
 ss3 = np.sqrt(s3**2 - s_t**2)
 otf = w1*mtf_pd+w2*np.exp(-2 * np.pi**2 * ss2**2* r**2 )+w3*np.exp(-2 * np.pi**2 *ss3**2 * r**2 )
 otf[yc,xc] = otf[yc,xc]+w4
 ind_ne_0 = np.where(mtf_pd!=0)
 ind_eq_0 = np.where(mtf_pd==0)
 otf[ind_ne_0] = otf[ind_ne_0]/mtf_pd[ind_ne_0]
 otf[ind_eq_0] = w1
 otf = np.fft.ifftshift(otf)
 Input_fft = np.fft.fft2(Input_emb)
 Input_deconv = np.real(np.fft.ifft2(Input_fft/otf))
 Input_deconv = Input_deconv[:,(h_in-w_in)/2:(h_in+w_in)/2]
 RMS = Input_deconv.std()/Input_deconv.mean()
 return 0.3-RMS

As you can notice, everytime minimize performs the least-square fitting, it should returns the computed parameters to the constaint function, Con to build the total PSF, decovolve the science image with it and test the returned RMS value of the corrected image if it's equal to 0.3 or not, so it is computationally expensive.

Now we construct the cons object:

cons = ({'type':'eq','fun':Con})

mini = minimize(Minimize,p0,args=(x,y,weights),constraints=cons,method='SLSQP',options={'xtol':1e-8,'gtol': 1e-8,'maxit':5000, 'disp': True})

Since the returned value by the objective function Minimize is the chi-squared value (check the expression of the objective function), we can compute the chi-squared value like this

chi_sq = Minimize(mini.x,x,y,weights)