Fitting complex model using Python and lmfit?

Question

I would like to fit ellipsometric data to complex model using LMFit. Two measured parameters, psi and delta, are variables in a complex function rho.

I could try with separating problem to real and imaginary part with shared parameters or piecewise approach, but is there any way to do it directly with complex function? Fitting only real part of function works beautifully, but when I define complex residual function I get:

TypeError: no ordering relation is defined for complex numbers.

Below is my code for real function fitting and my attempt at tackling complex fit problem:

    from __future__ import division
    from __future__ import print_function
    import numpy as np
    from pylab import *
    from lmfit import minimize, Parameters, Parameter, report_errors


    #=================================================================
    #             MODEL

    def r01_p(eps2, th):
        c=cos(th)
        s=(sin(th))**2

        stev= sqrt(eps2) * c - sqrt(1-(s / eps2))
        imen= sqrt(eps2) * c + sqrt(1-(s / eps2))
        return stev/imen

    def r01_s(eps2, th):
        c=cos(th)
        s=(sin(th))**2

        stev= c - sqrt(eps2) * sqrt(1-(s/eps2))
        imen= c + sqrt(eps2) * sqrt(1-(s/eps2))
        return stev/imen


    def rho(eps2, th):
        return r01_p(eps2, th)/r01_s(eps2, th)

    def psi(eps2, th):
        x1=abs(r01_p(eps2, th))
        x2=abs(r01_s(eps2, th))
        return np.arctan2(x1,x2)

    #=================================================================
    #                   REAL FIT
    #

    #%%

    # generate data from model  
    th=linspace(deg2rad(45),deg2rad(70),70-45)
    error=0.01
    var_re=np.random.normal(size=len(th), scale=error)
    data = psi(2,th) + var_re

    # residual function
    def residuals(params, th, data):
        eps2 = params['eps2'].value

        diff = psi(eps2, th) - data
        return diff

    # create a set of Parameters
    params = Parameters()
    params.add('eps2',   value= 1.0,  min=1.5, max=3.0)


    # do fit, here with leastsq model
    result = minimize(residuals, params, args=(th, data),method="leastsq")

    # calculate final result
    final = data + result.residual

    # write error report
    report_errors(params)


    # try to plot results
    th, data, final=rad2deg([th, data, final])
    try:
        import pylab
        clf()
        fig=plot(th, data, 'r o',
                 th, final, 'b')
        setp(fig,lw=2.)
        xlabel(r'$\theta$ $(^{\circ})$', size=20)
        ylabel(r'$\psi$ $(^{\circ})$',size=20)

        show()

    except:
        pass

    #%%
    #=================================================================
    #                   COMPLEX FIT

    # TypeError: no ordering relation is defined for complex numbers

    """
    # data from model with added noise   
    th=linspace(deg2rad(45),deg2rad(70),70-45)
    error=0.001
    var_re=np.random.normal(size=len(th), scale=error)
    var_im=np.random.normal(size=len(th), scale=error) * 1j

    data = rho(4-1j,th) + var_re + var_im


    # residual function
    def residuals(params, th, data):
        eps2 = params['eps2'].value

        diff = rho(eps2, th) - data
        return np.abs(diff)

    # create a set of Parameters
    params = Parameters()
    params.add('eps2',   value= 1.5+1j,  min=1+1j, max=3+3j)


    # do fit, here with leastsq model
    result = minimize(residuals, params, args=(th, data),method="leastsq")

    # calculate final result
    final = data + result.residual

    # write error report
    report_errors(params)
    """
    #=================================================================

Edit: I solved problem with separated variables for imaginary and real part. Data should be shaped as [[imaginary_data],[real_data]], objective function must return 1D array.

def objective(params, th_data, data):
    eps_re  = params['eps_re'].value
    eps_im  = params['eps_im'].value
    d       = params['d'].value

    residual_delta = data[0,:] - delta(eps_re - eps_im*1j, d, frac, lambd, th_data)
    residual_psi   = data[1,:] - psi(eps_re - eps_im*1j, d, frac, lambd, th_data)

    return np.append(residual_delta,residual_psi)

# create a set of Parameters
params = Parameters()
params.add('eps_re',   value= 1.5,  min=1.0,      max=5  )
params.add('eps_im',   value= 1.0,  min=0.0,      max=5  )
params.add('d',        value= 10.0,   min=5.0,    max=100.0   )


# All available methods
methods=['leastsq','nelder','lbfgsb','anneal','powell','cobyla','slsqp']
# Chosen method
#metoda='leastsq'

# run the global fit to all the data sets
result = minimize(objective, params, args=(th_data,data),method=metoda))

....

return ...

Answer 1

The lmfit FAQ suggests simply taking both real and imaginary parts by using numpy.ndarray.view, which means you don't need to go through the separation of the real and imaginary parts manually.

def residuals(params, th, data):
    eps2 = params['eps2'].value

    diff = rho(eps2, th) - data
    # The only change required is to use view instead of abs.
    return diff.view()