CODEWITHSUNDEEP
pythonoptimizationscipy
Navaneeth Mohan
Navaneeth Mohan

Reputation: 41

scipy optimize fmin syntax

numseq = ['0012000', '0112000', '0212000', '0312000', '1012000', '1112000',                                                                                   '1212000', '1312000', '2012000', '2112000', '2212000', '2312000', '3012000', '3112000',          '3212000', '3312000', '0002000', '0022000', '0032000', '1002000', '1022000', '1032000',     '2002000', '2022000', '2032000', '3002000', '3022000', '3032000', '0010000', '0011000', '0013000', '1010000', '1011000', '1013000', '2010000', '2011000', '2013000', '3010000', '3011000', '3013000', '0012100', '0012200', '0012300', '1012100', '1012200', '1012300', '2012100', '2012200', '2012300', '3012100']
prob = [-0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.78361598908750163, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212]

numseq and prob are lists of length 50 each. They are the experimental data that is collected. numseq corresponds to the X axis values, and prob corresponds to the Y axis values.

The function that I want to minimise is:

def residue(allparams, xdata, ydata):
    chi2 = 0.0
    for i in range(0,len(xdata)):
        x = xdata[i]
        y = 0
        for j in range(len(x)):
            y = y-allparams[int(x[j])][j]
            chi2 = chi2 + (ydata[i]-y)**2
return chi2

So:

chi2 is the squared difference between the experimental and model values. This is what has to be minimised.

The initial guess for the parameters is given by:

x0 = [[-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6]]

Now how do I call fmin on this function? I tried 

fmin(residue, x0, args=(numseq, prob))

but I keep getting an error:

Traceback (most recent call last):
  File "<pyshell#362>", line 1, in <module>
    fmin(residue, x0, args=(numseq, prob))
  File "C:\Python31\lib\site-packages\scipy\optimize\optimize.py", line 258, in fmin
    fsim[0] = func(x0)
  File "C:\Python31\lib\site-packages\scipy\optimize\optimize.py", line 177, in function_wrapper
    return function(x, *args)
  File "<pyshell#361>", line 7, in residue
    y = y-allparams[int(x[j])][j]
IndexError: invalid index to scalar variable.

why is this so? Is it because fmin can't accept 2D arrays as initial guesses? Then do i have to change my entire code to work on a 1D array of parameters?

Even if you can't explain this problem, could you at least tell me how exactly the fmin module works? i.e the syntax of how to implement fmin for the optimization of an N-dimensional array? Could you explain what args() is? I'm new to optimisation and I have no idea on how to implement it :(

Upvotes: 3

Views: 2965

Answers (1)

Joel Vroom
Joel Vroom

Reputation: 1691

The "fmin" routine can accept 2d arrays as an initial guess. But the first thing it does is flatten this array [(4,7) --> (28)]. So what happens is that your residue function takes a (4,7) array as an input and the "fmin" routine gives it a flattened "x0" of length 28. That's why you see the error:
y = y-allparams[int(x[j])][j]
IndexError: invalid index to scalar variable.

See the source code here.

So it appears that you will have to change your residue function to accept a vector instead of an array. This doesn't seem too bad however. I tried the following which seemed to work (Note: Please double check!)

def residue_alternative(allparams, inshape, xdata, ydata):
    m, n = inshape
    chi2 = 0.0
    for i in range(0,len(xdata)):
        x = xdata[i]
        y = 0
        for j in range(len(x)):
            idx = int(x[j]) * n +  j #Double check this to 
            y = y-allparams[idx]     #make sure it does what you want
            chi2 = chi2 + (ydata[i]-y)**2
    return chi2

I called it using:

x0 = -0.6 * np.ones((4,7), dtype=np.double)
[xopt, fopt, iter, funcalls, warnflag] = \
    fmin(residue_alternative, x0, args=(x0.shape, numseq, prob),
       maxiter = 100000,
       maxfun  = 100000,
       full_output=True,
       disp=True)

And received the following results:

Optimization terminated successfully.
         Current function value: 7.750523
         Iterations: 21570
         Function evaluations: 26076

>>>xopt
array([ 0.57669042, -0.21965861,  0.2635061 , -0.08284016, -0.0779489 ,
   -0.10358114,  0.14041582,  0.72469391, -0.43190214,  0.31269757,
   -0.0338726 , -0.14919739, -2.58314651,  2.74251214,  0.57695759,
   -0.49574628,  0.1490926 ,  0.04912353,  0.02420988,  1.17924051,
   -7.2147027 ,  0.57860843, -0.28386938,  0.2431877 , -0.22674694,
   -0.58308225, -6.05706775, -2.06350063])

Which you can reshape to a 4x7 array. Give this a try and let me know if it works/helps.

Upvotes: 2

Related Questions