diff --git a/deerlab/fitsignal.py b/deerlab/fitsignal.py index f1e66978d..61466142f 100644 --- a/deerlab/fitsignal.py +++ b/deerlab/fitsignal.py @@ -16,7 +16,7 @@ def fitsignal(Vexp, t, r, dd_model='P', bg_model=bg_hom3d, ex_model=ex_4pdeer, dd_par0=None, bg_par0=None, ex_par0=None, verbose=False, dd_lb=None, bg_lb=None, ex_lb=None, dd_ub=None, bg_ub=None, ex_ub=None, - weights=1, uqanalysis=True, regparam='aic', regtype = 'tikhonov'): + weights=1, uqanalysis=True, uq='covariance', regparam='aic', regtype = 'tikhonov'): r""" Fits a dipolar model to the experimental signal ``V`` with time axis ``t``, using distance axis ``r``. The model is specified by the distance distribution (dd), @@ -81,6 +81,14 @@ def fitsignal(Vexp, t, r, dd_model='P', bg_model=bg_hom3d, ex_model=ex_4pdeer, If a model does not require parameters or are to be determined automatically it can be omitted or specified as ``None`` (default). + uq : string or list, optional + Type of uncertainty quantification analysis. Any ``UncertQuant`` output returned by this function will + be adjusted accordingly. The options are: + + * ``'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]``. + weights : array_like, optional Array of weighting coefficients for the individual signals in global fitting, the default is all weighted equally. @@ -239,6 +247,18 @@ def fitsignal(Vexp, t, r, dd_model='P', bg_model=bg_hom3d, ex_model=ex_4pdeer, if len(ex_model)!=nSignals: ex_model = ex_model*nSignals + # Default bootstrap samples + bootsamples = 1000 + if isinstance(uq, str): + uq = [uq] + if uq[0]!='bootstrap' and uq[0]!='covariance': + raise KeyError("Uncertainty quantification must be either 'covariance' or 'bootstrap'.") + + if uq[0]=='bootstrap': + # OVerride default if user has specified bootstraped samples + if len(uq)>1: bootsamples = uq[1] + uq = uq[0] + # 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]] lb = [[] if lb_i is None else lb_i for lb_i in [dd_lb,bg_lb,ex_lb]] @@ -334,11 +354,11 @@ def multiPathwayModel(par): K_ = dl.dipolarkernel(t[iSignal],r,pathways,Bfcn) Ks.append(K_) Bs.append(B_) - + return Ks, Bs # ========================================================================= - def splituq(full_uq,scales,Kfit=None): + def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1): # ========================================================================= """ Uncertainty quantification @@ -392,6 +412,11 @@ def splituq(full_uq,scales,Kfit=None): # ---------------------------------- nonneg = np.zeros_like(r) if parametricDistribution: + # Prepare parametric model + if includeForeground: + 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: subcovmat = covmat[np.ix_(Pfreeidx,Pfreeidx)] @@ -455,14 +480,39 @@ def splituq(full_uq,scales,Kfit=None): else: Vmod_uq.append([None]) - return Vfit_uq,Pfit_uq,Bfit_uq,Vmod_uq,Vunmod_uq,paruq_bg,paruq_ex,paruq_dd + return Vfit_uq, Pfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd # ========================================================================= - OnlyRegularization = np.all(~parametricDistribution & ~includeExperiment & ~includeBackground) - OnlyParametric = not OnlyRegularization and (parametricDistribution or not includeForeground) + def calculate_Vmod_Vunmod(parfit,Vfit,Bfit,scales): + # ========================================================================= + " Calculation of the (un)modulated components of the dipolar signal" + + # Calculate the unmodulated contribution (Vunmod) + # -------------------------------------------------------- + Vunmod = [] + for j in range(nSignals): + if includeExperiment[j]: + Lam0 = ex_model[j](parfit[exidx[j]])[0][0] + if includeBackground[j]: + Vunmod.append(Lam0*np.array(Bfit[j])) + else: + Vunmod.append(np.full_like(t[j],scales[j]*Lam0)) + else: + Vunmod.append(np.zeros_like(t[j])) + + # Calculate the modulated contribution (Vmod) + # -------------------------------------------------------- + Vmod = [] + for j in range(nSignals): + Vmod.append(Vfit[i] - Vunmod[i]) - if OnlyRegularization: + return Vmod, Vunmod + # ========================================================================= + + def regularization_analysis(Vexp): + # ========================================================================= + " Analysis workflow for non-parametric models based on regularized least-squares" # Use basic dipolar kernel Ks = [dl.dipolarkernel(ts,r) for ts in t] @@ -472,20 +522,27 @@ def splituq(full_uq,scales,Kfit=None): Pfit = fit.P Pfit_uq = fit.uncertainty scales = np.atleast_1d(fit.scale) - alphaopt = fit.regparam # Get fitted models Vfit = [scale*K@Pfit for K,scale in zip(Ks,scales)] Bfit = [scale*np.ones_like(V) for V,scale in zip(Vexp,scales)] + Vmod, Vunmod = calculate_Vmod_Vunmod(None,Vfit,Bfit,scales) # No parameters - parfit_ = np.asarray([None]) - if uqanalysis: - Vfit_uq, Pfit_uq, Bfit_uq,Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(Pfit_uq,scales,Ks) + 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) + 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 + # ========================================================================= + + def nonlinear_lsq_analysis(Vexp): + # ========================================================================= + " Analysis workflow for fully parametric models based on nonlinear least-squares" - elif OnlyParametric: - # Prepare the full-parametric model if includeForeground: Pfcn = lambda par: dd_model(r,par[ddidx]) @@ -495,30 +552,37 @@ def splituq(full_uq,scales,Kfit=None): # Non-linear parametric fit fit = dl.fitparamodel(Vexp,Vmodel,par0,lb,ub,weights=weights,uqanalysis=uqanalysis) - parfit_ = fit.param + parfit = fit.param param_uq = fit.uncertainty - scales = fit.scale + scales = np.atleast_1d(fit.scale) alphaopt = None # Get fitted models - Vfit = Vmodel(parfit_) - _,Bfit = multiPathwayModel(parfit_) + Vfit = Vmodel(parfit) + _,Bfit = multiPathwayModel(parfit) if includeForeground: - Pfit = Pfcn(parfit_) + Pfit = Pfcn(parfit) else: Pfit = [] if type(Vfit) is not list: Vfit = [Vfit] - if type(scales) is not list: - scales = [scales] + if type(Bfit) is not list: + Bfit = [Bfit] Bfit = [scale*B for B,scale in zip(Bfit,scales)] - Vfit = [V*scale for scale,V in zip(scales,Vfit) ] + Vfit = [scale*V for V,scale in zip(Vfit,scales) ] + 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) + 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 + # ========================================================================= + + def separable_nonlinear_lsq_analysis(Vexp): + # ========================================================================= + " Analysis workflow for semiparametric models based on separable nonlinear least-squares" - if uqanalysis: - Vfit_uq, Pfit_uq, Bfit_uq,Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(param_uq,scales) - - else: - # Non-negativity constraint on distributions lbl = np.zeros_like(r) @@ -528,41 +592,78 @@ def splituq(full_uq,scales,Kfit=None): # Separable non-linear least squares (SNNLS) fit = dl.snlls(Vexp_,lambda par: multiPathwayModel(par)[0],par0,lb,ub,lbl, reg=True, regparam=regparam, uqanalysis=uqanalysis, weights=weights) - parfit_ = fit.nonlin + parfit = fit.nonlin Pfit = fit.lin snlls_uq = fit.uncertainty alphaopt = fit.regparam scales = [prescales[i]*np.trapz(Pfit,r) for i in range(nSignals)] # Get the fitted models - Kfit,Bfit = multiPathwayModel(parfit_) + Kfit,Bfit = multiPathwayModel(parfit) Bfit = [scale*B for B,scale in zip(Bfit,scales)] Vfit = [scale*K@Pfit for K,scale in zip(Kfit,scales)] + Vmod, Vunmod = calculate_Vmod_Vunmod(parfit,Vfit,Bfit,scales) - if uqanalysis: - Vfit_uq, Pfit_uq, Bfit_uq,Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(snlls_uq,scales,Kfit) - - # Calculate the unmodulated contribution (Vunmod) - # -------------------------------------------------------- - Vunmod = [] - for j in range(nSignals): - if includeExperiment[j]: - Lam0 = ex_model[j](parfit_[exidx[j]])[0][0] - if includeBackground[j]: - Vunmod.append(Lam0*np.array(Bfit[j]) ) - else: - print(ex_model[j](parfit_[exidx[j]])) - print(scales) - Vunmod.append(np.full_like(t[j],scales[j]*Lam0)) + if uqanalysis and uq=='covariance': + Vfit_uq, Pfit_uq, Bfit_uq, Vmod_uq, Vunmod_uq, paruq_bg, paruq_ex, paruq_dd = splituq(snlls_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: - Vunmod.append(np.zeros_like(t[j])) + return fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit, scales, alphaopt + # ========================================================================= - - # Calculate the modulated contribution (Vmod) - # -------------------------------------------------------- - Vmod = [] - for j in range(nSignals): - Vmod.append(Vfit[i] - Vunmod[i]) + # Analyze the data + # ---------------------- + + # Determine type of model + nonparametric = np.all(~parametricDistribution & ~includeExperiment & ~includeBackground) + fullparametric = not nonparametric and (parametricDistribution or not includeForeground) + semiparametric = not nonparametric and not fullparametric + + # Choose appropiate analysis for type of model + if nonparametric: + analysis = regularization_analysis + elif fullparametric: + analysis = nonlinear_lsq_analysis + elif semiparametric: + analysis = separable_nonlinear_lsq_analysis + + # Run the analysis + results = analysis(Vexp) + + # Unpack results + if uqanalysis and uq=='covariance': + 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 = results + else: + fit, Pfit, Vfit, Bfit, Vmod, Vunmod, parfit_, scales, alphaopt = results + + # Bootstrapping uncertainty quantification + # ----------------------------------------- + if uqanalysis and uq=='bootstrap': + + def bootstrapfcn(Vexp): + # ====================================================== + # Fit the data + _, Pfit_, Vfit_, Bfit_, Vmod_, Vunmod_, parfit, _, _ = analysis(Vexp) + # Extract the individual parameter subsets + parfit_bg = [parfit[bgidx[n]] for n in range(nSignals)] + parfit_ex = [parfit[exidx[n]] for n in range(nSignals)] + parfit_dd = parfit[ddidx] + + return Pfit_,*Vfit_,*Bfit_,*Vmod_,*Vunmod_,*parfit_bg,*parfit_ex,parfit_dd + # ====================================================== + + # Run bootstrapping + boot_uq = dl.bootan(bootstrapfcn,Vexp,Vfit,samples=bootsamples,verbose=verbose) + + # Unpack bootstrapping results + Pfit_uq = boot_uq[0] + Vfit_uq = [boot_uq[1+n] for n in range(nSignals)] + Bfit_uq = [boot_uq[1+nSignals+n] for n in range(nSignals)] + Vmod_uq = [boot_uq[1+2*nSignals+n] for n in range(nSignals)] + Vunmod_uq = [boot_uq[1+3*nSignals+n] for n in range(nSignals)] + paruq_bg = [boot_uq[1+4*nSignals+n] for n in range(nSignals)] + paruq_ex = [boot_uq[1+5*nSignals+n] for n in range(nSignals)] + paruq_dd = boot_uq[-1] # Normalize distribution # ----------------------- diff --git a/test/test_fitsignal.py b/test/test_fitsignal.py index 2d66362bc..a5101f1a0 100644 --- a/test/test_fitsignal.py +++ b/test/test_fitsignal.py @@ -332,28 +332,30 @@ def assert_confidence_intervals(pci50,pci95,pfit,lb,ub): assert not errors, "Errors occured:\n{}".format("\n".join(errors)) #---------------------------------------------------------------------- -def assert_confinter_param(subset): #---------------------------------------------------------------------- - exmodel = ex_4pdeer - ddmodel = dd_gauss - bgmodel = bg_exp +exmodel = ex_4pdeer +ddmodel = dd_gauss +bgmodel = bg_exp - r = np.linspace(2,6,40) - P = ddmodel(r,[4.5, 0.25]) +r = np.linspace(2,6,40) +P = ddmodel(r,[4.5, 0.25]) - info = exmodel() - parIn = info['Start'] - pathways = exmodel(parIn) +info = exmodel() +parIn = info['Start'] +pathways = exmodel(parIn) - kappa = 0.4 - Bmodel = lambda t,lam: bgmodel(t,kappa) +kappa = 0.4 +Bmodel = lambda t,lam: bgmodel(t,kappa) - t = np.linspace(0,5,100) - np.random.seed(0) - V = dipolarkernel(t,r,pathways,Bmodel)@P + whitegaussnoise(t,0.01) - - fit = fitsignal(V,t,r,ddmodel,bgmodel,exmodel,uqanalysis=True) +t = np.linspace(0,5,100) +np.random.seed(0) +V = dipolarkernel(t,r,pathways,Bmodel)@P + whitegaussnoise(t,0.01) +fit = fitsignal(V,t,r,ddmodel,bgmodel,exmodel,uqanalysis=True) +#---------------------------------------------------------------------- + +def assert_confinter_param(subset): +#---------------------------------------------------------------------- if subset == 'ex': info = exmodel() pfit = fit.exparam @@ -394,33 +396,37 @@ def test_confinter_ddparam(): assert_confinter_param('dd') # ====================================================================== -def assert_confinter_models(subset): #---------------------------------------------------------------------- - exmodel = ex_4pdeer - if subset == 'Pfitfree': - ddmodel = 'P' - subset = 'Pfit' - else: - ddmodel= dd_gauss - bgmodel = bg_exp +exmodel = ex_4pdeer +bgmodel = bg_exp - r = np.linspace(2,6,40) - P = dd_gauss(r,[4.5, 0.25]) +r = np.linspace(2,6,40) +P = dd_gauss(r,[4.5, 0.25]) - info = exmodel() - parIn = info['Start'] - pathways = exmodel(parIn) +info = exmodel() +parIn = info['Start'] +pathways = exmodel(parIn) - kappa = 0.4 - Bmodel = lambda t: bgmodel(t,kappa) +kappa = 0.4 +Bmodel = lambda t: bgmodel(t,kappa) - t = np.linspace(0,5,100) - np.random.seed(0) - V = dipolarkernel(t,r,pathways,Bmodel)@P + whitegaussnoise(t,0.03) - - fit = fitsignal(V,t,r,ddmodel,bgmodel,exmodel,uqanalysis=True) +t = np.linspace(0,5,100) +np.random.seed(0) +V = dipolarkernel(t,r,pathways,Bmodel)@P + whitegaussnoise(t,0.03) + +fit_Pparam = fitsignal(V,t,r,ddmodel,bgmodel,exmodel,uq='covariance') +fit_Pfree = fitsignal(V,t,r,'P',bgmodel,exmodel,uq='covariance') +#---------------------------------------------------------------------- - if subset == 'Pfit': +def assert_confinter_models(subset): +#---------------------------------------------------------------------- + + if subset=='Pfitfree': + fit = fit_Pfree + else: + fit = fit_Pparam + + if subset == 'Pfit' or subset == 'Pfitfree': modelfit = fit.P lb = np.zeros_like(r) ub = np.full_like(r,inf) @@ -466,6 +472,7 @@ def test_confinter_Bfit(): assert_confinter_models('Bfit') # ====================================================================== + def assert_confinter_noforeground(): # ====================================================================== "Check that the confidence inervals for a pure background fit are correct" @@ -674,3 +681,96 @@ def test_cost_value(): assert isinstance(fit.cost,float) and np.round(fit.cost/np.sum(fit.residuals**2),5)==1 # ====================================================================== + +# ---------------------------------------------------------------------- + exmodel = ex_4pdeer + bgmodel = bg_exp + + r = np.linspace(2,6,40) + P = dd_gauss(r,[4.5, 0.25]) + + info = exmodel() + parIn = info['Start'] + pathways = exmodel(parIn) + + kappa = 0.4 + Bmodel = lambda t: bgmodel(t,kappa) + + t = np.linspace(0,5,100) + np.random.seed(0) + V = dipolarkernel(t,r,pathways,Bmodel)@P + whitegaussnoise(t,0.03) + + fit = fitsignal(V,t,r,ddmodel,bgmodel,exmodel,uq=['bootstrap',2]) +# ---------------------------------------------------------------------- + +def assert_boot_ci(quantity): +# ---------------------------------------------------------------------- + + if quantity=='P': + ci = fit.Puncert.ci(95) + elif quantity=='V': + ci = fit.Vuncert.ci(95) + elif quantity=='Vmod': + ci = fit.VmodUncert.ci(95) + elif quantity=='Vunmod': + ci = fit.VunmodUncert.ci(95) + elif quantity=='B': + ci = fit.Buncert.ci(95) + elif quantity=='ddparam': + ci = fit.ddparamUncert.ci(95) + elif quantity=='bgparam': + ci = fit.bgparamUncert.ci(95) + elif quantity=='exparam': + ci = fit.exparamUncert.ci(95) + + assert np.all(ci[:,0]<=ci[:,1]) +# ---------------------------------------------------------------------- + + +def test_bootci_P(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('P') +# ====================================================================== + +def test_bootci_V(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('V') +# ====================================================================== + +def test_bootci_B(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('B') +# ====================================================================== + +def test_bootci_Vmod(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('Vmod') +# ====================================================================== + +def test_bootci_Vunmod(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('Vunmod') +# ====================================================================== + +def test_bootci_ddparam(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('ddparam') +# ====================================================================== + +def test_bootci_exparam(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('exparam') +# ====================================================================== + +def test_bootci_bparam(): +# ====================================================================== + "Check that the bootstrapped confidence intervals work" + assert_boot_ci('bgparam') +# ====================================================================== \ No newline at end of file