From 863eaa969af4a3f2e0a00f867374fb7a26e08b28 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Mon, 10 May 2021 09:48:50 +0200 Subject: [PATCH 1/2] snlls: allow `extrapenalty` to take linear parameters as input --- deerlab/fitmodel.py | 2 +- deerlab/snlls.py | 16 +++++++++------- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/deerlab/fitmodel.py b/deerlab/fitmodel.py index 9ebff6ac1..1425ceaf2 100644 --- a/deerlab/fitmodel.py +++ b/deerlab/fitmodel.py @@ -568,7 +568,7 @@ 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): + def scale_constraint(nonlinpar,linpar): # -------------------------------------------------------- penalty = np.zeros(nSignals) for i in range(nSignals): diff --git a/deerlab/snlls.py b/deerlab/snlls.py index 3628a1b2d..5656f9bc1 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -87,9 +87,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. + extrapenalty: callable + Custom penalty function to impose upon the solution. Must take two inputs, a vector of non-linear parameters + and a vector of linear parameters, and return a vector to be added to the residual vector (``pen = fcn(pnonlin,plin)``). + The square of the penalty is computed internally. alphareopt : float scalar, optional Relative parameter change threshold for reoptimizing the regularization parameter @@ -383,7 +384,7 @@ def ResidualsFcn(p): # Compute residual from custom penalty if callable(extrapenalty): - penres = extrapenalty(p) + penres = extrapenalty(p,linfit) penres = np.atleast_1d(penres) res = np.concatenate((res,penres)) @@ -434,10 +435,11 @@ def ResidualsFcn(p): Jnonlin = Jacobian(ResidualsFcn,nonlinfit,lb,ub) # Jacobian (linear part) - Jlin = np.zeros((len(res),len(linfit))) - Jlin[:len(y),:] = scales_vec[:,np.newaxis]*Amodel(nonlinfit) + Jlin = scales_vec[:,np.newaxis]*Amodel(nonlinfit) if includeRegularization: - Jlin[len(res)-Nlin:,:] = reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin)[1] + Jlin = np.concatenate((Jlin, reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin)[1])) + if callable(extrapenalty): + Jlin = np.concatenate((Jlin, Jacobian(lambda plin:extrapenalty(nonlinfit,plin),linfit,lbl,ubl))) # Full Jacobian J = np.concatenate((Jnonlin,Jlin),axis=1) From 1335808af5b67bd3ee92a7cdbd6e4c520a6af4f8 Mon Sep 17 00:00:00 2001 From: Luis Fabregas Date: Mon, 10 May 2021 12:54:44 +0200 Subject: [PATCH 2/2] snlls: fix order of concatenation of Jacobians --- deerlab/snlls.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/deerlab/snlls.py b/deerlab/snlls.py index 5656f9bc1..ab75dcca3 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -436,10 +436,10 @@ def ResidualsFcn(p): # Jacobian (linear part) Jlin = scales_vec[:,np.newaxis]*Amodel(nonlinfit) + if callable(extrapenalty): + Jlin = np.concatenate((Jlin, Jacobian(lambda plin: extrapenalty(nonlinfit,plin),linfit,lbl,ubl))) if includeRegularization: Jlin = np.concatenate((Jlin, reg_penalty(regtype, alpha, L, linfit, huberparam, Nnonlin)[1])) - if callable(extrapenalty): - Jlin = np.concatenate((Jlin, Jacobian(lambda plin:extrapenalty(nonlinfit,plin),linfit,lbl,ubl))) # Full Jacobian J = np.concatenate((Jnonlin,Jlin),axis=1)