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.

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

fitsignal: Merge uq and uqanalysis keyword arguments #98

Merged
merged 2 commits into from
Mar 18, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
18 changes: 10 additions & 8 deletions deerlab/fitsignal.py
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,9 @@ def fitsignal(Vexp, t, r, dd_model='P', bg_model=bg_hom3d, ex_model=ex_4pdeer,
* ``'covariance'`` - Covariance-based uncertainty quantification. Fast, but approximate.
* ``'bootstrap'`` - Bootstrapped uncertainty quantification. Slow, but accurate. By default, 1000 bootstrap
samples are used. Alternatively, a different number can be specified as follows ``uq=['bootstrap',Nsamples]``.
* ``None`` - Disable the uncertainty quantification analysis.

The default is ``'covariance'``.

weights : array_like, optional
Array of weighting coefficients for the individual signals in global fitting,
Expand Down Expand Up @@ -124,10 +127,6 @@ def fitsignal(Vexp, t, r, dd_model='P', bg_model=bg_hom3d, ex_model=ex_4pdeer,

verbose : boolean, optional
Enable/disable printing a table of fit results, by default is disabled

uqanalysis : boolean, optional
Enable/disable the uncertainty quantification analysis, by default it is enabled.



Returns
Expand Down Expand Up @@ -249,15 +248,18 @@ def fitsignal(Vexp, t, r, dd_model='P', bg_model=bg_hom3d, ex_model=ex_4pdeer,

# Default bootstrap samples
bootsamples = 1000
if isinstance(uq, str):
if isinstance(uq, str) or uq==None:
uq = [uq]
if uq[0]!='bootstrap' and uq[0]!='covariance':
raise KeyError("Uncertainty quantification must be either 'covariance' or 'bootstrap'.")

if uq[0]!='bootstrap' and uq[0]!='covariance'and uq[0]!=None:
raise KeyError("Uncertainty quantification must be either 'covariance', 'bootstrap', or None.")
if uq[0]=='bootstrap':
# OVerride default if user has specified bootstraped samples
if len(uq)>1: bootsamples = uq[1]
uq = uq[0]
if uq is None:
uqanalysis = False
else:
uqanalysis = True

# Combine input boundary and start conditions
par0 = [[] if par0_i is None else par0_i for par0_i in [dd_par0,bg_par0,ex_par0]]
Expand Down
46 changes: 23 additions & 23 deletions test/test_fitsignal.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ def assert_experiment_model(model):
K = dipolarkernel(t,r,pathways,Bmodel)
V = K@P

fit = fitsignal(V,t,r,'P',bg_exp,model,uqanalysis=False)
fit = fitsignal(V,t,r,'P',bg_exp,model,uq=None)

assert ovl(P,fit.P) > 0.90
# --------------------------------------------------------------------
Expand Down Expand Up @@ -88,7 +88,7 @@ def test_dipevo_function():
P = dd_gauss(r,[4.5, 0.25])
K = dipolarkernel(t,r)
V = K@P
fit = fitsignal(V,t,r,'P',None,None,uqanalysis=False)
fit = fitsignal(V,t,r,'P',None,None,uq=None)
assert ovl(P,fit.P) > 0.90
# ======================================================================

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

Expand All @@ -117,7 +117,7 @@ def test_full_parametric():
B = bg_exp(t,lam*kappa)
V = dipolarkernel(t,r,lam,B)@P

fit = fitsignal(V,t,r,dd_gauss,bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal(V,t,r,dd_gauss,bg_exp,ex_4pdeer,uq=None)

assert ovl(P,fit.P) > 0.99
# ======================================================================
Expand All @@ -131,7 +131,7 @@ def test_no_foreground():
k = 0.2
B = bg_exp(t,k)

fit = fitsignal(B,t,r,None,bg_exp,None,uqanalysis=False)
fit = fitsignal(B,t,r,None,bg_exp,None,uq=None)

assert abs(k-fit.bgparam) < 1
# ======================================================================
Expand All @@ -149,7 +149,7 @@ def test_start_values():
K = dipolarkernel(t,r,lam,Bmodel)
V = K@P

fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,bg_par0=0.5,ex_par0=0.5,uqanalysis=False)
fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,bg_par0=0.5,ex_par0=0.5,uq=None)

assert ovl(P,fit.P) > 0.95
# ======================================================================
Expand All @@ -167,7 +167,7 @@ def test_boundaries():
K = dipolarkernel(t,r,lam,Bmodel)
V = K@P

fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer, uqanalysis=False,
fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer, uq=None,
bg_par0=0.4, bg_lb=0.2, bg_ub=0.5,
ex_par0=0.4, ex_lb=0.2, ex_ub=0.5)

Expand All @@ -187,7 +187,7 @@ def test_boundaries_adjust_bg():
K = dipolarkernel(t,r,lam,Bmodel)
V = K@P

fit = fitsignal(V,t,r,'P',bg_hom3d,ex_4pdeer, uqanalysis=False, bg_lb=70, bg_ub=90)
fit = fitsignal(V,t,r,'P',bg_hom3d,ex_4pdeer, uq=None, bg_lb=70, bg_ub=90)

assert ovl(P,fit.P) > 0.95
# ======================================================================
Expand All @@ -205,7 +205,7 @@ def test_boundaries_adjust_ex():
K = dipolarkernel(t,r,lam,Bmodel)
V = K@P

fit = fitsignal(V,t,r,'P',bg_hom3d,ex_4pdeer, uqanalysis=False, ex_lb=0.55, ex_ub=0.65)
fit = fitsignal(V,t,r,'P',bg_hom3d,ex_4pdeer, uq=None, ex_lb=0.55, ex_ub=0.65)

assert ovl(P,fit.P) > 0.95
# ======================================================================
Expand All @@ -223,7 +223,7 @@ def test_boundaries_adjust_dd():
K = dipolarkernel(t,r,lam,Bmodel)
V = K@P

fit = fitsignal(V,t,r,dd_gauss,bg_hom3d,ex_4pdeer, uqanalysis=False, dd_lb=[4,0.3],dd_ub=[6,0.5])
fit = fitsignal(V,t,r,dd_gauss,bg_hom3d,ex_4pdeer, uq=None, dd_lb=[4,0.3],dd_ub=[6,0.5])

assert ovl(P,fit.P) > 0.95
# ======================================================================
Expand All @@ -248,7 +248,7 @@ def test_global_4pdeer():
t2 = np.linspace(0,5,250)
V2 = dipolarkernel(t2,r,pathways,Bmodel)@P

fit = fitsignal([V1,V2],[t1,t2],r,'P',bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal([V1,V2],[t1,t2],r,'P',bg_exp,ex_4pdeer,uq=None)

assert ovl(P,fit.P) > 0.90
# ======================================================================
Expand All @@ -267,7 +267,7 @@ def test_global_full_parametric():
V1 = dipolarkernel(t1,r,lam,B)@P
V2 = dipolarkernel(t2,r,lam,B)@P

fit = fitsignal([V1,V2],[t1,t2],r,dd_gauss,bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal([V1,V2],[t1,t2],r,dd_gauss,bg_exp,ex_4pdeer,uq=None)

assert ovl(P,fit.P) > 0.99
# ======================================================================
Expand All @@ -286,7 +286,7 @@ def test_global_mixed_backgrounds():
V1 = dipolarkernel(t1,r,lam,B)@P
V2 = dipolarkernel(t2,r,lam)@P

fit = fitsignal([V1,V2],[t1,t2],r,dd_gauss,[bg_exp,None],ex_4pdeer,uqanalysis=False)
fit = fitsignal([V1,V2],[t1,t2],r,dd_gauss,[bg_exp,None],ex_4pdeer,uq=None)

assert ovl(P,fit.P) > 0.90
# ======================================================================
Expand All @@ -303,7 +303,7 @@ def test_global_mixed_experiments():
V1 = dipolarkernel(t1,r,lam)@P
V2 = dipolarkernel(t2,r)@P

fit = fitsignal([V1,V2],[t1,t2],r,dd_gauss,None,[ex_4pdeer,None],uqanalysis=False)
fit = fitsignal([V1,V2],[t1,t2],r,dd_gauss,None,[ex_4pdeer,None],uq=None)

assert ovl(P,fit.P) > 0.9
# ======================================================================
Expand Down Expand Up @@ -542,7 +542,7 @@ def test_global_scale_4pdeer():
t2 = np.linspace(0,5,250)
V2 = scales[1]*dipolarkernel(t2,r,pathways,Bmodel)@P

fit = fitsignal([V1,V2],[t1,t2],r,'P',bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal([V1,V2],[t1,t2],r,'P',bg_exp,ex_4pdeer,uq=None)

assert max(abs(np.asarray(scales)/np.asarray(fit.scale) - 1)) < 1e-2
# ======================================================================
Expand All @@ -559,7 +559,7 @@ def test_V_scale_parametric():
scale = 1.54e6
V = scale*(K@P)

fit = fitsignal(V,t,r,dd_gauss,bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal(V,t,r,dd_gauss,bg_exp,ex_4pdeer,uq=None)

assert max(abs(1 - V/fit.V)) < 1e-4
# ======================================================================
Expand All @@ -575,7 +575,7 @@ def test_V_scale():
scale =1.54e6
V = scale*(K@P)

fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uq=None)

assert max(abs(1 - V/fit.V)) < 1e-4
# ======================================================================
Expand All @@ -590,7 +590,7 @@ def test_V_scale_regularized():
scale = 1.54e6
V = scale*(K@P)

fit = fitsignal(V,t,r,'P',None,None,uqanalysis=False)
fit = fitsignal(V,t,r,'P',None,None,uq=None)

assert max(abs(1 - V/fit.V)) < 1e-4
# ======================================================================
Expand All @@ -606,7 +606,7 @@ def test_plot():
K = dipolarkernel(t,r,0.4,Bmodel)
V = K@P

fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uq=None)

fig = fit.plot(show=False)
assert str(fig.__class__)=="<class 'matplotlib.figure.Figure'>"
Expand All @@ -625,7 +625,7 @@ def test_physical_bg_model():
V = K@P
V = V0*V

fit = dl.fitsignal(V,t,r,'P',dl.bg_hom3d,dl.ex_4pdeer,uqanalysis=False)
fit = dl.fitsignal(V,t,r,'P',dl.bg_hom3d,dl.ex_4pdeer,uq=None)

assert abs(fit.bgparam - 50)<1e-1 and abs(fit.exparam - 0.4)<1e-1
# ======================================================================
Expand All @@ -642,7 +642,7 @@ def test_phenomenological_bg_model():
V = K@P
V = V0*V

fit = dl.fitsignal(V,t,r,'P',dl.bg_exp,dl.ex_4pdeer,uqanalysis=False)
fit = dl.fitsignal(V,t,r,'P',dl.bg_exp,dl.ex_4pdeer,uq=None)

assert abs(fit.bgparam - 0.3)<1e-1 and abs(fit.exparam - 0.4)<1e-1
# ======================================================================
Expand All @@ -660,7 +660,7 @@ def test_Vunmod():
Bscaled = (1-lam)*dl.bg_hom3d(t,50,lam)
V = K@P

fit = dl.fitsignal(V,t,r,'P',dl.bg_exp,dl.ex_4pdeer,uqanalysis=False)
fit = dl.fitsignal(V,t,r,'P',dl.bg_exp,dl.ex_4pdeer,uq=None)

assert max(abs(Bscaled - fit.Vunmod))<1e-4
# ======================================================================
Expand All @@ -677,7 +677,7 @@ def test_cost_value():
K = dipolarkernel(t,r,0.4,Bmodel)
V = K@P

fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uqanalysis=False)
fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uq=None)

assert isinstance(fit.cost,float) and np.round(fit.cost/np.sum(fit.residuals**2),5)==1
# ======================================================================
Expand Down