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/ci_scheduled.yml

@@ -12,7 +12,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: [3.8, 3.9, "3.10", "3.11", "3.12"]
+        python-version: [3.8, 3.9, "3.10", "3.11", "3.12","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

VERSION

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

deerlab/dd_models.py

@@ -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""" 

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/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,7 +600,8 @@ def linear_problem(y,A,optimize_alpha,alpha):
 
 
         if optimize_alpha:
-            output = dl.selregparam((y-yfrozen)[mask], Ared[mask,:], linSolver, regparam, 
+            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, candidates=regparamrange, 
                                         noiselvl=noiselvl,searchrange=regparamrange,full_output=True)
             alpha = output[0]

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,24 @@ 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`` - Ongoing
+------------------------------------------
+- |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.
+- |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`).

docsrc/source/installation.rst

@@ -45,7 +45,7 @@ DeerLab installs the following packages:
 * `joblib <https://joblib.readthedocs.io/en/latest/>`_ - Library lightweight pipelining and parallelization.
 * `tqdm <https://github.com/tqdm/tqdm>`_ - A lightweight package for smart progress meters.
 * `dill <https://github.com/uqfoundation/dill>`_ - An extension of Python's pickle module for serializing and de-serializing python objects.
-* `quadprog <https://pypi.org/project/quadprog/>`_ - A quadratic programming solver
+* `quadprog <https://pypi.org/project/quadprog/>`_ - A quadratic programming solver (Only for Python versions < 3.11)
 
 Importing DeerLab
 ------------------
@@ -74,9 +74,6 @@ Any DeerLab version released after v0.10.0 can be installed via ``pip`` using th
 
 		python -m pip install deerlab==x.y.z
 
-DeerLab version prior to 0.10 are written in MATLAB and are still available from an `archived repository <https://github.com/JeschkeLab/DeerLab-Matlab>`_. 
-Download and installation instruction for the MATLAB environment are provided in the released documentation. MATLAB releases have been deprecated and no further support is provided.
-
 
 Installing from source
 *****************************
@@ -130,10 +127,18 @@ On a **Windows** computer, if you are trying to run a DeerLab function, you migh
 
     ImportError: DLL load failed: The specified module could not be found.
 
-This happens when the MKL libraries have not been properly linked in ``numpy``, ``scipy`` or ``cvxopt`` 
-installations (typically ``numpy``). This can happen due to firewall restrictions,
-user rights, or internet connection issues during the DeerLab installation. To solve this, the
-best solution is to manually install as follows. 
+This issue can occur when a specific python package is not installled properly. This typicaly happens when installing the MKL linked libaries but can 
+occur when installing packages normaly from PyPI. This can happen due to firewall restrictions,
+user rights, or internet connection issues during the DeerLab installation.
+
+If your where installing packages from PyPI (i.e. via pip) then please uninstall the package and reinstall it.
+
+.. code-block:: text
+
+    python -m pip uninstall <package_name>
+    python -m pip install <package_name>
+
+If you were trying to install the MKL linked libraries then please follow the instructions below.
 
 1) Go to https://www.lfd.uci.edu/~gohlke/pythonlibs/#numpy
 
@@ -167,7 +172,7 @@ During installation on certain systems (e.g. some computation clusters) using on
 the following error might be raised during the installation:
 
 .. code-block:: text
-
+    
     Error while finding module specification for 'setup.py'
     (ModuleNotFoundError: __path__ attribute not found on 'setup' while trying to find 'setup.py')
 

test/test_utils.py

@@ -0,0 +1,16 @@
+import numpy as np
+from deerlab import sophegrid
+
+
+# ======================================================================
+# Test sophegrid function
+
+def test_sophegrid():
+    # Test that the function returns the values and weights summing to 1
+    # Comparison taken from EasySpin 6.0
+    phi, theta, weights = sophegrid(4,np.pi*2,3)
+
+    assert np.allclose(weights.sum(),1)
+    assert np.allclose(phi, np.array([0,0,1.57079632679490,3.14159265358979,4.71238898038469,0,0.785398163397448,1.57079632679490,2.35619449019235,3.14159265358979,3.92699081698724,4.71238898038469,5.49778714378214]))
+    assert np.allclose(theta, np.array([0,0.785398163397448,0.785398163397448,0.785398163397448,0.785398163397448,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490]))
+    assert np.allclose(weights*4*np.pi, np.array([0.956558005801449,1.70021769237074,1.70021769237074,1.70021769237074,1.70021769237074,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346]))