Robust Jacobian implementation (#58)

* implemented robust Jacobian from scipy that accounts for boundaries

deerlab/classes.py

@@ -2,7 +2,7 @@
 # Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors.
 
 import numpy as np
-from deerlab.utils import fdJacobian
+from deerlab.utils import Jacobian
 from scipy.stats import norm
 from scipy.signal import fftconvolve
 import copy
@@ -343,7 +343,7 @@ def propagate(self,model,lbm=[],ubm=[]):
             raise IndexError ('The 2nd and 3rd input arguments must have the same number of elements as the model output.')
 
         # Get jacobian of model to be propagated with respect to parameters
-        J = fdJacobian(model,parfit)
+        J = Jacobian(model,parfit,self.__lb,self.__ub)
 
         # Clip at boundaries
         modelfit = np.maximum(modelfit,lbm)

deerlab/fitmultimodel.py

@@ -8,7 +8,7 @@
 import matplotlib.pyplot as plt
 from types import FunctionType
 import deerlab as dl
-from deerlab.utils import hccm, goodness_of_fit, fdJacobian
+from deerlab.utils import hccm, goodness_of_fit, Jacobian
 from deerlab.classes import FitResult
 
 def fitmultimodel(V, Kmodel, r, model, maxModels, method='aic', lb=None, ub=None, lbK=None, ubK=None,
@@ -342,7 +342,7 @@ def spread_within_box(lb,ub,Nmodels):
         nlin_ub_.append(nlin_ub)
         nlin_lb_.append(nlin_lb)   
         lin_ub_.append(lin_ub)
-        lin_lb_.append(lin_lb)
+        lin_lb_.append(lin_lb-1) # rescale to zero
         fits.append(fit)
     # Select the optimal model
     # ========================
@@ -366,16 +366,19 @@ def spread_within_box(lb,ub,Nmodels):
         Knonlin = lambda par: nonlinmodel(par,Nopt)
         res = weights*(Vfit - V)
 
+        lb_full = np.concatenate((nlin_lb, lin_lb))
+        ub_full = np.concatenate((nlin_ub, lin_ub))
+
         # Compute the Jacobian
-        Jnonlin = fdJacobian(lambda p: weights*(Knonlin(p)@plin),pnonlin)
+        Jnonlin = Jacobian(lambda p: weights*(Knonlin(p)@plin),pnonlin,nlin_lb,nlin_ub)
         Jlin = weights[:,np.newaxis]*Knonlin(pnonlin)
         J = np.concatenate((Jnonlin, Jlin),1)
 
         # Estimate the heteroscedasticity-consistent covariance matrix
         covmatrix = hccm(J,res,'HC1')
 
         # Construct uncertainty quantification structure for fitted parameters
-        paramuq = dl.UncertQuant('covariance',np.concatenate((pnonlin, plin)),covmatrix,np.concatenate((nlin_lb, lin_lb)),np.concatenate((nlin_ub, lin_ub)))
+        paramuq = dl.UncertQuant('covariance',np.concatenate((pnonlin, plin)),covmatrix,lb_full,ub_full)
 
         P_subset = np.arange(0, nKparam+nparam*Nopt)
         amps_subset = np.arange(P_subset[-1]+1, P_subset[-1]+1+Nopt)

deerlab/fitparamodel.py

@@ -4,7 +4,7 @@
 # Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors.
 
 import numpy as np
-from deerlab.utils import isempty, multistarts, hccm, parse_multidatasets, goodness_of_fit, fdJacobian
+from deerlab.utils import isempty, multistarts, hccm, parse_multidatasets, goodness_of_fit, Jacobian
 from deerlab.classes import UncertQuant, FitResult
 import matplotlib.pyplot as plt
 from scipy.optimize import least_squares
@@ -215,7 +215,7 @@ def lsqresiduals(p):
         # of the negative log-likelihood
         with warnings.catch_warnings():
             warnings.simplefilter('ignore')
-            J = fdJacobian(lsqresiduals,parfit)
+            J = Jacobian(lsqresiduals,parfit,lb,ub)
 
         # Estimate the heteroscedasticity-consistent covariance matrix
         if isempty(covmatrix):

deerlab/fitsignal.py

@@ -12,7 +12,7 @@
 from deerlab.classes import UncertQuant, FitResult
 from deerlab.bg_models import bg_hom3d
 from deerlab.ex_models import ex_4pdeer
-from deerlab.utils import isempty, goodness_of_fit, fdJacobian
+from deerlab.utils import isempty, goodness_of_fit, Jacobian
 
 def fitsignal(Vexp, t, r, dd_model='P', bg_model=bg_hom3d, ex_model=ex_4pdeer,
               par0=[None,None,None], lb=[None,None,None], ub=[None,None,None], verbose= False,
@@ -403,7 +403,7 @@ def splituq(full_uq,scales=1):
                 Vfit_uq.append( UncertQuant('covariance',Vfit[jj],Vcovmat,[],[]))
             elif includeForeground:
                 # Parametric signal with parameter-free distribution
-                Jnonlin = fdJacobian(lambda par: multiPathwayModel(par[paramidx])[0][jj]@Pfit,parfit_)
+                Jnonlin = Jacobian(lambda par: multiPathwayModel(par[paramidx])[0][jj]@Pfit,parfit_,lb[paramidx],ub[paramidx])
                 J = np.concatenate((Jnonlin, Kfit[jj]),1)
                 Vcovmat = J@covmat@J.T
                 Vfit_uq.append( UncertQuant('covariance',Vfit[jj],Vcovmat,[],[]))

deerlab/snlls.py

@@ -11,7 +11,7 @@
 
 # Import DeerLab depencies
 import deerlab as dl
-from deerlab.utils import goodness_of_fit, hccm, isempty, fdJacobian
+from deerlab.utils import goodness_of_fit, hccm, isempty, Jacobian
 from deerlab.nnls import cvxnnls, fnnls, nnlsbpp
 from deerlab.classes import UncertQuant, FitResult
 
@@ -364,7 +364,7 @@ def ResidualsFcn(p):
 
         # Compute the Jacobian for the linear and non-linear parameters
         fcn = lambda p: Amodel(p)@linfit
-        Jnonlin = fdJacobian(fcn,nonlinfit)
+        Jnonlin = Jacobian(fcn,nonlinfit,lb,ub)
 
         Jlin = Afit
         J = np.concatenate((Jnonlin, Jlin),1)

deerlab/utils/utils.py

@@ -1,18 +1,12 @@
-from builtins import str
-
-def isstring(x):
-    bool = isinstance(x,str)
-    return bool
-
-
 import warnings
-import matplotlib.pyplot as plt
 import numpy as np
 import cmath as math
 import scipy as scp
-import random
-from scipy.sparse import coo_matrix
+import scipy.optimize as opt
+
 from types import FunctionType 
+
+
 def parse_multidatasets(V,K,weights,precondition=False):
 #===============================================================================
 
@@ -410,25 +404,18 @@ def diagp(Y,X,k):
 #===============================================================================
 
 #===============================================================================
-def fdJacobian(fcn, param):
+def Jacobian(fcn, x0, lb, ub):
     """ 
-    Finite difference Jacobian
-    Estimates the Jacobian matrix of a function ``fcn`` at the point defined by ``param``. 
+    Finite difference Jacobian estimation 
+    Estimates the Jacobian matrix of a vector-valued function ``fcn`` at the 
+    point ``x0`` taking into consideration box-constraints defined by the lower
+    and upper bounds ``lb`` and ``ub``.
+
+    This is a wrapper around the ``scipy.optimize._numdiff.approx_derivative`` function.
 
     """
-    param = np.atleast_1d(param)
-    # Step size
-    h = 1e-5
-    # Central value
-    f0 = fcn(param)
-    f0 = np.atleast_1d(f0)
-    # Loop over parameters
-    J = np.zeros((len(f0),len(param)))
-    for i in range(len(param)):
-        rise = np.copy(param)
-        rise[i] = rise[i] + h
-        # Finite difference derivative of i-th parameter
-        J[:,i] = (fcn(rise) - f0)/h
+    J = opt._numdiff.approx_derivative(fcn,x0,method='2-point',bounds=(lb,ub))
+    J = np.atleast_2d(J)
     return J
 #===============================================================================