¿Ajustando un modelo complejo usando Python y lmfit?

Me gustaría ajustar los datos elipsométricos a un modelo complejo usando lmfit. Dos parámetros medidos, psi y delta , son variables en una función compleja rho .

Podría intentar separar el problema en partes reales e imaginarias con parámetros compartidos o enfoque picewise , pero ¿hay alguna forma de hacerlo directamente con una función compleja? Ajustar solo una parte real de la función funciona a la perfección, pero cuando defino una función residual compleja obtengo:

TypeError: no se define una relación de orden para los números complejos.

A continuación se encuentra mi código para el ajuste de la función real y mi bash de abordar el problema del ajuste complejo:

  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) """ #================================================================= 

Edición: resolví problema con variables separadas para parte imaginaria y real. Los datos deben configurarse como [[imaginary_data], [real_data]], la función objective debe devolver la matriz 1D.

 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 ... 

Las Preguntas Frecuentes de lmfit sugieren simplemente tomar partes reales e imaginarias usando numpy.ndarray.view , lo que significa que no necesita pasar la separación de las partes real e imaginaria manualmente.

 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()