title: Fitting complex model using Python and lmfit?
body: <p>I would like to fit <a href="http://en.wikipedia.org/wiki/Ellipsometry" rel="nofollow noreferrer">ellipsometric</a> data to complex model using LMFit. Two measured parameters, <code>psi</code> and <code>delta</code>, are variables in a complex function <code>rho</code>.</p>
  
  <p>I could try with separating problem to real and imaginary part with <a href="https://stackoverflow.com/questions/20339234/python-and-lmfit-how-to-fit-multiple-datasets-with-shared-parameters">shared parameters</a> or <a href="https://stackoverflow.com/questions/11129812/fitting-piecewise-function-in-python">piecewise</a> 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: </p>
  
  <blockquote>
    <p>TypeError: no ordering relation is defined for complex numbers.</p>
  </blockquote>
  
  <p>Below is my code for real function fitting and my attempt at tackling complex fit problem:</p>
  
  <pre><code>    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)
      """
      #=================================================================
  </code></pre>
  
  <p><strong>Edit:</strong>
  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.</p>
  
  <pre><code>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 ...
  </code></pre>