Merge pull request #901 from lmfit/complex_eval_uncertain

Complex eval uncertainty

lmfit/model.py

@@ -927,7 +927,7 @@ def eval(self, params=None, **kwargs):
         values may be Python `float`, `int`, or  `complex` values.
 
         """
-        return self.func(**self.make_funcargs(params, kwargs))
+        return coerce_arraylike(self.func(**self.make_funcargs(params, kwargs)))
 
     @property
     def components(self):
@@ -1585,18 +1585,20 @@ def eval_uncertainty(self, params=None, sigma=1, **kwargs):
 
         nvarys = self.nvarys
         # ensure fjac and df2 are correct size if independent var updated by kwargs
-        ndata = self.model.eval(params, **userkws).size
+        feval = self.model.eval(params, **userkws)
+        ndata = len(feval.view('float64'))        # allows feval to be complex
         covar = self.covar
         if any(p.stderr is None for p in params.values()):
             return np.zeros(ndata)
 
-        fjac = {'0': np.zeros((nvarys, ndata))}  # '0' signify 'Full', an invalid prefix
-        df2 = {'0': np.zeros(ndata)}
+        # '0' would be an invalid prefix, here signifying 'Full'
+        fjac = {'0': np.zeros((nvarys, ndata), dtype='float64')}
+        df2 = {'0': np.zeros(ndata, dtype='float64')}
 
         for comp in self.components:
             label = comp.prefix if len(comp.prefix) > 1 else comp._name
-            fjac[label] = np.zeros((nvarys, ndata))
-            df2[label] = np.zeros(ndata)
+            fjac[label] = np.zeros((nvarys, ndata), dtype='float64')
+            df2[label] = np.zeros(ndata, dtype='float64')
 
         # find derivative by hand!
         pars = params.copy()
@@ -1614,7 +1616,8 @@ def eval_uncertainty(self, params=None, sigma=1, **kwargs):
 
             pars[pname].value = val0
             for key in fjac:
-                fjac[key][i] = (res1[key] - res2[key]) / (2*dval)
+                fjac[key][i] = (res1[key].view('float64')
+                                - res2[key].view('float64')) / (2*dval)
 
         for i in range(nvarys):
             for j in range(nvarys):
@@ -1627,6 +1630,12 @@ def eval_uncertainty(self, params=None, sigma=1, **kwargs):
             prob = erf(sigma/np.sqrt(2))
 
         scale = t.ppf((prob+1)/2.0, self.ndata-nvarys)
+
+        # for complex data, convert back to real/imag pairs
+        if feval.dtype in ('complex64', 'complex128'):
+            for key in fjac:
+                df2[key] = df2[key][0::2] + 1j * df2[key][1::2]
+
         self.dely = scale * np.sqrt(df2.pop('0'))
 
         self.dely_comps = {}

tests/test_model.py

@@ -1496,3 +1496,42 @@ def test_make_params_valuetypes():
 
     with pytest.raises(TypeError):
         pars = mod.make_params(amplitude={}, frequency=2, shift=7)
+
+
+def test_complex_model_eval_uncertainty():
+    """Github #900"""
+    def cmplx(f, omega, areal, aimag, off, sigma):
+        return (areal*np.cos(f*omega + off) + 1j*aimag*np.sin(f*omega + off))*np.exp(-f/sigma)
+
+    f = np.linspace(0, 10, 501)
+    dat = cmplx(f, 4, 10, 5, 0.2, 4.5) + (0.1 + 0.2j)*np.random.normal(scale=0.25, size=len(f))
+    mod = Model(cmplx)
+    params = mod.make_params(omega=5, areal=5, aimag=5,
+                             off={'value': 0.5, 'min': -2, 'max': 2},
+                             sigma={'value': 3, 'min': 1.e-5, 'max': 1000})
+
+    result = mod.fit(dat, params=params, f=f)
+    dfit = result.eval_uncertainty()
+    assert len(dfit) == len(f)
+    assert dfit.dtype == 'complex128'
+
+
+def test_compositemodel_returning_list():
+    """Github #875"""
+    def lin1(x, k):
+        return [k*x1 for x1 in x]
+
+    def lin2(x, k):
+        return [k*x1 for x1 in x]
+
+    y = np.linspace(0, 100, 100)
+    x = np.linspace(0, 100, 100)
+
+    Model1 = Model(lin1, independent_vars=["x"], prefix="m1_")
+    Model2 = Model(lin2, independent_vars=["x"], prefix="m2_")
+    ModelSum = Model1 + Model2
+    pars = Parameters()
+    pars.add('m1_k', value=0.5)
+    pars.add('m2_k', value=0.5)
+    result = ModelSum.fit(y, pars, x=x)
+    assert len(result.best_fit) == len(x)