# Constrained Model fitting

## Curve fitting - Constraints¶

#### 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.¶

In [15]:
from scipy.optimize import minimize
import scipy
import math
from scipy import special
import pyfits

In [2]:
## 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

In [3]:
## 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)

In [4]:
## 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]))

In [5]:
import numpy as np
import matplotlib.pyplot as plt

In [6]:
path = "/home/fatima/Desktop/project_3/"
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.¶

In [8]:
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
[ 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:¶

In [9]:
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.¶

In [11]:
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.¶

In [13]:
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
[ 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



#### 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).¶

In [19]:
## 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
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)

In [20]:
## 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


#### Now we construct the cons object:¶

In [21]:
cons = ({'type':'eq','fun':Con})

In [ ]:
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¶

In [ ]:
chi_sq = Minimize(mini.x,x,y,weights)