Advanced usage of Expressions in lmfit

[Jingzhong Yang]

Hi,

   I have a question about the advanced usage of expression in constraints because it seems like the example listed in the end does not detail too much for me.

   My question is that if I have several parameters need to be fit such as a, b, c, d……(they have their own constraints) and I also have a parameter (eg: y) which has an extremely complicated relationship with previous ones and a simple constraint for it as well, so that it is unable for me to define it in the expression shortly. In this case, I think I need to define a specific function to express the relationship and cite this function in 'y' parameters expression. The thing is that when I try to operate the example code, it gave me an error called '' __init__() takes at least 3 arguments (1 given)'' 

  So, could you please explain me more about the usage of this advanced expression for constraints

  Bests

  Jingzhong 

[Matt Newville]

Hi Jingzhong,

Please ask an actual question that includes a "?" symbol and that can be answered. 

It is always helpful to come up with and post a minimal example that shows any problem you are asking about.  It is helpful to us because then we will know what you are trying to do, and it may be helpful to you because it will help you distill the problem.

 
For sure, if you're going to create a Minimizer(), you'll have to pass in an objective function and parameters object.  Is that your question?  I do not know.

You should be able to do something like this:

    from lmfit import Parameters

    def mylorentzian(x, amp, cen, wid):
        return (amp / (1 + ((x-cen) / wid)**2))

    params = Parameters()
    params.add('a1', value=10)
    params.add('c1', value=2)
    params.add('s1', value=0.5)

    params._asteval.symtable['fcn'] = mylorentzian
    params.add('x12', expr='fcn(12.0, a1, c1, s1)')

    for name, par in params.items():
        print(name, par)

which will print
   a1 <Parameter 'a1', 10, bounds=[-inf:inf]>
   c1 <Parameter 'c1', 2, bounds=[-inf:inf]>
   s1 <Parameter 's1', 0.5, bounds=[-inf:inf]>
   x12 <Parameter 'x', 0.02493765586034913, bounds=[-inf:inf], expr='fcn(12.0, a1, c1, s1)'>

Hope that helps.

--Matt

[Jingzhong Yang]

Hello Matt,

  Thank you very much for your reply and your code above definitely helps me a lot. 

  Now I still have a similar question about the advanced usage of the fitting constraint like the example below, in which I do not add the data. 

  My question is that if I also want the mylorentzian() constructed by variables 'a','b','c','d' always satisfies some boundary condition like [1,100) during the fitting. How can I realize this constraint in the code?

from lmfit import Parameters

def residuals(params,s,data):

    a=params['a']

    b=params['b']

    c=params['c']

    d=params['d']

    model=a*s[:,0]+b*s[:,1]+c*s[:,2]+d*s[:,3]

    return model-data

def mylorentzian(a,b,c,d):

    return (a/ (1 + ((b-c) / d)**2))

s=s_original

data=s_polution

params = Parameters()

params.add('a', value=10,vary=True,min=0,max=15)

params.add('b', value=2,vary=True,min=0,max=3)

params.add('c', value=0.5,vary=True,min=0,max=1)

params.add('d',value=2,vary=True,min=0,max=4)

out=Minimizer(residuals,params,args=(s,data))

 

Best regards

Jingzhong

[Matt Newville]

Hi Jingzhong,

On Thu, Apr 5, 2018 at 10:44 AM, Jingzhong Yang <yangjingzh...@gmail.com> wrote:

Hello Matt,

  Thank you very much for your reply and your code above definitely helps me a lot. 

  Now I still have a similar question about the advanced usage of the fitting constraint like the example below, in which I do not add the data. 

  My question is that if I also want the mylorentzian() constructed by variables 'a','b','c','d' always satisfies some boundary condition like [1,100) during the fitting. How can I realize this constraint in the code?

Well, what did you try?   There are a couple of ways that could work:

  a) modify your function to impose the bounds, perhaps as

    def mylorentzian(a, b, c, d, minval=0, maxval=10):

     

        val = (a/ (1 + ((b-c) / d)**2))

        val = min(maxval, max(minval, val))

        return val

  b) explicitly set bounds on the parameter:

      params._asteval.symtable['fcn'] = mylorentzian
      params.add('x12', expr='fcn(12.0, a1, c1, s1)')

      params['x12'].min = 10

      params['x12'].max = 90

--Matt

[Jingzhong Yang]

Dear Matt,

  According to your suggestion, I  program the code to fit the variables a, b, c, d and these four variables are necessary for the residual calculation function and they are marked with red below. The second function (mylorentzian()) needs to be imposed as the boundary condition for a,b,c,d but the variable 'x12' is not used in the residual minimization. Is it possible for this code to fit out the a,b,c,d which satisfy the second type boundary condition  10<=(a/ (1 + ((b-c) / d)**2))<=90 ?

from lmfit import Parameters

def residuals(params,s,data):

    a=params['a']

    b=params['b']

    c=params['c']

    d=params['d']

    model=a*s[:,0]+b*s[:,1]+c*s[:,2]+d*s[:,3]

    return model-data

def mylorentzian(a, b, c, d):

      val = (a/ (1 + ((b-c) / d)**2))

      return val

s=s_original

data=s_polution

params = Parameters()

params.add('a', value=10,vary=True,min=0,max=15)

params.add('b', value=2,vary=True,min=0,max=3)

params.add('c', value=0.5,vary=True,min=0,max=1)

params.add('d',value=2,vary=True,min=0,max=4)

params._asteval.symtable['fcn'] = mylorentzian

params.add('x12', value=20,vary=True, min=10, max=90, expr='fcn(a,b,c,d)')

out=minimize(residuals,params,args=(s,data))

Bests

Jingzhong

[Matt Newville]

Hi Jingzhong,

On Sun, Apr 8, 2018 at 10:39 AM, Jingzhong Yang <yangjingzh...@gmail.com> wrote:

Dear Matt,

  According to your suggestion, I  program the code to fit the variables a, b, c, d and these four variables are necessary for the residual calculation function and they are marked with red below. The second function (mylorentzian()) needs to be imposed as the boundary condition for a,b,c,d but the variable 'x12' is not used in the residual minimization. Is it possible for this code to fit out the a,b,c,d which satisfy the second type boundary condition  10<=(a/ (1 + ((b-c) / d)**2))<=90 ?

Well, did you try this?  I'm pretty sure I asked you this once before -- you didn't answer then either. 

As I said before, use

      params._asteval.symtable['fcn'] = mylorentzian
      params.add('x12', expr='fcn(12.0, a1, c1, s1)')

      params['x12'].min = 10

      params['x12'].max = 90

Write a test script changing parameter values to test if that works the way you want it to.  If it doesn't, post the full test script and the output you get and what you expected.

--Matt

[Jingzhong Yang]

Hi Matt

  Please check this code for instance and what I would expect from it is that the variables should obey  1.5<m+n<2. However, it still does not work under this constraint.

import numpy as np

from lmfit import Parameters, minimize,fit_report

def residuals(params,s,data):

    m=params['m']

    n=params['n']

    residual=m*s+n-data

    return residual

def myfunc(m,n):

    y=m+n

    return y

s=np.arange(20,dtype=float)

data=np.linspace(0,19,20)#+np.random.rand(20)

params=Parameters()

params.add('m',value=0,vary=True)

params.add('n',value=0,vary=True)

params._asteval.symtable['fcn'] = myfunc

params.add('y1',vary=True, expr='fcn(m,n)') 

params['y1'].min=1.5

params['y1'].max=2

out=minimize(residuals,params,args=(s,data))

report=fit_report(out)

print report

[[Fit Statistics]]

    # fitting method   = leastsq

    # function evals   = 9

    # data points      = 20

    # variables        = 2

    chi-square         = 1.2252e-28

    reduced chi-square = 6.8067e-30

    Akaike info crit   = -1341.30022

    Bayesian info crit = -1339.30876

[[Variables]]

    m:   1.00000000 +/- 5.2495e-17 (0.00%) (init = 0)

    n:   9.1065e-24 +/- 2.6090e-15 (28649540321.09%) (init = 0)

    y1:  1.50000000 +/- 5.2495e-17 (0.00%) == 'fcn(m,n)'

[Matt Newville]

Jingzhong,

On Mon, Apr 9, 2018 at 7:38 AM, Jingzhong Yang <yangjingzh...@gmail.com> wrote:

Hi Matt

  Please check this code for instance and what I would expect from it is that the variables should obey  1.5<m+n<2. However, it still does not work under this constraint.

You did not actually set a constraint that 1.5 < m+n < 2.   You set a constraint that 1.5 < y1 < 2.   And that did indeed work.

If you want m + n to be between 1.5 and 2 you could do this (see https://lmfit.github.io/lmfit-py/constraints.html#using-inequality-constraints):
 

   params=Parameters()

   params.add('m',value=0,vary=True)

   params.add('delta', value=1.75,  min=1.50, max=2.00)

   params.add('n',  expr='delta - m')

Now `m` varies, and `n` may take any value that satisfies  1.5 < m + n < 2.

 

--
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+unsubscribe@googlegroups.com.
To post to this group, send email to lmfi...@googlegroups.com.
Visit this group at https://groups.google.com/group/lmfit-py.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/9a5be506-dd9f-40bd-9f01-6180c877b3d6%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Jingzhong Yang]

Hi Matt,

  Thank you for your reply.

  The problem is that in my actual code, there are up to 16 variables need to be fitted and they are mixed together with extremely sophisticated relationship instead of easy format like 1.5<m+n<2. So, is there another way to impose the complicate boundary condition in the fitting? Thank you in advanced.

  Best regards

  Jingzhong

To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.

To post to this group, send email to lmfi...@googlegroups.com.
Visit this group at https://groups.google.com/group/lmfit-py.

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbpNDkMQqVGXXTzbf6644%2BmdK8ej9NuET4G5vJmw_E1Cew%40mail.gmail.com.

[Matt Newville]

Jingzhong,

On Mon, Apr 9, 2018 at 8:24 AM, Jingzhong Yang <yangjingzh...@gmail.com> wrote:

Hi Matt,

  Thank you for your reply.

  The problem is that in my actual code, there are up to 16 variables need to be fitted and they are mixed together with extremely sophisticated relationship instead of easy format like 1.5<m+n<2. So, is there another way to impose the complicate boundary condition in the fitting? Thank you in advanced.

  Best regards

  Jingzhong

Hm, well there are other ways.  You could write your residual function to impose them, for example.  Or you could write out complicated expressions. It's hard to give more specific advice than that without knowing more details. 

--Matt

[Jingzhong Yang]

Hi Matt,

  Please check the codes and files attached, from which you can get a code for the fitting and a code to check the fitting result. My aim is to guarantee that the boundary conditions applied can work on the 16 variables (from a to p) ready to be fitted in residuals function. 

  The constraints for 16 variables:

Md_det_abs_f > 0

0<P_f<1

0<D_f<1 (equivalent to ‘a**2>sqrt(b**2+c**2+d**2)')

m_r_delta_det_abs_f > 0

p_delta_f > 0

0<dep_power_f<1

0< R_f <pi

[Matt Newville]

Jingzhong,

On Wed, Apr 11, 2018 at 4:08 AM, Jingzhong Yang <yangjingzhong.optics@gmail.com> wrote:

Hi Matt,

  Please check the codes and files attached, from which you can get a code for the fitting and a code to check the fitting result. My aim is to guarantee that the boundary conditions applied can work on the 16 variables (from a to p) ready to be fitted in residuals function. 

You keep asking to "check the code", and you keep giving no indication that you have actually run the code yourself.  I keep asking if you tried this, and you keep not responding to my direct questions.   Have you run this code?  What results and/or error messages do you get?   These are not idle questions.

Your code you posted is far too complex and hard to read to give any prediction about whether they do what you intend.  Normally, I would trust that someone has tested the various functions and convinced themselves that the results are correct.  If you have not done that, you should.

Your parameter definitions look OK to me, at least in the sense that they ought to do something...  whether that something is actually what you intend is harder to say, but it is testable.   Again, if you have not tested them by trying different values for the variable parameters and seeing what values the constrained parameters take, I strongly recommend you do that.  I don't think anyone else is going to do that for you.

I do notice that your objective function calculates `residual` as  sqrt (model1 - data1)**2 + (model2 - data2)**2 + ...).   I doubt that this is what you want to do.  Normally, one would return

        np.concatenate( (model1 - data1), (model2 - data2), ...) )

That is, not taking the square root of the sum of squares yourself.

Good luck,

--Matt