From b6584be197f3145ceb5429682bb9c232de44b7b3 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Thu, 18 Mar 2021 20:55:36 +0100 Subject: [PATCH 01/12] snlls: remove lines leading to double smoothness penalization --- deerlab/snlls.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/deerlab/snlls.py b/deerlab/snlls.py index 95417a25e..afeea3fec 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -332,9 +332,7 @@ def ResidualsFcn(p): # 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) return res From 03700cb80f6113737be0106eced0b4a76f1e4b1c Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Fri, 19 Mar 2021 18:51:49 +0100 Subject: [PATCH 02/12] 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 --- deerlab/fitmodel.py | 45 +++++++------ deerlab/snlls.py | 158 ++++++++++++++++++++++++++------------------ test/test_snlls.py | 24 +++---- 3 files changed, 131 insertions(+), 96 deletions(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index ca0fa0d3b..b6c9de83a 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -371,13 +371,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 # ------------------------------- @@ -421,8 +424,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 # ----------------------- @@ -459,10 +461,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]) @@ -471,14 +470,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]) @@ -588,7 +580,8 @@ 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)] # Separable non-linear least squares (SNNLS) @@ -596,9 +589,21 @@ def separable_nonlinear_lsq_analysis(Vexp): regparam=regparam, uqanalysis=uqanalysis, weights=weights) 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 = [prescales[i]*np.trapz(Pfit,r) for i in range(nSignals)] + #scales = [prescales[i]*np.trapz(Pfit,r) for i in range(nSignals)] + scales = fit.scale + # Normalize distribution + # ----------------------- + Pscale = np.trapz(Pfit,r) + Pfit /= Pscale + scales = [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/snlls.py b/deerlab/snlls.py index afeea3fec..0f78ea424 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, custom_penalty=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,7 +193,7 @@ 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) @@ -201,7 +201,8 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx 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': @@ -247,9 +248,9 @@ 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))) + # Use an arbitrary axis + ax = np.arange(1, Nlin+1) if includePenalty: - # Use an arbitrary axis - ax = np.arange(1, Nlin+1) # Get regularization operator regorder = np.minimum(Nlin-1, regorder) L = dl.regoperator(ax, regorder) @@ -284,6 +285,26 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx regparam_prev = 0 par_prev = [0]*len(par0) + + def linear_problem(A,optimize_alpha): + #=========================================================================== + + global alpha + + # Regularization components + 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 + #=========================================================================== + def ResidualsFcn(p): #=========================================================================== """ @@ -293,44 +314,61 @@ def ResidualsFcn(p): non-linear least-squares solver. """ - nonlocal par_prev, check, regparam_prev, linfit, alpha + nonlocal par_prev, check, regparam_prev, scales, linfit + global alpha # Non-linear model evaluation A = Amodel(p) + # Regularization components if includePenalty: 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 = linear_problem(A,optimize_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 + + scales = [] + 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 = np.zeros_like(yfit) + for subset,scale in zip(subsets,scales): + scales_vec[subset] = scale + # Compute residual vector - res = weights*(yfit - y) + res = weights*(scales_vec*(Amodel(p)@linfit) - y) + + # Compute residual from custom penalty + if callable(custom_penalty): + penres = custom_penalty(p) + penres = np.atleast_1d(penres) + res = np.concatenate((res,penres)) + if includePenalty: # Augmented residual res, _ = _augment(res, [], regtype, alpha, L, linfit, huberparam, Nnonlin) @@ -361,54 +399,28 @@ def ResidualsFcn(p): 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) - - Jlin = Afit - J = np.concatenate((Jnonlin, Jlin),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) + # Compute the residual and Jacobian + res = ResidualsFcn(nonlinfit) + J = Jacobian(ResidualsFcn,nonlinfit,lb,ub) # 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_nonlin = UncertQuant('covariance', nonlinfit, covmatrix, lb, ub) + paramuq_lin = paramuq_nonlin.propagate(lambda p: linear_problem(Amodel(p),optimize_alpha=True),lbl,ubl) else: - paramuq = [] + paramuq_nonlin = [] + paramuq_lin = [] # Goodness-of-fit # -------------------------------------- @@ -420,15 +432,33 @@ 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 res_lambda(p,y,linfit,Amodel,weights,L,scales_vec=1): + res = weights*(scales_vec*(Amodel(p)@linfit) - y) + + lams = p[[1,2,3]] + if np.sum(lams)<1: + pen = 0 + else: + pen = np.atleast_1d(max(y)*np.linalg.norm(np.sum(lams) - 1))**2 + pen = np.full_like(res,pen) + + res = np.concatenate((res,pen)) + res = np.concatenate((res,L@P)) + + return np.concatenate((res,pen)) + def _augment(res, J, regtype, alpha, L, x, eta, Nnonlin): # =========================================================================================== 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) #======================================================================= From 9c355d40444cf6e7eb8c002aa3d5de420c6d5654 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Fri, 19 Mar 2021 19:03:30 +0100 Subject: [PATCH 03/12] fitmodel: implement constraint to ensure identifiability of Lam0 and V0 --- deerlab/fitmodel.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index b6c9de83a..f3a498307 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -573,6 +573,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" @@ -584,9 +585,21 @@ def separable_nonlinear_lsq_analysis(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: + penalty[i] = np.atleast_1d(0) + else: + penalty[i] = np.atleast_1d(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,custom_penalty=scale_constraint) parfit = fit.nonlin Pfit = fit.lin param_uq = fit.nonlinUncert From 42d536e9dab62478fb6c001bc4c9702ff3528e4c Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Sat, 20 Mar 2021 21:13:25 +0100 Subject: [PATCH 04/12] snlls: correct behaviour of the uncertainty quantification --- deerlab/fitmodel.py | 14 +++++--- deerlab/snlls.py | 83 ++++++++++++++++++++++++++++----------------- 2 files changed, 61 insertions(+), 36 deletions(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index f3a498307..558200080 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -586,17 +586,21 @@ def separable_nonlinear_lsq_analysis(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: - penalty[i] = np.atleast_1d(0) + if np.sum(lams)<1 or ex_model[i].__name__=='ex_4pdeer': + penalty[i] = 0 else: - penalty[i] = np.atleast_1d(max(Vexp[i])*(np.sum(lams) - 1)) + 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,custom_penalty=scale_constraint) @@ -606,8 +610,8 @@ def scale_constraint(nonlinpar): Pfit_uq = fit.linUncert snlls_uq = [param_uq,Pfit_uq] alphaopt = fit.regparam - #scales = [prescales[i]*np.trapz(Pfit,r) for i in range(nSignals)] scales = fit.scale + # Normalize distribution # ----------------------- Pscale = np.trapz(Pfit,r) diff --git a/deerlab/snlls.py b/deerlab/snlls.py index 0f78ea424..7c5121d49 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -200,7 +200,6 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx 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 @@ -288,7 +287,12 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx def linear_problem(A,optimize_alpha): #=========================================================================== - + """ + Linear problem + ------------------ + Solves the linear subproblem of the SNLLS objective function via linear LSQ + constrained, unconstrained, or regularized. + """ global alpha # Regularization components @@ -352,12 +356,11 @@ def ResidualsFcn(p): yfit = A@linfit 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 = np.zeros_like(yfit) - for subset,scale in zip(subsets,scales): scales_vec[subset] = scale # Compute residual vector @@ -371,12 +374,12 @@ def ResidualsFcn(p): if includePenalty: # Augmented residual - res, _ = _augment(res, [], regtype, alpha, L, linfit, huberparam, Nnonlin) + 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.') @@ -406,24 +409,45 @@ def ResidualsFcn(p): yfit = scales_vec*(Afit@linfit) # Uncertainty analysis - #-------------------------------------------------------- + #--------------------- if uqanalysis: - # Compute the residual and Jacobian + # Compute the fit residual res = ResidualsFcn(nonlinfit) - J = Jacobian(ResidualsFcn,nonlinfit,lb,ub) + + # Jacobian (non-linear part) + Jnonlin = Jacobian(ResidualsFcn,nonlinfit,lb,ub) + + # 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] + + # 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_nonlin = UncertQuant('covariance', nonlinfit, covmatrix, lb, ub) - paramuq_lin = paramuq_nonlin.propagate(lambda p: linear_problem(Amodel(p),optimize_alpha=True),lbl,ubl) + 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_nonlin = [] paramuq_lin = [] # Goodness-of-fit - # -------------------------------------- + # --------------- stats = [] for subset in subsets: Ndof = len(y[subset]) - Nnonlin @@ -444,23 +468,26 @@ def plotfcn(show=False): stats=stats, cost=fvals, residuals=sol.fun, success=sol.success, scale=scales) # =========================================================================================== -def res_lambda(p,y,linfit,Amodel,weights,L,scales_vec=1): - res = weights*(scales_vec*(Amodel(p)@linfit) - y) - lams = p[[1,2,3]] - if np.sum(lams)<1: - pen = 0 - else: - pen = np.atleast_1d(max(y)*np.linalg.norm(np.sum(lams) - 1))**2 - pen = np.full_like(res,pen) +def uq_subset(uq_full,subset): +#=========================================================================== + "Wrapper around the CI function handle of the uncertainty structure" + uq_subset = copy.deepcopy(uq_full) - res = np.concatenate((res,pen)) - res = np.concatenate((res,L@P)) + 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) - return np.concatenate((res,pen)) + # 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 _augment(res, J, regtype, alpha, L, x, eta, Nnonlin): +def reg_penalty(regtype, alpha, L, x, eta, Nnonlin): # =========================================================================================== """ LSQ residual and Jacobian augmentation @@ -486,13 +513,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 # =========================================================================================== From 21324d51ce1a690e0fe8f00aec2202b27b82b598 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Sat, 20 Mar 2021 21:16:56 +0100 Subject: [PATCH 05/12] tests: relax tight assertion for one test --- test/test_fitmodel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_fitmodel.py b/test/test_fitmodel.py index b300c5fe2..63ecb6bb0 100644 --- a/test/test_fitmodel.py +++ b/test/test_fitmodel.py @@ -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(): # ====================================================================== From 6938a9142bef82f2b5e192fcba5ee12310128f46 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Sat, 20 Mar 2021 22:03:19 +0100 Subject: [PATCH 06/12] hccm: fix minor bug in internal function --- deerlab/utils/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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))) From 7a5e220b33743c33485d40f991f47fced5efb53d Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Tue, 23 Mar 2021 09:59:13 +0100 Subject: [PATCH 07/12] fitmodel: remove deprecated keyword argument --- deerlab/fitmodel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index 558200080..1efd982a7 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'): 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), From 6899c8106a930d5eed1fabce4b308277550ac60c Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Tue, 23 Mar 2021 10:22:16 +0100 Subject: [PATCH 08/12] snlls: rename keyword custom_penalty to extrapenalty, fix some tests --- deerlab/fitmodel.py | 2 +- deerlab/snlls.py | 16 +++++++++------- test/test_fitmodel.py | 6 +++--- 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index 1efd982a7..3cbcd8718 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -603,7 +603,7 @@ def scale_constraint(nonlinpar): # 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,custom_penalty=scale_constraint) + regparam=regparam, uqanalysis=uqanalysis, weights=weights,extrapenalty=scale_constraint) parfit = fit.nonlin Pfit = fit.lin param_uq = fit.nonlinUncert diff --git a/deerlab/snlls.py b/deerlab/snlls.py index 7c5121d49..f7976311d 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, custom_penalty=None, + 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 @@ -295,12 +295,13 @@ def linear_problem(A,optimize_alpha): """ global alpha - # Regularization components + # 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) @@ -324,7 +325,7 @@ def ResidualsFcn(p): A = Amodel(p) - # Regularization components + # Check whether optimization of the regularization parameter is needed if includePenalty: if type(regparam) is str: # If the parameter vector has not changed by much... @@ -355,8 +356,8 @@ def ResidualsFcn(p): # Evaluate full model residual yfit = A@linfit - scales = [] - scales_vec = np.zeros_like(yfit) + # 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]) @@ -367,8 +368,8 @@ def ResidualsFcn(p): res = weights*(scales_vec*(Amodel(p)@linfit) - y) # Compute residual from custom penalty - if callable(custom_penalty): - penres = custom_penalty(p) + if callable(extrapenalty): + penres = extrapenalty(p) penres = np.atleast_1d(penres) res = np.concatenate((res,penres)) @@ -396,6 +397,7 @@ 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] diff --git a/test/test_fitmodel.py b/test/test_fitmodel.py index 63ecb6bb0..ba1de158c 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 From 735d9dd96af38d0cb00780444f619560dd7aa496 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Wed, 24 Mar 2021 08:22:02 +0100 Subject: [PATCH 09/12] snlls: code refactoring --- deerlab/lsqcomponents.py | 5 ++++- deerlab/snlls.py | 27 +++++++++++++-------------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/deerlab/lsqcomponents.py b/deerlab/lsqcomponents.py index de79e806c..0be8dc52d 100644 --- a/deerlab/lsqcomponents.py +++ b/deerlab/lsqcomponents.py @@ -1,11 +1,14 @@ 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 nr = np.shape(K)[1] + if L is None: + L = np.eye(nr) + # Compute residual terms KtV = weights*K.T@V KtK = weights*K.T@K diff --git a/deerlab/snlls.py b/deerlab/snlls.py index f7976311d..27352f957 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -205,9 +205,9 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx # 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 # ---------------------------------------------------------- @@ -249,12 +249,12 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx # Use an arbitrary axis ax = np.arange(1, Nlin+1) - if includePenalty: + 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,9 +283,9 @@ 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): + def linear_problem(A,optimize_alpha,alpha): #=========================================================================== """ Linear problem @@ -293,7 +293,6 @@ def linear_problem(A,optimize_alpha): Solves the linear subproblem of the SNLLS objective function via linear LSQ constrained, unconstrained, or regularized. """ - global alpha # Optimiza the regularization parameter only if needed if optimize_alpha: @@ -307,7 +306,7 @@ def linear_problem(A,optimize_alpha): linfit = parseResult(result) linfit = np.atleast_1d(linfit) - return linfit + return linfit, alpha #=========================================================================== def ResidualsFcn(p): @@ -319,14 +318,14 @@ def ResidualsFcn(p): non-linear least-squares solver. """ - nonlocal par_prev, check, regparam_prev, scales, linfit - global alpha + nonlocal par_prev, check, regparam_prev, scales, linfit, alpha + # Non-linear model evaluation A = Amodel(p) # Check whether optimization of the regularization parameter is needed - if includePenalty: + 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): @@ -350,7 +349,7 @@ def ResidualsFcn(p): optimize_alpha = False alpha = 0 - linfit = linear_problem(A,optimize_alpha) + linfit,alpha = linear_problem(A,optimize_alpha,alpha) regparam_prev = alpha # Evaluate full model residual @@ -373,11 +372,11 @@ def ResidualsFcn(p): penres = np.atleast_1d(penres) res = np.concatenate((res,penres)) - if includePenalty: + if includeRegularization : # Augmented residual res_reg, _ = reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin) res = np.concatenate((res,res_reg)) - + return res #=========================================================================== From 27a232b7af006124997e00464019aa8dfe6d1a96 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Wed, 24 Mar 2021 08:52:04 +0100 Subject: [PATCH 10/12] fitmodel: fix merge error --- deerlab/fitmodel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index 716969b4f..4e4225e34 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -609,7 +609,7 @@ def scale_constraint(nonlinpar): # 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,extrapenalty=scale_constraint) + regparam=regparam, uqanalysis=uqanalysis, weights=weights,extrapenalty=scale_constraint, nonlin_tol=tol,nonlin_maxiter=maxiter) parfit = fit.nonlin Pfit = fit.lin From ac83e6874cf137b8ba75797e191bd958d82c5f61 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Wed, 24 Mar 2021 08:59:55 +0100 Subject: [PATCH 11/12] fitmodel: fix another merge error --- deerlab/fitmodel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index 4e4225e34..883c3fd5d 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -627,7 +627,7 @@ def scale_constraint(nonlinpar): 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)/Pscal + Pfit_uq.ci = lambda p: Pfit_uq_.ci(p)/Pscale # Get the fitted models Kfit,Bfit = multiPathwayModel(parfit) From be11bd93a71aae62a2611af3ad088b9b4082c058 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Wed, 24 Mar 2021 19:49:08 +0100 Subject: [PATCH 12/12] lsqcomponents: refactor code --- deerlab/lsqcomponents.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/deerlab/lsqcomponents.py b/deerlab/lsqcomponents.py index 0be8dc52d..d2123509f 100644 --- a/deerlab/lsqcomponents.py +++ b/deerlab/lsqcomponents.py @@ -6,13 +6,13 @@ def lsqcomponents(V,K,L=None,alpha=0,weights=1, regtype = 'tikh',huberparam = 1. # Prepare nr = np.shape(K)[1] - if L is None: - L = np.eye(nr) - # Compute residual terms 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)