Release V1.1.3

.github/workflows/ci_PR.yml

@@ -12,7 +12,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ['3.10','3.11']
+        python-version: ['3.11','3.12']
 
     steps:
       - uses: actions/checkout@v3

.github/workflows/package_upload.yml

@@ -11,10 +11,10 @@ jobs:
     steps:
     - name: Checkout
       uses: actions/checkout@v1
-    - name: Set up Python 3.10
+    - name: Set up Python 3.11
       uses: actions/setup-python@v1
       with:
-        python-version: "3.10"
+        python-version: "3.11"
     - name: Install pypa/build
       run: >-
         python -m pip install build --user

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2019-2023 Luis Fabregas, Stefan Stoll, Gunnar Jeschke, and other contributors
+Copyright (c) 2019-2024 Luis Fabregas, Stefan Stoll, Gunnar Jeschke, and other contributors
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

VERSION

@@ -1 +1 @@
-v1.1.0
+v1.1.3

deerlab/classes.py

@@ -381,7 +381,7 @@ def percentile(self,p):
             cdf = np.cumsum(pdf)
             cdf /= max(cdf)
             # Eliminate duplicates
-            cdf, index = np.lib.arraysetops.unique(cdf,return_index=True)
+            cdf, index = np.unique(cdf,return_index=True)
             # Interpolate requested percentile
             x[n] = np.interp(p/100,cdf,values[index])
 

deerlab/dd_models.py

@@ -110,7 +110,7 @@ def _multigaussfun(r,r0,sig):
     P = np.sqrt(1/(2*np.pi))*1/sig*np.exp(-0.5*((r.T-r0)/sig)**2)
     if not np.all(P==0):
         # Normalization
-        P = np.squeeze(P)/np.sum([np.trapz(c,r) for c in P.T])
+        P = np.squeeze(P)/np.sum([np.trapezoid(c,r) for c in P.T])
     else: 
         P = np.squeeze(P)
     return P
@@ -125,7 +125,7 @@ def _multirice3dfun(r,nu,sig):
     r = r.T
     s2 = sig**2
     I_scaled = spc.ive(n/2-1, nu*r/s2)
-    P = nu**(n/2-1)/s2*r**(n/2)*np.exp(-(r**2+nu**2)/(2*s2)+nu*r/s2)*I_scaled
+    P = nu**(1-n/2)/s2*r**(n/2)*np.exp(-(r**2+nu**2)/(2*s2)+nu*r/s2)*I_scaled
     P[P<0] = 0
 
     # Normalization

deerlab/dipolarkernel.py

@@ -367,7 +367,7 @@ def dipolarkernel(t, r, *, pathways=None, mod=None, bg=None, method='fresnel', e
 
     if tinterp is not None:
         # Construct interpolator, this way elementarykernel_twospin is executed only once independently of how many pathways there are
-        Kinterpolator = [interp1d(tinterp,elementarykernel_twospin(tinterp,r_q,method,excbandwidth,gridsize,g,orisel,complex),axis=0,kind='cubic') for r_q in r]
+        Kinterpolator = _elementarykernel_twospin_interp(tinterp,r,method,excbandwidth,gridsize,g,orisel,complex)
         withinInterpolation = lambda tdip: np.all((np.max(tinterp) >= np.max(tdip)) & (np.min(tinterp) <= np.min(tdip)))
 
     # Define kernel matrix auxiliary functions
@@ -435,8 +435,14 @@ def K0_3spin(tdip):
 
     return K
 #==============================================================================
-
-
+@cached(max_size=100)
+def _elementarykernel_twospin_interp(tinterp,r,method,excbandwidth,gridsize,g,orisel,complex):
+    """
+    Construct interpolator, this way elementarykernel_twospin is executed only once independently of how many pathways there are
+    Cached for performance reasons, interpolation is slow.
+    """
+    Kinterpolator = [interp1d(tinterp,elementarykernel_twospin(tinterp,r_q,method,excbandwidth,gridsize,g,orisel,complex),axis=0,kind='cubic') for r_q in r]
+    return Kinterpolator
 
 #==============================================================================
 #   TWO-SPIN ELEMENTARY DIPOLAR KERNEL
@@ -468,14 +474,15 @@ def elementarykernel_twospin(tdip,r,method,ωex,gridsize,g,Pθ,complex):
     ωr = (μ0/2)*μB**2*g[0]*g[1]/h*1e21/(r**3)  # rad μs^-1
 
     # Orientation selection
-    orientationselection = Pθ is not None 
-    if orientationselection:
-        # Ensure zero-derivatives at [0,π/2]
-        θ = np.linspace(0,π/2,50) # rad
-        Pθ_ = make_interp_spline(θ, Pθ(θ),bc_type="clamped")
-        # Ensure normalization of probability density function (∫P(cosθ)dcosθ=1)
-        Pθnorm,_ = quad(lambda cosθ: Pθ_(np.arccos(cosθ)),0,1,limit=1000)
-        Pθ = lambda θ: Pθ_(θ)/Pθnorm
+    if method != 'fresnel':
+        orientationselection = Pθ is not None 
+        if orientationselection:
+            # Ensure zero-derivatives at [0,π/2]
+            θ = np.linspace(0,π/2,50) # rad
+            Pθ_ = make_interp_spline(θ, Pθ(θ),bc_type="clamped")
+            # Ensure normalization of probability density function (∫P(cosθ)dcosθ=1)
+            Pθnorm,_ = quad(lambda cosθ: Pθ_(np.arccos(cosθ)),0,1,limit=1000)
+            Pθ = lambda θ: Pθ_(θ)/Pθnorm
 
     def elementarykernel_fresnel(tdip):
     #------------------------------------------------------------------------

deerlab/diststats.py

@@ -7,7 +7,7 @@
 import warnings
 import copy
 from scipy.signal import find_peaks
-from scipy.integrate import cumtrapz
+from scipy.integrate import cumulative_trapezoid as cumtrapz
 
 def diststats(r, P, Puq=None, verbose=False, threshold=None):
     r""" 
@@ -109,7 +109,7 @@ def normalize(P):
     # Percentile function
     def pctile(r,P,p):
         cdf = cumtrapz(normalize(P),r,initial=0)
-        cdf, index = np.lib.arraysetops.unique(cdf,return_index=True)
+        cdf, index = np.unique(cdf,return_index=True)
         rpctile = np.interp(p/100,cdf,r[index])
         return rpctile
     # Expectation operator function

deerlab/fitresult.py

@@ -42,6 +42,8 @@ class FitResult(dict):
         * ``stats['aic']`` - Akaike information criterion
         * ``stats['aicc']`` - Corrected Akaike information criterion
         * ``stats['bic']`` - Bayesian information criterion
+        * ``stats['autocorr']`` - Autocorrelation based on Durbin–Watson statistic
+
 
     nonlin : ndarray
         Fitted non-linear parameters. [:ref:`snlls` specific attribute]

deerlab/selregparam.py

@@ -8,8 +8,8 @@
 import math as m
 import deerlab as dl
 
-def selregparam(y, A, solver, method='aic', algorithm='brent', noiselvl=None,
-                searchrange=[1e-8,1e2],regop=None, weights=None, full_output=False, candidates=None):
+def selregparam(y, A, solver, method='aic', algorithm='auto', noiselvl=None,
+                searchrange=[1e-8,1e2],regop=None, weights=None, full_output=False):
     r"""
     Selection of optimal regularization parameter based on a selection criterion.
 
@@ -52,6 +52,8 @@ def selregparam(y, A, solver, method='aic', algorithm='brent', noiselvl=None,
         * ``'ncp'`` - Normalized Cumulative Periodogram (NCP)
         * ``'gml'`` - Generalized Maximum Likelihood (GML)
         * ``'mcl'`` - Mallows' C_L (MCL)
+
+        If  ``'lr'`` or ``'lc'`` is specified, the search algorithm is automatically set to ``'grid'``.
 
     weights : array_like, optional
         Array of weighting coefficients for the individual datasets in global fitting.
@@ -60,18 +62,16 @@ def selregparam(y, A, solver, method='aic', algorithm='brent', noiselvl=None,
     algorithm : string, optional
         Search algorithm: 
 
+        * ``'auto'`` - Automatically, selects algrothium based on the searchrange. If the searchrange has two elements its set to ``'brent'`` otherwise to ``'grid'``.
         * ``'grid'`` - Grid-search, slow.
         * ``'brent'`` - Brent-algorithm, fast.
 
-        The default is ``'brent'``.
-
-    searchrange : two-element list, optional 
-        Search range for the optimization of the regularization parameter with the ``'brent'`` algorithm.
-        If not specified the default search range defaults to ``[1e-8,1e2]``.
+        The default is ``'auto'``.
 
-    candidates : list, optional
-        List or array of candidate regularization parameter values to be evaluated with the ``'grid'`` algorithm. 
-        If not specified, these are automatically computed from a grid within ``searchrange``.
+    searchrange : list, optional 
+        Either the search range for the optimization of the regularization parameter with the ``'brent'`` algorithm.
+        Or if more than two values are specified, then it is interpreted as candidates for the ``'grid'`` algorithm.
+        If not specified the default search range defaults to ``[1e-8,1e2]`` and the ``'brent'`` algorithm.
 
     regop : 2D array_like, optional
         Regularization operator matrix, the default is the second-order differential operator.
@@ -108,6 +108,13 @@ def selregparam(y, A, solver, method='aic', algorithm='brent', noiselvl=None,
     # If multiple datasets are passed, concatenate the datasets and kernels
     y, A, weights,_,__, noiselvl = dl.utils.parse_multidatasets(y, A, weights, noiselvl)
 
+    if algorithm == 'auto' and len(searchrange) == 2:
+        algorithm = 'brent'
+    elif algorithm == 'auto' and len(searchrange) > 2:
+        algorithm = 'grid'
+    elif algorithm  == 'auto' and len(searchrange) < 2:
+        raise ValueError("`searchrange` must have at least two elements if `algorithm` is set to `'auto'")
+
     # The L-curve criteria require a grid-evaluation
     if method == 'lr' or method == 'lc':
         algorithm = 'grid'
@@ -121,9 +128,11 @@ def selregparam(y, A, solver, method='aic', algorithm='brent', noiselvl=None,
     evalalpha = lambda alpha: _evalalpha(alpha, y, A, L, solver, method, noiselvl, weights)
 
     # Evaluate functional over search range, using specified search method
-    if algorithm == 'brent':
-                    
+    if algorithm == 'brent':    
+
         # Search boundaries
+        if len(searchrange) != 2:
+            raise ValueError("Search range must have two elements for the 'brent' algorithm.")
         lga_min = m.log10(searchrange[0])
         lga_max = m.log10(searchrange[1])
 
@@ -148,10 +157,10 @@ def register_ouputs(optout):
     elif algorithm=='grid':
 
         # Get range of potential alpha values candidates
-        if candidates is None:
+        if len(searchrange) == 2:
             alphaCandidates = 10**np.linspace(np.log10(searchrange[0]),np.log10(searchrange[1]),60)
         else: 
-            alphaCandidates = np.atleast_1d(candidates)
+            alphaCandidates = np.atleast_1d(searchrange)
 
         # Evaluate the full grid of alpha-candidates 
         functional,residuals,penalties,alphas_evaled = tuple(zip(*[evalalpha(alpha) for alpha in alphaCandidates]))
@@ -176,6 +185,10 @@ def register_ouputs(optout):
 
         # Find minimum of the selection functional              
         alphaOpt = alphaCandidates[np.argmin(functional)]
+        functional = np.array(functional)
+        residuals = np.array(residuals)
+        penalties = np.array(penalties)
+        alphas_evaled = np.array(alphas_evaled)
     else: 
         raise KeyError("Search method not found. Must be either 'brent' or 'grid'.")
 

deerlab/solvers.py

@@ -426,6 +426,11 @@ def snlls(y, Amodel, par0=None, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver
             * ``0`` : Work silently (default).
             * ``1`` : Display progress including the non-linear least-squares' solver termination report.
             * ``2`` : Display progress including the non-linear least-squares' solver iteration details.
+
+        .. caution::
+
+            The verbose output from the non-linear least-squares solver uses a different definition of the cost function than DeerLab.
+            DeerLab uses the sum of squares of the residuals divided by the number of data points, whereas the non-linear least-squares solver uses the sum of squares of the residuals divided by 2.
 
     Returns
     -------
@@ -457,6 +462,7 @@ def snlls(y, Amodel, par0=None, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver
         * ``stats['aic']`` - Akaike information criterion
         * ``stats['aicc']`` - Corrected Akaike information criterion
         * ``stats['bic']`` - Bayesian information criterion
+        * ``stats['autocorr']`` - Autocorrelation based on Durbin–Watson statistic
     success : bool
         Whether or not the optimizer exited successfully.
     cost : float
@@ -594,9 +600,10 @@ def linear_problem(y,A,optimize_alpha,alpha):
 
 
         if optimize_alpha:
-            output = dl.selregparam((y-yfrozen)[mask], Ared[mask,:], linSolver, regparam, 
-                                        weights=weights[mask], regop=L, candidates=regparamrange, 
-                                        noiselvl=noiselvl,searchrange=regparamrange,full_output=True)
+            linsolver_result = lambda AtA, Aty: parseResult(linSolver(AtA, Aty))
+            output = dl.selregparam((y-yfrozen)[mask], Ared[mask,:], linsolver_result, regparam, 
+                                        weights=weights[mask], regop=L, noiselvl=noiselvl,
+                                        searchrange=regparamrange,full_output=True)
             alpha = output[0]
             alpha_stats['alphas_evaled'] = output[1]
             alpha_stats['functional'] = output[2]

deerlab/utils.py

@@ -282,6 +282,7 @@ def goodness_of_fit(x,xfit,Ndof,noiselvl):
             stats['aic'] - Akaike information criterion
             stats['aicc'] - Corrected Akaike information criterion
             stats['bic'] - Bayesian information criterion
+            stats['autocorr'] - Autocorrelation based on Durbin–Watson statistic
     """
     sigma = noiselvl
     Ndof = np.maximum(Ndof,1)
@@ -833,15 +834,15 @@ def sophegrid(octants,maxphi,size):
         weights = np.zeros(nOrientations)
 
         sindth2 = np.sin(dtheta/2)
-        w1 = 0.5
+        w1 = 1.0
 
         # North pole (z orientation)
         phi[0] = 0
         theta[0] = 0
         weights[0] = maxphi*(1 - np.cos(dtheta/2))
 
         # All but equatorial slice
-        Start = 2
+        Start = 1
         for iSlice in np.arange(2,size):
             nPhi = nOct*(iSlice-1) + 1
             dPhi = maxphi/(nPhi - 1)
@@ -854,13 +855,13 @@ def sophegrid(octants,maxphi,size):
         # Equatorial slice
         nPhi = nOct*(size - 1) + 1
         dPhi = maxphi/(nPhi - 1)
-        idx = Start + (np.arange(0,nPhi) - 1)
+        idx = Start + np.arange(0,nPhi)
         phi[idx] = np.linspace(0,maxphi,nPhi)
         theta[idx] = np.pi/2
         weights[idx] = sindth2*dPhi*np.concatenate([[w1], np.ones(nPhi-2), [0.5]])
 
         # Border removal
-        rmv = np.cumsum(nOct*np.arange(1,size-1)+1)+1 
+        rmv = np.cumsum(nOct*np.arange(1,size)+1)
         phi = np.delete(phi,rmv)
         theta = np.delete(theta,rmv)
         weights = np.delete(weights,rmv)

docsrc/source/changelog.rst

@@ -24,6 +24,27 @@ Release Notes
 - |fix| : Something which was not working as expected or leading to errors has been fixed.
 - |api| : This will require changes in your scripts or code.
 
+Release ``v1.1.3`` - July 2024
+------------------------------------------
+- |fix| : Removes unnecessary files from the docs
+- |efficiency| : Improves the performance of the ``dipolarkernel`` function by 10-30% (:pr:`473`), by caching the interpolation of he effective dipolar evolution time vector.
+- |api| : Removes the `candidates` input from `selregparam` and integrates its function into `regparamrange`.
+- |fix| : Adds support for Numpy 2.0
+- |fix| : Add support for Python 3.12
+
+
+Release ``v1.1.2`` - November 2023
+------------------------------------------
+- |fix| : Fixes an issue with sophgrid (:pr:`463`).
+- |fix| : Fixes an error in the normalisation of the rice models (:pr:`459`).
+- |fix| : Removes the broken three spin example (:pr:`427`).
+- |fix| : Fixes an error in the linear solver for linearly constrained not non-negative problems.
+
+Release ``v1.1.1`` - August 2023
+------------------------------------------
+- |fix| : Fixes an error in the `FitResult.evaluate` function. (:pr:`454`)
+
+
 Release ``v1.1.0`` - August 2023
 ------------------------------------------
 - |api| : The definition of the dipolar time axis for RIDME has changed to match the one for 4-pulse DEER (:pr:`436`).