Skip to content

Commit

Permalink
fix ModelResult.eval_uncertainty() for complex model
Browse files Browse the repository at this point in the history
  • Loading branch information
newville committed Jul 5, 2023
1 parent d0a5936 commit 5f7af07
Showing 1 changed file with 15 additions and 6 deletions.
21 changes: 15 additions & 6 deletions lmfit/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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):
Expand All @@ -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 = {}
Expand Down

0 comments on commit 5f7af07

Please sign in to comment.