Installing a complex model using Python and lmfit?

Question

I would like to put ellipsometric data in a complex model using lmfit. The two measured parameters, psi and delta , are variables in the complex function rho .

I could try to separate the problem into the real and imaginary parts using general parameters or picewise , but is there a way to do this directly with a complex function? Setting only the real part of the function works beautifully, but when I define a complex residual function, I get:

TypeError: An order relation is not defined for complex numbers.

Below is my code for installing a real function and my attempt to solve a complex 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 the problem with separated variables for the imaginary and real parts. The data must be in the form [[imaginary_data], [real_data]], the target function must return a 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 ...