lmfit -py using arrays for parameter optimization

Question

Situation:
I'm trying to optimise parameters for a natural creek where gases degas or ingas at a certain rate according to reasonably well established equations. We have the measured concentrations at certain distances downstream and would like to use optimization techniques to determine the value of some unknown parameters in the model.

I'm trying to use lmfit lmfit-py (github) to optimize parameters using the code pasted below.

• I would like some of the unknown parameters to be allowed to have a different value for each different sample point along the stream (this reflects the actual situation). In the script these are termed 'local parameters'. As such, the results for these parameters should be a list for each location along the stream.

• The other parameters should be optimized to have the same value for all locations along the stream, termed 'global parameters'. As such the optimization results for these parameters should be a single value.

I'm using Python 3.2.

Problems:

1. At the moment I'm getting only single value results for all parameters.
2. Can lmfit be used to produce optimized arrays?
3. I really think I've not set the script up properly and is probably not calling the bounds and initial guesses correctly or something???

Script I'm Using:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import minimize, Parameters, report_fit
import scipy as sp


# import data to be fitted
x,dc_dx,E,lmbd,c,h,th,theta,d,w,c_min,\
    c_max,Q_min,Q_max,w_min,w_max,d_min,d_max,theta_min,theta_max,\
    I_min,I_max,k_min,k_max,h_min,h_max,th_min,th_max,c_ini,Q_ini,\
    w_ini,d_ini,theta_ini,I_ini,k_ini,h_ini,th_ini,cigl_min,\
    cigl_max,gammagl_min,gammagl_max,\
            cigl_ini,gammagl_ini,kgl_min,kgl_max,kgl_ini,\
            temp,scond,econd,disol,O2per,O2,pH,ORP,\
            c_atmdef,dO2_dx,kO2_ini,kO2_min,kO2_max,O2i=\
        np.genfromtxt("Rn2011DryInput_1.csv",unpack=True,\
                        skip_header =2,delimiter=",")

'''
'Global parameters' are those which when optimised result in the
same value along the whole transect.
'Local parameters' are those parameters optimised where the value can
vary at each sampling site.
(global parameters should be packed at the beginning)
              for globals unpack each one e.g.   p1=par[0]
                                                 p2 = par[1]

'''

# Define smoothing error
def ratio(a,b):
    error = 0.0
    if b > sys.float_info.epsilon:
        error = ((a-b)/b)*2
    else:
        if a > sys.float_info.epsilon:
            error = ((a-b)/a)**2

# Try reducing the error for the regularisation
    error = 0.5*error
    return(error)


# Define model
def radon_f(par, *args):
    nargs = 12  #  number of arguments
    npl = 4     #  number of local parameters
    npg = 2     #  number of global parametrs
    ci=par['ci'].value
    gamma=par['gamma'].value

    n = int(len(par)/npl)     
    error = 0.0
    for i in range(n):
        Q=par['Q'].value
        I=par['I'].value
        k=par['k'].value
        kO2=par['kO2'].value
        dc_dx=args[i*nargs]
        E=args[i*nargs+1]
        lmbd=args[i*nargs+2]
        c=args[i*nargs+3]
        h=args[i*nargs+4]
        th=args[i*nargs+5]
        theta=args[i*nargs+6]
        d=args[i*nargs+7]
        w=args[i*nargs+8]
        O2i=args[i*nargs+9]
        O2=args[i*nargs+10]
        c_atmdef=args[i*nargs+10]


        measurements = dc_dx
        model=(I*(ci-c)+w*E*c-k*w*c-d*w*lmbd*c+\
           gamma*h*w*theta/(1+lmbd*th)-lmbd*h*w*theta*c/(1+lmbd*th))/Q

        error= error + ((measurements - model)/measurements)**2

        measurements=dO2_dx
        model = (I*(O2i-O2)+w*E*O2-kO2*w*c_atmdef)/Q

        error= error + ((measurements - model)/measurements)**2

        if i > 0:    # apply regularization
            Qp,Ip,kp,kO2p=par[(npg+npl*(i-1)):(npg+npl*i)]
            error = error + ratio(Q,Qp)
            error = error + ratio(I,Ip)
            error = error + ratio(k,kp)
            error = error + ratio(kO2,kO2p)
    return(error)


# create a set of parameters
par = Parameters()
for i in range(int(len(x))):    
    par.add('Q',   value=Q_ini[i],  vary=True,  min= Q_min[i],    max=Q_max[i])
    par.add('I',   value=I_ini[i],     vary=True,  min= I_min[i], max=I_max[i])
    par.add('k',   value=k_ini[i],      vary=True,  min= 0.1,     max=6.2    )
    par.add('ci',  value=cigl_ini[i],  vary=False, min= 6000,     max=25000  )
    par.add('gamma',value=gammagl_ini[i],   vary=False,  min=200,   max=3000 )
    par.add('kO2',  value=kO2_ini[i],   vary=True, min=2.5,       max=10     )

# Define arguements
args=[]
for i in range(len(E)):
    args.append(dc_dx[i])
    args.append(E[i])
    args.append(lmbd[i])
    args.append(c[i])
    args.append(h[i])
    args.append(th[i])
    args.append(theta[i])
    args.append(d[i])
    args.append(w[i])
    args.append(O2i[i])

# run the fit
result = minimize(radon_f, par, args=args)

print("all results")
print(result)

Q=par['Q'].value
I=par['I'].value
k=par['k'].value
kO2=par['kO2'].value

print(Q)
print(I)
print(k)
print(kO2)


# print results of global parameters
print('ci', 'gamma')
print(result[0][0], result[0][3])

for i in range(len(c_min)):
    # Print out the results of local parameters
     print(result[0][npg+npl*i],result[0][npg+npl*i+1],\
           result[0][npg+npl*i+2],'\t')

Data I'm using:

0   1   2   3   4   5   6   7   8   9   10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59
x   dc_dx   E   lmbd    c   h   th  theta   d   w   c_min   c_max   Q_min   Q_max   w_min   w_max   d_min   d_max   theta_min   theta_max   I_min   I_max   k_min   k_max   h_min   h_max   th_min  th_max  c_ini   Q_ini   w_ini   d_ini   theta_ini   I_ini   k_ini   h_ini   th_ini  cigl_min    cigl_max    gammagl_min gammagl_max cigl_ini    gammagl_ini kgl_min kgl_max kgl_ini temp    scond   econd   disol   O2per   O2  pH  ORP c_atmdef    dO2_dx  kO2_ini kO2_min kO2_max O2i
514 -6.763285345    0.0045  0.181   9572    0.7 0.5 0.3 0.5 12.2    8614.762109 10529.15369 0   14200   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   9572    14000   3   0.4 0.3 50  2   0.6 0.6 6000    25000   200 3000    7000    700 1   7   2   26.93   752.3   780.1   489 43  3.42    7.27    8   -267.48 0.012840237 5   2.5 10  1.3
683 -5.490998291    0.0045  0.181   8429    0.7 0.5 0.3 0.5 12.2    7586.066408 9271.858943 0   14200   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   8429    14000   3   0.4 0.3 50  2   0.6 0.6                                     26.61   750.4   773.4   487.7   69.8    5.59    7.33    12  -265.31 0.005837412 5   2.5 10  1.3
949 -4.22793979 0.0045  0.181   7307    0.7 0.5 0.3 0.5 12.2    6576.106938 8037.464035 0   14200   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   7307    14000   3   0.4 0.3 50  2   0.6 0.6                                     26.3    750.7   769.4   488 65.6    5.28    7.38    26  -265.62 0.000472201 5   2.5 10  1.3
1287    -2.085993575    0.0045  0.181   5875    0.7 0.5 0.3 0.5 12.2    5287.160328 6462.084846 0   14200   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   5875    14000   3   0.4 0.3 50  2   0.6 0.6                                     26.15   750 766 487 71.6    5.78    7.47    19.5    -265.12 0.00183869  5   2.5 10  1.3
1623    -1.048348565    0.0045  0.181   5897    0.7 0.5 0.3 0.5 12.2    5306.871121 6486.175814 12822.3 15671.7 3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   5897    14247   3   0.4 0.3 50  2   0.6 0.6                                     25.98   748.1   762.2   486.3   80.5    6.52    7.64    27  -264.38 0.001575445 5   2.5 10  1.3
1992    3.150847718 0.0045  0.181   5099    0.7 0.5 0.3 0.5 12.2    4588.91133  5608.669404 14500   16500   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   5099    16707   8   1   0.3 50  2   0.6 0.6                                     26.03   750 765 487 85  6.87    7.56    26  -264.03 -0.003467278    5   2.5 10  1.3
2488    17.92239204 0.0045  0.181   9297    0.7 0.5 0.3 0.5 12.2    8367.050656 10226.39525 35500   59500   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   9297    49279   8   1   0.3 100 2   0.6 0.6                                     29.06   726 783 472 38.5    2.96    7.14    -83.4   -267.94 -0.010477451    5   2.5 10  1.3
2569    10.03067057 0.0045  0.181   11515   0.7 0.5 0.3 0.5 12.2    10363.14089 12666.06109 44500   60000   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   11515   59739   8   1   0.3 750 2   0.6 0.6                                     29.5    730 793 474 43  3.28    7.07    -1.3    -267.62 0.002783579 5   2.5 10  1.3
2835    -6.606394762    0.0045  0.181   9568    0.7 0.5 0.3 0.5 12.2    8610.764203 10524.26736 44500   60000   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   9568    58576   8   1   0.3 250 2   0.6 0.6                                     29.3    727 787 472 48  3.71    7.12    11.1    -267.19 0.001373823 5   2.5 10  1.3
3224    -4.694876493    0.0045  0.181   7275    0.7 0.5 0.3 0.5 12.2    6547.652797 8002.686751 44500   60000   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   7275    56874   8   1   0.3 150 2   0.6 0.6                                     28.29   729.9   775.7   474.4   53.4    4.15    7.27    21  -266.75 0.001830971 5   2.5 10  1.3
3609    -4.654680461    0.0045  0.181   5929    0.7 0.5 0.3 0.5 12.2    5336.000281 6521.778121 44500   60000   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   5929    55190   10  1   0.3 150 2   0.6 0.6                                     29.403142   734 702 477 59.6    5.13    7.22    -28 -265.77 0.001991543 5   2.5 10  1.3
4082    -3.263752353    0.0045  0.181   3180    0.7 0.5 0.3 0.5 12.2    2861.606999 3497.519665 44500   60000   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   3180    53121   10  1.5 0.3 150 2   0.6 0.6                                     28.164949   729 681 474 66  5.81    7.5 -33 -265.09 0.001000506 5   2.5 10  1.3
5218    1.167013246 0.0045  0.181   2367    0.7 0.5 0.3 0.5 12.2    2130.615084 2604.085102 44500   60000   3   15  0.1 2   0.2 0.45    0   1000    1   6.17    0.4 1   0   2   2367    48152   10  1.5 0.3 50  2   0.6 0.6                                     26.425339   684 616 444 70.8    6.45    7.83    -41 -264.45 0.00048454  5   2.5 10  1.3
5506    1.045825103 0.0045  0.181   3245    0.7 0.5 0.3 0.5 12.2    2920.916644 3570.009232 44500   60000   3   15  0.1 2   0.2 0.45    0   1000    0   6.17    0.4 1   0   2   3245    46893   10  1.5 0.3 150 2   0.6 0.6                                     26.527669   681 615 443 75.4    6.85    7.83    -52 -264.05 0.000191651 5   2.5 10  1.3
6043    -0.741605916    0.0045  0.181   2731    0.7 0.5 0.3 0.5 12.2    2458.228071 3004.500975 40089.6 48998.4 3   15  0.1 2   0.2 0.45    0   1000    0   6.17    0.4 1   0   2   2731    44544   10  1.5 0.3 100 2   0.6 0.6                                     24.900622   690 602 448 67.2    6.31    7.75    -38 -264.59 0.000307351 5   2.5 10  1.3
7851    -0.366090159    0.0045  0.181   1781    0.7 0.5 0.3 0.5 12.2    1602.550136 1958.672388 44500   53500   3   15  0.1 2   0.2 0.45    0   1000    0   6.17    0.4 1   0   2   1781    46298   10  1.5 0.3 50  2   0.6 0.6                                     24.675496   679 590 441 71.8    6.77    7.85    -26 -264.13 0.000430647 5   2.5 10  1.3
15277   -0.206321214    0.0045  0.181   248 0.7 0.5 0.3 0.5 12.2    223.6229375 273.3169236 48150   58850   3   15  0.1 2   0.2 0.45    0   1000    0   6.17    0.4 1   0   2   248 53500   10  1.5 0.3 25  2   0.6 0.6                                     22.97   621.9   597.8   404.2   114.9   9.84    8.02    11  -261.06 0.000413412 5   2.5 10  1.3

Thanks heaps!

Relevant Background Info:
http://lmfit.github.io/lmfit-py/

https://pypi.python.org/pypi/lmfit/

https://github.com/lmfit/lmfit-py

http://cars9.uchicago.edu/software/python/lmfit/

Answer 1

LMFIT cannot optimize arrays of values. The only way to do a fit to multiple "data sets" is to write your objective function to do that -- it looks like you have that. But, your creation of the parameters with

# create a set of parameters
par = Parameters()
for i in range(int(len(x))):    
    par.add('Q',   value=Q_ini[i],  vary=True,  min= Q_min[i], max=Q_max[i])
    par.add('I',   value=I_ini[i],  vary=True,  min= I_min[i], max=I_max[i])
    ....

is just overwriting the definitions for Parameters Q and I. You probably want to index those per data set, with something like

# create a set of parameters
par = Parameters()
for i in range(int(len(x))):    
    par.add('Q_%i'%(i+1),   value=Q_ini[i],  vary=True,  min= Q_min[i], max=Q_max[i])
    par.add('I_%i'%(i+1),   value=I_ini[i],  vary=True,  min= I_min[i], max=I_max[i])
    ....

and then use those scalar 'per-data-set' parameters in your objective function. You can still have parameters global to all data sets, of course.