From fdb09af9c472da129c4179625e9469df8907524e Mon Sep 17 00:00:00 2001 From: Luis Fabregas <48292540+luisfabib@users.noreply.github.com> Date: Wed, 24 Mar 2021 21:14:27 +0100 Subject: [PATCH] Implement arbitrary penalties for SNLLS & Pathway amplitudes contstraint (#108) * snlls: remove lines leading to double smoothness penalization * snlls: redesign of residual function and uncertainty quantification - add option to pass a custom penalty function - snlls new returns separate uncertainty quantification objects for the linear and nonlinear components - adapted fitisignal to new outputs of snlls * fitmodel: implement constraint to ensure identifiability of Lam0 and V0 * snlls: correct behaviour of the uncertainty quantification * tests: relax tight assertion for one test * hccm: fix minor bug in internal function * fitmodel: remove deprecated keyword argument * snlls: rename keyword custom_penalty to extrapenalty, fix some tests * snlls: code refactoring * fitmodel: fix merge error * fitmodel: fix another merge error * lsqcomponents: refactor code --- deerlab/fitmodel.py | 66 ++++++++---- deerlab/lsqcomponents.py | 5 +- deerlab/snlls.py | 212 ++++++++++++++++++++++++--------------- deerlab/utils/utils.py | 4 +- test/test_fitmodel.py | 10 +- test/test_snlls.py | 24 ++--- 6 files changed, 198 insertions(+), 123 deletions(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index 5eac545c2..883c3fd5d 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -16,7 +16,7 @@ def fitmodel(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, uq='covariance', regparam='aic', regtype = 'tikhonov', + weights=1, uq='covariance', regparam='aic', regtype = 'tikhonov', tol=1e-10,maxiter=1e8): r""" Fits a dipolar model to the experimental signal ``V`` with time axis ``t``, using @@ -377,13 +377,16 @@ 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 - covmat = full_uq.covmat + if isinstance(full_uq,list): + paruq = full_uq[0] + else: + paruq = full_uq + covmat = paruq.covmat Nparam = len(parfit_) paramidx = np.arange(Nparam) - Pfreeidx = np.arange(Nparam,np.shape(covmat)[0]) # Full parameter set uncertainty # ------------------------------- @@ -427,8 +430,7 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1): 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)] - Pfit_uq = UncertQuant('covariance',Pfit,subcovmat,nonneg,[]) + Pfit_uq = full_uq[1] # Background uncertainty # ----------------------- @@ -465,10 +467,7 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1): elif includeForeground: # Parametric signal with parameter-free distribution Vmodel = lambda par: scales[jj]*multiPathwayModel(par[paramidx])[0][jj]@Pfit - Jnonlin = Jacobian(Vmodel,parfit_,lb[paramidx],ub[paramidx]) - J = np.concatenate((Jnonlin, scales[jj]*Kfit[jj]),1) - Vcovmat = J@covmat@J.T - Vfit_uq.append( UncertQuant('covariance',Vfit[jj],Vcovmat)) + Vfit_uq.append( paruq.propagate(Vmodel) ) else: Vfit_uq.append([None]) @@ -477,14 +476,7 @@ def splituq(full_uq,Pfit,Vfit,Bfit,parfit_,Kfit,scales=1): for jj in range(nSignals): if includeForeground: Vmod_fcn = lambda par: Vmodel(par) - Vunmod_fcn(par) - Vunmod_uq.append( paruq.propagate(lambda par:Vmod_fcn(par)) ) - Jnonlin = Jacobian(Vmod_fcn,parfit_,lb[paramidx],ub[paramidx]) - if parametricDistribution: - J = Jnonlin - else: - J = np.concatenate((Jnonlin, scales[jj]*Kfit[jj]),1) - Vmod_covmat = J@covmat@J.T - Vmod_uq.append( UncertQuant('covariance',Vmod_fcn(parfit_),Vmod_covmat)) + Vmod_uq.append( paruq.propagate(Vmod_fcn) ) else: Vmod_uq.append([None]) @@ -587,6 +579,7 @@ def nonlinear_lsq_analysis(Vexp): 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" @@ -594,18 +587,47 @@ def separable_nonlinear_lsq_analysis(Vexp): # Non-negativity constraint on distributions lbl = np.zeros_like(r) - prescales = [max(V) for V in Vexp] + #prescales = [max(V) for V in Vexp] + prescales = [1 for V in Vexp] Vexp_ = [Vexp[i]/prescales[i] for i in range(nSignals)] + def scale_constraint(nonlinpar): + # -------------------------------------------------------- + penalty = np.zeros(nSignals) + for i in range(nSignals): + ex_par = nonlinpar[exidx[i]] + pathways = ex_fcn[i](ex_par) + lams = [pathway[0] for pathway in pathways] + if np.sum(lams)<1 or ex_model[i].__name__=='ex_4pdeer': + penalty[i] = 0 + else: + penalty[i] = max(Vexp[i])*(np.sum(lams) - 1) + + return penalty + # -------------------------------------------------------- + + # 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, + regparam=regparam, uqanalysis=uqanalysis, weights=weights,extrapenalty=scale_constraint, nonlin_tol=tol,nonlin_maxiter=maxiter) parfit = fit.nonlin Pfit = fit.lin - snlls_uq = fit.uncertainty + param_uq = fit.nonlinUncert + Pfit_uq = fit.linUncert + snlls_uq = [param_uq,Pfit_uq] alphaopt = fit.regparam - scales = np.atleast_1d([prescales[i]*np.trapz(Pfit,r) for i in range(nSignals)]) + scales = fit.scale + + # Normalize distribution + # ----------------------- + Pscale = np.trapz(Pfit,r) + Pfit /= Pscale + scales = np.atleast_1d([scale*Pscale for scale in scales]) + if uqanalysis and uq=='covariance': + # scale CIs accordingly + Pfit_uq_ = copy.deepcopy(Pfit_uq) # need a copy to avoid infite recursion on next step + Pfit_uq.ci = lambda p: Pfit_uq_.ci(p)/Pscale # Get the fitted models Kfit,Bfit = multiPathwayModel(parfit) diff --git a/deerlab/lsqcomponents.py b/deerlab/lsqcomponents.py index de79e806c..d2123509f 100644 --- a/deerlab/lsqcomponents.py +++ b/deerlab/lsqcomponents.py @@ -1,6 +1,6 @@ import numpy as np -def lsqcomponents(V,K,L,alpha,weights, regtype = 'tikh',huberparam = 1.35): +def lsqcomponents(V,K,L=None,alpha=0,weights=1, regtype = 'tikh',huberparam = 1.35): # ============================================================================================== # Prepare @@ -10,6 +10,9 @@ def lsqcomponents(V,K,L,alpha,weights, regtype = 'tikh',huberparam = 1.35): KtV = weights*K.T@V KtK = weights*K.T@K + if L is None: + return KtK, KtV + def optimizeregterm(fun,threshold): # ============================================================================ P = np.zeros(nr) diff --git a/deerlab/snlls.py b/deerlab/snlls.py index 95417a25e..27352f957 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -16,7 +16,7 @@ from deerlab.classes import UncertQuant, FitResult def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx', reg='auto', weights=1, - regtype='tikhonov', regparam='aic', multistart=1, regorder=2, alphareopt=1e-3, + regtype='tikhonov', regparam='aic', multistart=1, regorder=2, alphareopt=1e-3, extrapenalty=None, nonlin_tol=1e-9, nonlin_maxiter=1e8, lin_tol=1e-15, lin_maxiter=1e4, huberparam=1.35, uqanalysis=True): r""" Separable Non-linear Least Squares Solver @@ -76,6 +76,10 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx The regularization parameter can be manually specified by passing a scalar value instead of a string. The default ``'aic'``. + custom_penalty: callable + Custom penalty function to impose upon the solution. Must return a vector to be + added to the residual vector. + alphareopt : float scalar, optional Relative parameter change threshold for reoptimizing the regularization parameter when using a selection method, the default is 1e-3. @@ -118,14 +122,10 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx Fitted non-linear parameters lin : ndarray Fitted linear parameters - paramuq : :ref:`UncertQuant` - Uncertainty quantification of the joined parameter - set (linear + non-linear parameters). The confidence intervals - of the individual subsets can be requested via: - - * ``paramuq.ci(n)`` - n%-CI of the full parameter set - * ``paramuq.ci(n,'lin')`` - n%-CI of the linear parameter set - * ``paramuq.ci(n,'nonlin')`` - n%-CI of the non-linear parameter set + nonlinUncert : :ref:`UncertQuant` + Uncertainty quantification of the non-linear parameter set. + linUncert : :ref:`UncertQuant` + Uncertainty quantification of the linear parameter set. regparam : scalar Regularization parameter value used for the regularization of the linear parameters. plot : callable @@ -193,21 +193,21 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx par0 = np.atleast_1d(par0) # Parse multiple datsets and non-linear operators into a single concatenated vector/matrix - y, Amodel, weights, subsets, _ = dl.utils.parse_multidatasets(y, Amodel, weights, precondition=True) + y, Amodel, weights, subsets = dl.utils.parse_multidatasets(y, Amodel, weights, precondition=False) # Get info on the problem parameters and non-linear operator A0 = Amodel(par0) Nnonlin = len(par0) Nlin = np.shape(A0)[1] linfit = np.zeros(Nlin) - alpha = 0 - + scales = [1 for _ in subsets] + prescales= [1 for _ in subsets] # Determine whether to use regularization penalty illConditioned = np.linalg.cond(A0) > 10 if reg == 'auto': - includePenalty = illConditioned + includeRegularization = illConditioned else: - includePenalty = reg + includeRegularization = reg # Checks for bounds constraints # ---------------------------------------------------------- @@ -247,14 +247,14 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx nonNegativeOnly = (np.all(lbl == 0)) and (np.all(np.isinf(ubl))) - if includePenalty: - # Use an arbitrary axis - ax = np.arange(1, Nlin+1) + # Use an arbitrary axis + ax = np.arange(1, Nlin+1) + if includeRegularization : # Get regularization operator regorder = np.minimum(Nlin-1, regorder) L = dl.regoperator(ax, regorder) else: - L = np.eye(Nlin, Nlin) + L = None # Prepare the linear solver # ---------------------------------------------------------- @@ -283,6 +283,31 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx check = False regparam_prev = 0 par_prev = [0]*len(par0) + alpha = None + + def linear_problem(A,optimize_alpha,alpha): + #=========================================================================== + """ + Linear problem + ------------------ + Solves the linear subproblem of the SNLLS objective function via linear LSQ + constrained, unconstrained, or regularized. + """ + + # Optimiza the regularization parameter only if needed + if optimize_alpha: + alpha = dl.selregparam(y, A, ax, regtype, regparam, regorder=regorder) + + # Components for linear least-squares + AtA, Aty = dl.lsqcomponents(y, A, L, alpha, weights, regtype=regtype) + + # Solve the linear least-squares problem + result = linSolver(AtA, Aty) + linfit = parseResult(result) + linfit = np.atleast_1d(linfit) + + return linfit, alpha + #=========================================================================== def ResidualsFcn(p): #=========================================================================== @@ -293,54 +318,68 @@ def ResidualsFcn(p): non-linear least-squares solver. """ - nonlocal par_prev, check, regparam_prev, linfit, alpha + nonlocal par_prev, check, regparam_prev, scales, linfit, alpha + # Non-linear model evaluation A = Amodel(p) - # Regularization components - if includePenalty: + + # Check whether optimization of the regularization parameter is needed + if includeRegularization : if type(regparam) is str: # If the parameter vector has not changed by much... if check and all(abs(par_prev-p)/p < alphareopt): # ...use the alpha optimized in the previous iteration + optimize_alpha = False alpha = regparam_prev else: # ...otherwise optimize with current settings - alpha = dl.selregparam(y, A, ax, regtype, regparam, regorder=regorder) + alpha = regparam + optimize_alpha = True check = True else: # Fixed regularization parameter alpha = regparam - + optimize_alpha = False # Store current iteration data for next one par_prev = p regparam_prev = alpha - else: # Non-linear operator without penalty + optimize_alpha = False alpha = 0 - # Components for linear least-squares - AtA, Aty = dl.lsqcomponents(y, A, L, alpha, weights, regtype=regtype) + linfit,alpha = linear_problem(A,optimize_alpha,alpha) + regparam_prev = alpha - # Solve the linear least-squares problem - result = linSolver(AtA, Aty) - linfit = parseResult(result) - linfit = np.atleast_1d(linfit) # Evaluate full model residual yfit = A@linfit + + # Optimize the scale yfit + scales, scales_vec = [], np.zeros_like(yfit) + for subset in subsets: + yfit_,y_ = (np.atleast_2d(y[subset]) for y in [yfit, y]) # Rescale the subsets corresponding to each signal + scale = np.squeeze(np.linalg.lstsq(yfit_.T,y_.T,rcond=None)[0]) + scales.append(scale) # Store the optimized scales of each signal + scales_vec[subset] = scale + # Compute residual vector - res = weights*(yfit - y) - if includePenalty: - penalty = alpha*L@linfit - # Augmented residual - res = np.concatenate((res, penalty)) - res, _ = _augment(res, [], regtype, alpha, L, linfit, huberparam, Nnonlin) + res = weights*(scales_vec*(Amodel(p)@linfit) - y) + # Compute residual from custom penalty + if callable(extrapenalty): + penres = extrapenalty(p) + penres = np.atleast_1d(penres) + res = np.concatenate((res,penres)) + + if includeRegularization : + # Augmented residual + res_reg, _ = reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin) + res = np.concatenate((res,res_reg)) + return res #=========================================================================== - # Preprare multiple start global optimization if requested if multistart > 1 and not nonLinearConstrained: raise TypeError('Multistart optimization cannot be used with unconstrained non-linear parameters.') @@ -357,63 +396,59 @@ def ResidualsFcn(p): linfits.append(linfit) fvals.append(2*sol.cost) # least_squares uses 0.5*sum(residual**2) sols.append(sol) + # Find global minimum from multiple runs globmin = np.argmin(fvals) linfit = linfits[globmin] nonlinfit = nonlinfits[globmin] sol = sols[globmin] Afit = Amodel(nonlinfit) - yfit = Afit@linfit + + scales_vec = np.zeros_like(y) + for subset,scale in zip(subsets,scales): + scales_vec[subset] = scale + yfit = scales_vec*(Afit@linfit) # Uncertainty analysis - #-------------------------------------------------------- + #--------------------- if uqanalysis: - # Compue the residual vector - res = weights*(yfit - y) - - # Compute the Jacobian for the linear and non-linear parameters - fcn = lambda p: Amodel(p)@linfit - Jnonlin = Jacobian(fcn,nonlinfit,lb,ub) + # Compute the fit residual + res = ResidualsFcn(nonlinfit) + + # Jacobian (non-linear part) + Jnonlin = Jacobian(ResidualsFcn,nonlinfit,lb,ub) - Jlin = Afit - J = np.concatenate((Jnonlin, Jlin),1) + # Jacobian (linear part) + Jlin = np.zeros((len(res),len(linfit))) + Jlin[:len(y),:] = Amodel(nonlinfit) + Jlin[len(res)-Nlin:,:] = reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin)[1] - # Augment the residual and Jacobian with the regularization penalty on the linear parameters - res, J = _augment(res, J, regtype, regparam_prev, L, linfit, huberparam, Nnonlin) + # Full Jacobian + J = np.concatenate((Jnonlin,Jlin),axis=1) # Calculate the heteroscedasticity consistent covariance matrix covmatrix = hccm(J, res, 'HC1') - + # Get combined parameter sets and boundaries parfit = np.concatenate((nonlinfit, linfit)) lbs = np.concatenate((lb, lbl)) ubs = np.concatenate((ub, ubl)) # Construct the uncertainty quantification object - paramuq_ = UncertQuant('covariance', parfit, covmatrix, lbs, ubs) - paramuq = copy.deepcopy(paramuq_) - - def ci(coverage,ptype='full'): - #=========================================================================== - "Wrapper around the CI function handle of the uncertainty structure" - # Get requested confidence interval of joined parameter set - paramci = paramuq_.ci(coverage) - if ptype == 'nonlin': - # Return only confidence intervals on non-linear parameters - paramci = paramci[range(Nnonlin), :] - elif ptype == 'lin': - # Return only confidence intervals on linear parameters - paramci = paramci[Nnonlin:, :] - return paramci - #=========================================================================== - - # Add the function to the confidence interval function call - paramuq.ci = ci + paramuq = UncertQuant('covariance', parfit, covmatrix, lbs, ubs) + + # Split the uncertainty quantification of nonlinear/linear parts + nonlin_subset = np.arange(0,Nnonlin) + lin_subset = np.arange(Nnonlin,Nnonlin+Nlin) + paramuq_nonlin = uq_subset(paramuq,nonlin_subset) + paramuq_lin = uq_subset(paramuq,lin_subset) + else: - paramuq = [] + paramuq_nonlin = [] + paramuq_lin = [] # Goodness-of-fit - # -------------------------------------- + # --------------- stats = [] for subset in subsets: Ndof = len(y[subset]) - Nnonlin @@ -422,17 +457,38 @@ def ci(coverage,ptype='full'): stats = stats[0] fvals = fvals[0] + for i in range(len(subsets)): + scales[i] *= prescales[i] + # Display function def plotfcn(show=False): fig = _plot(subsets,y,yfit,show) return fig - return FitResult(nonlin=nonlinfit, lin=linfit, uncertainty=paramuq, regparam=alpha, plot=plotfcn, - stats=stats, cost=fvals, residuals=sol.fun, success=sol.success) + return FitResult(nonlin=nonlinfit, lin=linfit, nonlinUncert=paramuq_nonlin, linUncert=paramuq_lin, regparam=alpha, plot=plotfcn, + stats=stats, cost=fvals, residuals=sol.fun, success=sol.success, scale=scales) # =========================================================================================== -def _augment(res, J, regtype, alpha, L, x, eta, Nnonlin): +def uq_subset(uq_full,subset): +#=========================================================================== + "Wrapper around the CI function handle of the uncertainty structure" + uq_subset = copy.deepcopy(uq_full) + + uq_subset.mean = uq_subset.mean[subset] + uq_subset.median = uq_subset.median[subset] + uq_subset.std = uq_subset.std[subset] + uq_subset.covmat = uq_subset.covmat[np.ix_(subset,subset)] + uq_subset.nparam = len(subset) + + # Get requested confidence interval of joined parameter set + uq_subset.ci = lambda coverage: uq_full.ci(coverage)[subset, :] + uq_subset.percentile = lambda p: uq_full.percentile(p)[subset] + + return uq_subset +#=========================================================================== + +def reg_penalty(regtype, alpha, L, x, eta, Nnonlin): # =========================================================================================== """ LSQ residual and Jacobian augmentation @@ -458,13 +514,7 @@ def _augment(res, J, regtype, alpha, L, x, eta, Nnonlin): resreg = alpha*resreg Jreg = alpha*Jreg - # Augment jacobian and residual - res = np.concatenate((res, resreg)) - if np.size(J) != 0: - Jreg = np.concatenate((np.zeros((np.shape(L)[0],Nnonlin)), Jreg),1) - J = np.concatenate((J, Jreg)) - - return res, J + return resreg, Jreg # =========================================================================================== diff --git a/deerlab/utils/utils.py b/deerlab/utils/utils.py index 62f237cac..fe29f4677 100644 --- a/deerlab/utils/utils.py +++ b/deerlab/utils/utils.py @@ -190,14 +190,14 @@ def hccm(J,*args): elif mode.upper() == 'HC4': # Cribari-Neto,(2004),[4] # Compute discount factor - delta = min(4,n*h/k) + delta = np.minimum(4,n*h/k) # Estimate the data covariance matrix V = np.diag(res**2./((1 - h)**delta)) elif mode.upper() == 'HC5': # Cribari-Neto,(2007),[5] # Compute inflation factor k = 0.7 - alpha = min(max(4,k*max(h)/np.mean(h)),h/np.mean(h)) + alpha = np.minimum(np.maximum(4,k*max(h)/np.mean(h)),h/np.mean(h)) # Estimate the data covariance matrix V = np.diag(res**2./(np.sqrt((1 - h)**alpha))) diff --git a/test/test_fitmodel.py b/test/test_fitmodel.py index cfb908727..bd7dbd6fe 100644 --- a/test/test_fitmodel.py +++ b/test/test_fitmodel.py @@ -351,7 +351,7 @@ def assert_confidence_intervals(pci50,pci95,pfit,lb,ub): np.random.seed(0) V = dipolarkernel(t,r,pathways,Bmodel)@P + whitegaussnoise(t,0.01) -fit = fitmodel(V,t,r,ddmodel,bgmodel,exmodel,uqanalysis=True) +fit = fitmodel(V,t,r,ddmodel,bgmodel,exmodel,uq='covariance') #---------------------------------------------------------------------- def assert_confinter_param(subset): @@ -489,7 +489,7 @@ def assert_confinter_noforeground(): np.random.seed(0) V = dipolarkernel(t,r,lam,Bmodel)@P + whitegaussnoise(t,0.01) - fit = fitmodel(V,t,r,None,bgmodel,ex_4pdeer,uqanalysis=True) + fit = fitmodel(V,t,r,None,bgmodel,ex_4pdeer,uq='covariance') Bfit = fit.B Buq = fit.Buncert @@ -512,7 +512,7 @@ def assert_confinter_dipevofun(): np.random.seed(0) V = dipolarkernel(t,r)@P + whitegaussnoise(t,0.01) - fit = fitmodel(V,t,r,'P',None,None,uqanalysis=True) + fit = fitmodel(V,t,r,'P',None,None,uq='covariance') Pfit = fit.P Puq = fit.Puncert @@ -544,9 +544,9 @@ def test_global_scale_4pdeer(): fit = fitmodel([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 + assert max(abs(np.asarray(scales)/np.asarray(fit.scale) - 1)) < 5e-2 # ====================================================================== - +test_global_scale_4pdeer() def test_V_scale_parametric(): # ====================================================================== diff --git a/test/test_snlls.py b/test/test_snlls.py index 226a8e3ce..7acb761ae 100644 --- a/test/test_snlls.py +++ b/test/test_snlls.py @@ -1,5 +1,5 @@ import numpy as np -from deerlab import dipolarkernel,dd_gauss,dd_gauss2,snlls +from deerlab import dipolarkernel,dd_gauss,dd_gauss2,snlls,whitegaussnoise from deerlab.bg_models import bg_exp from deerlab.utils import ovl @@ -119,7 +119,7 @@ def assert_confidence_intervals(pci50,pci95,pfit,lb,ub): errors = [] if not np.all(p95lb <= pfit) and not np.all(p50lb <= pfit): errors.append("Some fitted values are below the lower bound of the confidence intervals.") - if not np.all(p95ub >= pfit) and not np.all(p50lb >= pfit): + if not np.all(p95ub >= pfit) and not np.all(p50ub >= pfit): errors.append("Some fitted values are over the upper bound of the confidence intervals.") if not np.all(p95lb <= p50lb): errors.append("The 50%-CI has lower values than the 95%-CI") @@ -136,13 +136,13 @@ def test_confinter_linear(): "Check that the confidence intervals of the linear parameters are correct" # Prepare test data - r = np.linspace(1,8,80) + r = np.linspace(1,8,150) t = np.linspace(0,4,200) lam = 0.25 K = dipolarkernel(t,r,lam) parin = [3.5, 0.4, 0.6, 4.5, 0.5, 0.4] P = dd_gauss2(r,parin) - V = K@P + V = K@P + whitegaussnoise(t,0.05,seed=1) # Non-linear parameters # nlpar = [lam] @@ -153,11 +153,11 @@ def test_confinter_linear(): lbl = np.zeros(len(r)) ubl = np.full(len(r), np.inf) # Separable LSQ fit - fit = snlls(V,lambda lam: dipolarkernel(t,r,lam),nlpar0,lb,ub,lbl,ubl) - Pfit = fit.lin - uq = fit.uncertainty - Pci50 = uq.ci(50,'lin') - Pci95 = uq.ci(95,'lin') + fit = snlls(V,lambda lam: dipolarkernel(t,r,lam),nlpar0,lb,ub,lbl) + Pfit = np.round(fit.lin,6) + uq = fit.linUncert + Pci50 = np.round(uq.ci(50),6) + Pci95 = np.round(uq.ci(95),6) assert_confidence_intervals(Pci50,Pci95,Pfit,lbl,ubl) #======================================================================= @@ -187,9 +187,9 @@ def test_confinter_nonlinear(): fit = snlls(V,lambda lam: dipolarkernel(t,r,lam),nlpar0,lb,ub,lbl,ubl) parfit = fit.nonlin - uq = fit.uncertainty - parci50 = uq.ci(50,'nonlin') - parci95 = uq.ci(95,'nonlin') + uq = fit.nonlinUncert + parci50 = uq.ci(50) + parci95 = uq.ci(95) assert_confidence_intervals(parci50,parci95,parfit,lb,ub) #=======================================================================