From 2546b4290fe4530575462fb343182e29bba069a3 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Fri, 3 Nov 2023 12:20:13 +0100 Subject: [PATCH 01/10] Push patch into V1.1 (#467) * Increase version number and update changelog (#455) * Bugfixes 4th sep (#460) * Fix normalisation in rice model Closes #459 * Remove three spin anaysis Closes #427 * Improved installation instructions * Caution about difference in definition of cost function Closes #450 * Bump Version * Bug for non linearly constrained problems When a problem is not linearly constrained and not non-negative, the linear solver outputs a result class not the solution. * Keeping changelog up-to date * Fixing Sophgrid bug (#464) * Fixing Sophgrid bug * Add unit test for sophgrid * Update changelog.rst for new release (#466) --- VERSION | 2 +- deerlab/dd_models.py | 2 +- deerlab/solvers.py | 8 +- deerlab/utils.py | 8 +- docsrc/source/changelog.rst | 12 +++ docsrc/source/installation.rst | 23 +++-- .../advanced/ex_long_threespin_analysis.py | 95 ------------------- test/test_utils.py | 16 ++++ 8 files changed, 55 insertions(+), 111 deletions(-) delete mode 100644 examples/advanced/ex_long_threespin_analysis.py create mode 100644 test/test_utils.py diff --git a/VERSION b/VERSION index 795460fc..0f1acbd5 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v1.1.0 +v1.1.2 diff --git a/deerlab/dd_models.py b/deerlab/dd_models.py index afcd381d..dc32cbaa 100644 --- a/deerlab/dd_models.py +++ b/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 diff --git a/deerlab/solvers.py b/deerlab/solvers.py index 8c1c6a29..1052d7b3 100644 --- a/deerlab/solvers.py +++ b/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 ------- @@ -594,7 +599,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] diff --git a/deerlab/utils.py b/deerlab/utils.py index 5079bd5b..3b7f69d1 100644 --- a/deerlab/utils.py +++ b/deerlab/utils.py @@ -833,7 +833,7 @@ 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 @@ -841,7 +841,7 @@ def sophegrid(octants,maxphi,size): 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 +854,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) diff --git a/docsrc/source/changelog.rst b/docsrc/source/changelog.rst index 9da8425e..894c269f 100644 --- a/docsrc/source/changelog.rst +++ b/docsrc/source/changelog.rst @@ -24,6 +24,18 @@ 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.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`). diff --git a/docsrc/source/installation.rst b/docsrc/source/installation.rst index dd974a46..95257bc3 100644 --- a/docsrc/source/installation.rst +++ b/docsrc/source/installation.rst @@ -45,7 +45,7 @@ DeerLab installs the following packages: * `joblib `_ - Library lightweight pipelining and parallelization. * `tqdm `_ - A lightweight package for smart progress meters. * `dill `_ - An extension of Python's pickle module for serializing and de-serializing python objects. -* `quadprog `_ - A quadratic programming solver +* `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 `_. -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 + python -m pip install + +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') diff --git a/examples/advanced/ex_long_threespin_analysis.py b/examples/advanced/ex_long_threespin_analysis.py deleted file mode 100644 index cdbde826..00000000 --- a/examples/advanced/ex_long_threespin_analysis.py +++ /dev/null @@ -1,95 +0,0 @@ -#%% -""" -Analyzing 4-pulse DEER data acquired on three-spin systems -============================================================================ - -As in the publication referenced below, this example will take two 4-pulse DEER signals acquired -on the same protein sample conatining three nitroxide spins with different attenuation levels of the pump -pulse power. - -For the original model and more information on these systems please refer to: -L. Fábregas Ibáñez, M. H. Tessmer, G. Jeschke, and S. Stoll. -Dipolar pathways in multi-spin and multi-dimensional dipolar EPR spectroscopy -Phys. Chem. Chem. Phys., 24 2022, 22645-22660 -""" -#%% -import numpy as np -import deerlab as dl - -# Load experimental data -files = [f'../data/triradical_protein_deer_{dB}dB.DTA' for dB in [0,6,9]] - -# Experiment information -t0 = 0.280 # Start time, μs -tau1 = 0.40 # First interpulse delay, μs -tau2 = 9.00 # Second interpulse delay, μs - -# Construct 4-pulse DEER experiment model -my4pdeer = dl.ex_4pdeer(tau1,tau2,pathways=[1]) - -# Loop over the different datasets -Vmodels,Vexps,ts,Vexps_sub,ts_sub = [],[],[],[],[] -for n,file in enumerate(files): - # Load the dataset - t,Vexp, descriptor = dl.deerload(file,full_output=True) - t = t[:-80] - Vexp = Vexp[:-80] - # Adjust the Start time - t = t - t[0] + t0 - - # Pre-processing - Vexp = dl.correctphase(Vexp) - Vexp /= np.max(Vexp) - - # Store the pre-processed datasets in a list - Vexps.append(Vexp) - ts.append(t) - - # Subsampling - # (required for efficient analysis in densely sampled datasets) - sampling = np.arange(0,len(t),4) # Take every 4th point - t_sub = t[sampling] - Vexp_sub = Vexp[sampling] - - # Store the subsampled datasets in a list - Vexps_sub.append(Vexp_sub) - ts_sub.append(t_sub) - - # Construct the three-spin dipolar model - Vmodel = dl.dipolarmodel(t_sub,spins=3,experiment=my4pdeer, minamp=0.01) - - # Add dipolar model to list of models - Vmodels.append(Vmodel) - -# Construct a global dipolar model describing all datasets -Vglobal = dl.merge(*Vmodels) -Vglobal = dl.link(Vglobal, - rmean1=[f'rmean1_{n+1}' for n in range(len(Vmodels))], - rmean2=[f'rmean2_{n+1}' for n in range(len(Vmodels))], - rmean3=[f'rmean3_{n+1}' for n in range(len(Vmodels))], - chol11=[f'chol11_{n+1}' for n in range(len(Vmodels))], - chol22=[f'chol22_{n+1}' for n in range(len(Vmodels))], - chol33=[f'chol33_{n+1}' for n in range(len(Vmodels))], - chol21=[f'chol21_{n+1}' for n in range(len(Vmodels))], - chol31=[f'chol31_{n+1}' for n in range(len(Vmodels))], - chol32=[f'chol32_{n+1}' for n in range(len(Vmodels))], - conc=[f'conc_{n+1}' for n in range(len(Vmodels))], - reftime1=[f'reftime1_{n+1}' for n in range(len(Vmodels))], -) -# Freeze the Cholesky-factors accounting for the correlation coefficients -# to zero (not always applicable) -Vglobal.chol21.freeze(0) -Vglobal.chol31.freeze(0) -Vglobal.chol32.freeze(0) - -# Fit the model to the data -results = dl.fit(Vglobal, Vexps_sub, reg=False, ftol=1e-5) - -# %% -# Plot the fitted datasets -results.plot(axis=ts_sub,xlabel='time (μs)') - -# Print the fit summary -print(results) - -# %% diff --git a/test/test_utils.py b/test/test_utils.py new file mode 100644 index 00000000..2df9a899 --- /dev/null +++ b/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])) \ No newline at end of file From 59a6584c8345084fb4ffda00a619ab3917f97ac0 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Sun, 14 Apr 2024 02:15:34 +0100 Subject: [PATCH 02/10] Removes katex and others (#471) * Fixing Sophgrid bug * Add unit test for sophgrid * Minor doc update * Bump version Number * Remove unnecessary doc files * Update changelog --- VERSION | 2 +- deerlab/fitresult.py | 2 ++ deerlab/solvers.py | 1 + deerlab/utils.py | 1 + docsrc/package-lock.json | 19 ------------------- docsrc/source/changelog.rst | 4 ++++ 6 files changed, 9 insertions(+), 20 deletions(-) delete mode 100644 docsrc/package-lock.json diff --git a/VERSION b/VERSION index 0f1acbd5..99a4aef0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v1.1.2 +v1.1.3 diff --git a/deerlab/fitresult.py b/deerlab/fitresult.py index 80c99bcc..4841f9e6 100644 --- a/deerlab/fitresult.py +++ b/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] diff --git a/deerlab/solvers.py b/deerlab/solvers.py index 1052d7b3..523afbe1 100644 --- a/deerlab/solvers.py +++ b/deerlab/solvers.py @@ -462,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 diff --git a/deerlab/utils.py b/deerlab/utils.py index 3b7f69d1..eb06134d 100644 --- a/deerlab/utils.py +++ b/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) diff --git a/docsrc/package-lock.json b/docsrc/package-lock.json deleted file mode 100644 index 02c9db5b..00000000 --- a/docsrc/package-lock.json +++ /dev/null @@ -1,19 +0,0 @@ -{ - "requires": true, - "lockfileVersion": 1, - "dependencies": { - "commander": { - "version": "2.20.0", - "resolved": "https://registry.npmjs.org/commander/-/commander-2.20.0.tgz", - "integrity": "sha512-7j2y+40w61zy6YC2iRNpUe/NwhNyoXrYpHMrSunaMG64nRnaf96zO/KMQR4OyN/UnE5KLyEBnKHd4aG3rskjpQ==" - }, - "katex": { - "version": "0.11.0", - "resolved": "https://registry.npmjs.org/katex/-/katex-0.11.0.tgz", - "integrity": "sha512-RQsU3HSMjLW9AdPpi2zaBwM123goCbUcUbBJfmjcAdA982RgtEUNMmrf+3y8anGjgJLantcNLa/VSK73xztQBg==", - "requires": { - "commander": "^2.19.0" - } - } - } -} diff --git a/docsrc/source/changelog.rst b/docsrc/source/changelog.rst index 894c269f..252ef0fc 100644 --- a/docsrc/source/changelog.rst +++ b/docsrc/source/changelog.rst @@ -24,6 +24,10 @@ 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 + Release ``v1.1.2`` - November 2023 ------------------------------------------ - |fix| : Fixes an issue with sophgrid (:pr:`463`). From e5cb1e2a3eb64e7be2d255f78c507cb4261f70cb Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 12 Jun 2024 19:31:44 +0200 Subject: [PATCH 03/10] Dipolarkernal speedup (#473) * Seperate Kinterpolator into its own function Interpolation in Scipy is very slow, currently the same interpolation is being rerun for every calculation of the dipolarkernal. This is now cached to speed it up. * Only run orientation selection on grid and integral based kernals This hunk of code is not needed when fresnel integrals are used so does not need to be evaluated * Update changelog --- deerlab/dipolarkernel.py | 29 ++++++++++++++++++----------- docsrc/source/changelog.rst | 1 + 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/deerlab/dipolarkernel.py b/deerlab/dipolarkernel.py index 33d815f1..8ee3d387 100644 --- a/deerlab/dipolarkernel.py +++ b/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): #------------------------------------------------------------------------ diff --git a/docsrc/source/changelog.rst b/docsrc/source/changelog.rst index 252ef0fc..3aa92232 100644 --- a/docsrc/source/changelog.rst +++ b/docsrc/source/changelog.rst @@ -27,6 +27,7 @@ Release Notes 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. Release ``v1.1.2`` - November 2023 ------------------------------------------ From 56689462e484504ce5561e19659ba58fcd499d3c Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 1 Jul 2024 17:21:18 +0200 Subject: [PATCH 04/10] Adding support for scipy 1.14 and tests for Python 3.12 (#474) * Update for 3.12 * Updated workflows and changelog * Upload to 3.12 --- .github/workflows/ci_PR.yml | 2 +- .github/workflows/ci_scheduled.yml | 2 +- .github/workflows/package_upload.yml | 4 ++-- deerlab/diststats.py | 2 +- docsrc/source/changelog.rst | 1 + 5 files changed, 6 insertions(+), 5 deletions(-) diff --git a/.github/workflows/ci_PR.yml b/.github/workflows/ci_PR.yml index 63de2eec..88374ad6 100644 --- a/.github/workflows/ci_PR.yml +++ b/.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 diff --git a/.github/workflows/ci_scheduled.yml b/.github/workflows/ci_scheduled.yml index 0659c65f..b101659c 100644 --- a/.github/workflows/ci_scheduled.yml +++ b/.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 diff --git a/.github/workflows/package_upload.yml b/.github/workflows/package_upload.yml index 69be3f25..4a3b9977 100644 --- a/.github/workflows/package_upload.yml +++ b/.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 diff --git a/deerlab/diststats.py b/deerlab/diststats.py index 9368de5b..3497e7b5 100644 --- a/deerlab/diststats.py +++ b/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""" diff --git a/docsrc/source/changelog.rst b/docsrc/source/changelog.rst index 3aa92232..d8edd8c8 100644 --- a/docsrc/source/changelog.rst +++ b/docsrc/source/changelog.rst @@ -28,6 +28,7 @@ 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 ------------------------------------------ From 6a0b992516dd53b6a8af79300f2520b3d5db14aa Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Tue, 2 Jul 2024 07:57:54 +0200 Subject: [PATCH 05/10] Numpy2 support (#479) * Update for 3.12 * Updated workflows and changelog * Upload to 3.12 * Numpy 2.0 compatibility updates * Update Changelog --- deerlab/classes.py | 2 +- deerlab/dd_models.py | 2 +- deerlab/diststats.py | 2 +- docsrc/source/changelog.rst | 4 +++- test/test_dipolarkernel.py | 2 +- 5 files changed, 7 insertions(+), 5 deletions(-) diff --git a/deerlab/classes.py b/deerlab/classes.py index 7e9b3e14..39063cbf 100644 --- a/deerlab/classes.py +++ b/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]) diff --git a/deerlab/dd_models.py b/deerlab/dd_models.py index dc32cbaa..1c4d6d3f 100644 --- a/deerlab/dd_models.py +++ b/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 diff --git a/deerlab/diststats.py b/deerlab/diststats.py index 3497e7b5..22d5d022 100644 --- a/deerlab/diststats.py +++ b/deerlab/diststats.py @@ -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 diff --git a/docsrc/source/changelog.rst b/docsrc/source/changelog.rst index d8edd8c8..c1aa84ca 100644 --- a/docsrc/source/changelog.rst +++ b/docsrc/source/changelog.rst @@ -28,7 +28,9 @@ 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 +- |fix| : Add support for Python 3.12 +- |fix| : Adds support for Numpy 2.0 + Release ``v1.1.2`` - November 2023 ------------------------------------------ diff --git a/test/test_dipolarkernel.py b/test/test_dipolarkernel.py index 8f7c505a..2d198541 100644 --- a/test/test_dipolarkernel.py +++ b/test/test_dipolarkernel.py @@ -1,6 +1,6 @@ import numpy as np -from numpy import pi, inf, NaN +from numpy import pi, inf from deerlab.bg_models import bg_hom3d,bg_exp from deerlab.dd_models import dd_gauss from deerlab.dipolarkernel import dipolarkernel,elementarykernel_twospin From 19a491fbe0783d61faeda5a37e0c6820074737c4 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 15 Jul 2024 13:50:49 +0200 Subject: [PATCH 06/10] Regparam grid bug fix (#477) * Update for 3.12 * Updated workflows and changelog * Upload to 3.12 * Fix bug in regparam grid search Regparam would never build the grid correctly. Now using grid or Brent is automatically determined from number of elements in the regparamrange. * Add extra error messages * Update changelog * Updated Example * Updated test The test has been updated. The previous convergence criteria was unreliable and only worked based on a coincidence. * Prepare For Release * Remove duplicate python version --- .github/workflows/ci_scheduled.yml | 2 +- LICENSE | 2 +- deerlab/selregparam.py | 41 ++++++++++++------- deerlab/solvers.py | 4 +- docsrc/source/changelog.rst | 5 ++- examples/intermediate/ex_selregparam.py | 54 ++++++++++++++++++++++++- test/test_selregparam.py | 9 +++-- 7 files changed, 91 insertions(+), 26 deletions(-) diff --git a/.github/workflows/ci_scheduled.yml b/.github/workflows/ci_scheduled.yml index b101659c..0659c65f 100644 --- a/.github/workflows/ci_scheduled.yml +++ b/.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","3.12"] + python-version: [3.8, 3.9, "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v3 diff --git a/LICENSE b/LICENSE index a2fb8b56..1c241f0e 100644 --- a/LICENSE +++ b/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 diff --git a/deerlab/selregparam.py b/deerlab/selregparam.py index 27a800cf..fffe61f9 100644 --- a/deerlab/selregparam.py +++ b/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'.") diff --git a/deerlab/solvers.py b/deerlab/solvers.py index 523afbe1..5fa94013 100644 --- a/deerlab/solvers.py +++ b/deerlab/solvers.py @@ -602,8 +602,8 @@ def linear_problem(y,A,optimize_alpha,alpha): if optimize_alpha: 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) + 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] diff --git a/docsrc/source/changelog.rst b/docsrc/source/changelog.rst index c1aa84ca..cbd4db7c 100644 --- a/docsrc/source/changelog.rst +++ b/docsrc/source/changelog.rst @@ -24,12 +24,13 @@ 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 +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. -- |fix| : Add support for Python 3.12 +- |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 diff --git a/examples/intermediate/ex_selregparam.py b/examples/intermediate/ex_selregparam.py index 417c284d..929d851c 100644 --- a/examples/intermediate/ex_selregparam.py +++ b/examples/intermediate/ex_selregparam.py @@ -33,7 +33,7 @@ t = t + tmin # Distance vector -r = np.linspace(1.5,7,50) # nm +r = np.linspace(1.5,7,100) # nm # Construct the model Vmodel = dl.dipolarmodel(t,r, experiment=dl.ex_4pdeer(tau1,tau2, pathways=[1])) @@ -101,7 +101,7 @@ # %% # Over and Under selection of the regularisation parameter -# -------------------------------------------------------- +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Here we will demonstrate the effect of selecting a regularisation parameter # that is either too small or too large. @@ -145,8 +145,58 @@ # As we can see when the regularisation parameter is too small we still get a high # quality fit in the time domain, however, our distance domain data is now way too # spikey and non-physical. + # In contrast when the regularisation parameter is too large we struggle to get # a good fit, however, we get a much smoother distance distribution. # This could have been seen from the selection functional above. The effect of # lower regularisation parameter had a smaller effect on the functional than the # effect of going to a larger one. + +# Plotting the full L-Curve +# ++++++++++++++++++++++++++ +# Normally the selection of the regularisation parameter is done through a Brent optimisation +# algorithm. This results in the L-Curve being sparsely sampled. However, the full L-Curve can be +# generated by evaluating the functional at specified regularisation paramters. + +# This has the disadvantage of being computationally expensive, however, it can be useful. + +regparamrange = 10**np.linspace(np.log10(1e-6),np.log10(1e1),60) # Builds the range of regularisation parameters +results_grid= dl.fit(Vmodel,Vexp,regparam='bic',regparamrange=regparamrange) +print(results_grid) + +fig, axs =plt.subplots(1,3, figsize=(9,4),width_ratios=(1,1,0.1)) +alphas = results_grid.regparam_stats['alphas_evaled'][:] +funcs = results_grid.regparam_stats['functional'][:] + + +idx = np.argsort(alphas) + +axs[0].semilogx(alphas[idx], funcs[idx],marker='.',ls='-') +axs[0].set_title(r"$\alpha$ selection functional"); +axs[0].set_xlabel("Regularisation Parameter") +axs[0].set_ylabel("Functional Value ") + +# Just the final L-Curve +x = results_grid.regparam_stats['residuals'] +y = results_grid.regparam_stats['penalties'] +idx = np.argsort(x) + + +axs[1].loglog(x[idx],y[idx]) + +n_points = results_grid.regparam_stats['alphas_evaled'].shape[-1] +lams = results_grid.regparam_stats['alphas_evaled'] +norm = mpl.colors.LogNorm(vmin=lams[:].min(), vmax=lams.max()) +for i in range(n_points): + axs[1].plot(x[i], y[i],marker = '.', ms=8, color=cmap(norm(lams[i]))) + +i_optimal = np.argmin(np.abs(lams - results_grid.regparam)) +axs[1].annotate(fr"$\alpha =$ {results_grid.regparam:.2g}", xy = (x[i_optimal],y[i_optimal]),arrowprops=dict(facecolor='black', shrink=0.05, width=5), xytext=(20, 20),textcoords='offset pixels') +fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap),cax=axs[2]) +axs[1].set_ylabel("Penalties") +axs[2].set_ylabel("Regularisation Parameter") +axs[1].set_xlabel("Residuals") +axs[1].set_title("L-Curve"); +fig.tight_layout() + +# %% diff --git a/test/test_selregparam.py b/test/test_selregparam.py index 30a64ee9..74ad7fd2 100644 --- a/test/test_selregparam.py +++ b/test/test_selregparam.py @@ -102,10 +102,11 @@ def test_unconstrained(dataset,design_matrix,regularization_matrix): #======================================================================= def test_manual_candidates(dataset,design_matrix,regularization_matrix): "Check that the alpha-search range can be manually passed" - alphas = np.linspace(-8,2,60) - alpha_manual = np.log10(selregparam(dataset,design_matrix,cvxnnls,method='aic',candidates=alphas,regop=regularization_matrix)) - alpha_auto = np.log10(selregparam(dataset,design_matrix,cvxnnls,method='aic',regop=regularization_matrix)) - assert abs(alpha_manual-alpha_auto)<1e-4 + alphas = np.logspace(-8,3,60,base=10) + alpha_manual = np.log10(selregparam(dataset,design_matrix,cvxnnls,method='aic',searchrange=alphas,regop=regularization_matrix)) + alpha_auto = np.log10(selregparam(dataset,design_matrix,cvxnnls,method='aic',searchrange=(1e-8,1e-3), regop=regularization_matrix)) + + assert abs(alpha_manual-alpha_auto)<1 #======================================================================= #======================================================================= From dbae0ba2047623845cb318dd136ac64f83c5a205 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 2 Sep 2024 21:59:52 +0200 Subject: [PATCH 07/10] Remove 3.8 Require Numpy 2.0 --- .github/workflows/ci_scheduled.yml | 2 +- README.md | 4 ++-- setup.py | 8 ++++---- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/.github/workflows/ci_scheduled.yml b/.github/workflows/ci_scheduled.yml index 0659c65f..65542133 100644 --- a/.github/workflows/ci_scheduled.yml +++ b/.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.9, "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v3 diff --git a/README.md b/README.md index 3ade50f1..2442dea3 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ The early versions of DeerLab (up to version 0.9.2) are written in MATLAB. The o ## Requirements -DeerLab is available for Windows, Mac and Linux systems and requires **Python 3.8**, **3.9**, **3.10**, or **3.11**. +DeerLab is available for Windows, Mac and Linux systems and requires **Python 3.9**, **3.10**, **3.11**, or **3.12**. All additional dependencies are automatically downloaded and installed during the setup. @@ -58,4 +58,4 @@ Here is the citation in bibtex format: DeerLab is licensed under the [MIT License](LICENSE). -Copyright © 2019-2023: Luis Fábregas Ibáñez, Stefan Stoll, Gunnar Jeschke, and [other contributors](https://github.com/JeschkeLab/DeerLab/contributors). +Copyright © 2019-2024: Luis Fábregas Ibáñez, Stefan Stoll, Gunnar Jeschke, and [other contributors](https://github.com/JeschkeLab/DeerLab/contributors). diff --git a/setup.py b/setup.py index ad8befbe..a0747967 100644 --- a/setup.py +++ b/setup.py @@ -11,7 +11,7 @@ 'Documentation': 'https://jeschkelab.github.io/DeerLab/', 'Source': 'https://github.com/JeschkeLab/DeerLab', }, - python_requires='>=3.8', + python_requires='>=3.9', license='LICENSE.txt', include_package_data = True, keywords='data analysis modeling least-squares EPR spectroscopy DEER PELDOR'.split(), @@ -19,9 +19,9 @@ long_description=open('README.md', encoding='utf-8').read(), long_description_content_type="text/markdown", install_requires = [ - 'numpy>=1.22.4', + 'numpy>=2.0', 'cvxopt>=1.0.0', - 'scipy>=1.6.3', + 'scipy>=1.11', 'joblib>=1.0.0', 'dill>=0.3.0', 'tqdm>=4.51.0', @@ -39,10 +39,10 @@ 'Operating System :: Microsoft :: Windows', 'Operating System :: MacOS', 'Operating System :: POSIX :: Linux', - 'Programming Language :: Python :: 3.8', 'Programming Language :: Python :: 3.9', 'Programming Language :: Python :: 3.10', 'Programming Language :: Python :: 3.11', + 'Programming Language :: Python :: 3.12', 'Topic :: Scientific/Engineering :: Chemistry', ] ) From ecfb6b13df28a558b37e26d204b7683d45f215d0 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 2 Sep 2024 22:25:53 +0200 Subject: [PATCH 08/10] Test 3.13 --- .github/workflows/ci_scheduled.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci_scheduled.yml b/.github/workflows/ci_scheduled.yml index 65542133..3d8839a9 100644 --- a/.github/workflows/ci_scheduled.yml +++ b/.github/workflows/ci_scheduled.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: [3.9, "3.10", "3.11", "3.12"] + python-version: [3.9, "3.10", "3.11", "3.12","3.13"] steps: - uses: actions/checkout@v3 From d633d6d9eb813d4f92f43c3695d490d98ec07ed1 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 2 Sep 2024 22:28:52 +0200 Subject: [PATCH 09/10] roll bacm --- .github/workflows/ci_scheduled.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci_scheduled.yml b/.github/workflows/ci_scheduled.yml index 3d8839a9..65542133 100644 --- a/.github/workflows/ci_scheduled.yml +++ b/.github/workflows/ci_scheduled.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: [3.9, "3.10", "3.11", "3.12","3.13"] + python-version: [3.9, "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v3 From 8d609c9d2af877aeb367c849254aaa23a8bc1828 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 2 Sep 2024 22:29:25 +0200 Subject: [PATCH 10/10] Roll back 3.13 Can't add 3.13 until official release