Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fitmodel: fix uncertainty quantification when fitting dipolar evolution functions #172

Merged
merged 4 commits into from
May 4, 2021
Merged

fitmodel: fix uncertainty quantification when fitting dipolar evolution functions #172

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0c6323d
fitmodel: fix errors in uncertainty quantification for dipolar evolut…
luisfabib May 3, 2021
0e65ef3
modify tests that missed the error
luisfabib May 3, 2021
4553113
fitmodel: fixed error found by tests
luisfabib May 3, 2021
1bb1e8c
Merge branch 'main' into fix/fitmodel_dipevofcn
luisfabib May 4, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 18 additions & 34 deletions deerlab/fitmodel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -358,7 +358,7 @@ def multiPathwayModel(par):
return Ks, Bs
# =========================================================================

def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1):
def splituq(param_uq,Pfit_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1):
# =========================================================================
"""
Uncertainty quantification
Expand All @@ -370,26 +370,11 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1):
# Pre-allocation
paruq_bg,paruq_ex,Bfit_uq,Vmod_uq,Vunmod_uq,Vfit_uq = ([],[],[],[],[],[])

# Retrieve full covariance matrix
if isinstance(full_uq,list):
paruq = full_uq[0]
else:
paruq = full_uq
covmat = paruq.covmat

Nparam = len(parfit_)
paramidx = np.arange(Nparam)

# Full parameter set uncertainty
# -------------------------------
subcovmat = covmat[np.ix_(paramidx,paramidx)]
paruq = UQResult('covariance',parfit_,subcovmat,lb,ub)

# Background parameters uncertainty
# ---------------------------------
for jj in range(nSignals):
if includeBackground[jj]:
bgsubcovmat = paruq.covmat[np.ix_(bgidx[jj],bgidx[jj])]
bgsubcovmat = param_uq.covmat[np.ix_(bgidx[jj],bgidx[jj])]
paruq_bg.append( UQResult('covariance',parfit_[bgidx[jj]],bgsubcovmat,lb[bgidx[jj]],ub[bgidx[jj]]))
else:
paruq_bg.append([None])
Expand All @@ -398,15 +383,15 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1):
# ----------------------------------
for jj in range(nSignals):
if includeExperiment[jj]:
exsubcovmat = paruq.covmat[np.ix_(exidx[jj],exidx[jj])]
exsubcovmat = param_uq.covmat[np.ix_(exidx[jj],exidx[jj])]
paruq_ex.append( UQResult('covariance',parfit_[exidx[jj]],exsubcovmat,lb[exidx[jj]],ub[exidx[jj]]))
else:
paruq_ex.append([None])

# Distribution parameters uncertainty
# ------------------------------------
if parametricDistribution:
ddsubcovmat = paruq.covmat[np.ix_(ddidx,ddidx)]
ddsubcovmat = param_uq.covmat[np.ix_(ddidx,ddidx)]
paruq_dd = UQResult('covariance',parfit_[ddidx],ddsubcovmat,lb[ddidx],ub[ddidx])
else:
paruq_dd = [None]
Expand All @@ -420,15 +405,13 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1):
Pfcn = lambda par: dd_model(r,par[ddidx])
else:
Pfcn = lambda _: np.ones_like(r)/np.trapz(np.ones_like(r),r)
Pfit_uq = paruq.propagate(Pfcn,nonneg)
else:
Pfit_uq = full_uq[1]
Pfit_uq = param_uq.propagate(Pfcn,nonneg)

# Background uncertainty
# -----------------------
for jj in range(nSignals):
if includeExperiment[jj]:
Bfit_uq.append( paruq.propagate(lambda par:scales[jj]*multiPathwayModel(par)[1][jj]) )
Bfit_uq.append( param_uq.propagate(lambda par:scales[jj]*multiPathwayModel(par)[1][jj]) )
else:
Bfit_uq.append([None])

Expand All @@ -439,7 +422,7 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1):
Lam0fcn = lambda par: ex_model[jj](par)[0]
Bfcn = lambda par: scales[jj]*multiPathwayModel(par)[1][jj]
Vunmod_fcn = lambda par: Lam0fcn(par[exidx[jj]])*Bfcn(par)
Vunmod_uq.append( paruq.propagate(lambda par:Vunmod_fcn(par)) )
Vunmod_uq.append( param_uq.propagate(lambda par:Vunmod_fcn(par)) )
else:
Vunmod_uq.append([None])

Expand All @@ -449,26 +432,28 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1):
if includeForeground and parametricDistribution:
# Full parametric signal
Vmodel = lambda par: scales[jj]*multiPathwayModel(par)[0][jj]@Pfcn(par[ddidx])
Vfit_uq.append( paruq.propagate(Vmodel))
Vfit_uq.append( param_uq.propagate(Vmodel))
elif includeForeground and np.all(~includeExperiment & ~includeBackground):
Vmodel = lambda _: Kfit[jj]@Pfit
# Dipola evolution function
J = Kfit[jj]
Vcovmat = J@covmat@J.T
Vcovmat = J@Pfit_uq.covmat@J.T
Vfit_uq.append( UQResult('covariance',Vfit[jj],Vcovmat))
elif includeForeground:
# Parametric signal with parameter-free distribution
Vmodel = lambda par: scales[jj]*multiPathwayModel(par[paramidx])[0][jj]@Pfit
Vfit_uq.append( paruq.propagate(Vmodel) )
Vmodel = lambda par: scales[jj]*multiPathwayModel(par)[0][jj]@Pfit
Vfit_uq.append( param_uq.propagate(Vmodel) )
else:
Vfit_uq.append([None])

# Modulated contribution uncertainty
# -----------------------------
for jj in range(nSignals):
if includeForeground:
if includeForeground and np.all(~includeExperiment & ~includeBackground):
Vmod_uq.append(Vfit_uq)
elif includeForeground:
Vmod_fcn = lambda par: Vmodel(par) - Vunmod_fcn(par)
Vmod_uq.append( paruq.propagate(Vmod_fcn) )
Vmod_uq.append( param_uq.propagate(Vmod_fcn))
else:
Vmod_uq.append([None])

Expand Down Expand Up @@ -525,7 +510,7 @@ def regularization_analysis(Vexp):
parfit = np.asarray([None])

if uqanalysis and uq=='covariance':
Vfit_uq, Pfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(Pfit_uq,Pfit,Vfit,Bfit,parfit,Ks,scales)
Vfit_uq, Pfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(None,Pfit_uq,Pfit,Vfit,Bfit,parfit,Ks,scales)
return fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit, Pfit_uq, Vfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd, scales, alphaopt
else:
return fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit, scales, alphaopt
Expand Down Expand Up @@ -565,7 +550,7 @@ def nonlinear_lsq_analysis(Vexp):
Vmod, Vunmod = calculate_Vmod_Vunmod(parfit,Vfit,Bfit,scales)

if uqanalysis and uq=='covariance':
Vfit_uq, Pfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(param_uq,Pfit,Vfit,Bfit,parfit,None, scales)
Vfit_uq, Pfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(param_uq,None,Pfit,Vfit,Bfit,parfit,None, scales)
return fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit, Pfit_uq, Vfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd,scales,alphaopt
else:
return fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit, scales, alphaopt
Expand Down Expand Up @@ -607,7 +592,6 @@ def scale_constraint(nonlinpar):
Pfit = fit.lin
param_uq = fit.nonlinUncert
Pfit_uq = fit.linUncert
snlls_uq = [param_uq,Pfit_uq]
alphaopt = fit.regparam
scales = fit.scale

Expand All @@ -628,7 +612,7 @@ def scale_constraint(nonlinpar):
Vmod, Vunmod = calculate_Vmod_Vunmod(parfit,Vfit,Bfit,scales)

if uqanalysis and uq=='covariance':
Vfit_uq, _, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(snlls_uq, Pfit, Vfit, Bfit, parfit, Kfit, scales)
Vfit_uq, _, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(param_uq,Pfit_uq, Pfit, Vfit, Bfit, parfit, Kfit, scales)
return fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit, Pfit_uq, Vfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd,scales,alphaopt
else:
return fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit, scales, alphaopt
Expand Down
4 changes: 2 additions & 2 deletions test/test_fitmodel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ def test_dipevo_function():
P = dd_gauss(r,[4.5, 0.25])
K = dipolarkernel(t,r)
V = K@P
fit = fitmodel(V,t,r,'P',None,None,uq=None)
fit = fitmodel(V,t,r,'P',None,None)
assert ovl(P,fit.P) > 0.90
# ======================================================================

Expand All @@ -100,7 +100,7 @@ def test_form_factor():
P = dd_gauss(r,[4.5, 0.25])
K = dipolarkernel(t,r,mod=0.3)
V = K@P
fit = fitmodel(V,t,r,'P',None,ex_4pdeer,uq=None)
fit = fitmodel(V,t,r,'P',None,ex_4pdeer)
assert ovl(P,fit.P) > 0.90
# ======================================================================

Expand Down