Prepare1.2

.github/workflows/ci_PR.yml

@@ -12,13 +12,13 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ['3.11','3.12']
+        python-version: ['3.12','3.13']
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
           cache: 'pip'

.github/workflows/ci_scheduled.yml

@@ -12,12 +12,12 @@ 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.10", "3.11", "3.12","3.13"]
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
           cache: 'pip'

.github/workflows/deploy_ghpages.yml

@@ -10,12 +10,12 @@ jobs:
   deploy:
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v4
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v5
       with:
-        python-version: '3.10'
+        python-version: '3.12'
         cache: 'pip'
         cache-dependency-path: setup.py
     - name: Install dependencies

.github/workflows/docs_PR.yml

@@ -14,12 +14,12 @@ jobs:
   docsbuild:
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v4
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v5
       with:
-        python-version: '3.10'
+        python-version: '3.12'
         cache: 'pip'
         cache-dependency-path: setup.py   
     - name: Install dependencies

.github/workflows/examples_PR.yml

@@ -14,12 +14,12 @@ jobs:
   deploy:
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v4
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v5
       with:
-        python-version: '3.10'
+        python-version: '3.12'
         cache: 'pip'
         cache-dependency-path: setup.py
 

.github/workflows/package_upload.yml

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

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.9**, **3.10**, **3.11**, or  **3.12**.
+DeerLab is available for Windows, Mac and Linux systems and requires **Python 3.9** to **3.13**
 
 All additional dependencies are automatically downloaded and installed during the setup.
 

VERSION

@@ -1 +1 @@
-v1.1.4
+v1.2.0

deerlab/bg_models.py

@@ -188,7 +188,7 @@ def _hom3dex(t,conc,rex,lam):
     # Averaging integral
     z = np.linspace(0,1,1000)[np.newaxis,:]
     Dt = D*t[:,np.newaxis]*1e-6
-    Is =  4*np.pi/3*np.trapz(Dt*(1-3*z**2)*sici((Dt*(1-3*z**2))/((rex*1e-9)**3))[0],z,axis=1)
+    Is =  4*np.pi/3*np.trapezoid(Dt*(1-3*z**2)*sici((Dt*(1-3*z**2))/((rex*1e-9)**3))[0],z,axis=1)
 
     # Background function
     C_k = -Vex + Is + np.squeeze(Vex*(dipolarkernel(t,rex,integralop=False)))
@@ -242,7 +242,7 @@ def _hom3dex_phase(t,conc,rex,lam):
     ξ = 8*pi**2/9/np.sqrt(3)*(np.sqrt(3)+np.log(2-np.sqrt(3)))/np.pi*D
     z = np.linspace(0,1,1000)[np.newaxis,:]
     Dt = D*t[:,np.newaxis]*1e-6
-    Ic =  -ξ*(t*1e-6) + 4*np.pi/3*np.trapz(Dt*(1-3*z**2)*sici((Dt*np.abs(1-3*z**2))/((rex*1e-9)**3))[1],z,axis=1)
+    Ic =  -ξ*(t*1e-6) + 4*np.pi/3*np.trapezoid(Dt*(1-3*z**2)*sici((Dt*np.abs(1-3*z**2))/((rex*1e-9)**3))[1],z,axis=1)
 
     # Background function
     C_k = - Ic - np.squeeze(Vex*(dipolarkernel(t,rex,integralop=False,complex=True)).imag)

deerlab/classes.py

@@ -151,8 +151,8 @@ def __init__(self,uqtype,data=None,covmat=None,lb=None,ub=None,threshold=None,pr
         elif uqtype == 'profile':
             xs = [self.pardist(n)[0] for n in range(nParam)]
             pardists = [self.pardist(n)[1] for n in range(nParam)]
-            means = [np.trapz(pardist*x,x) for x,pardist in zip(xs,pardists)]
-            std = [np.sqrt(np.trapz(pardist*(x-mean)**2,x)) for x,pardist,mean in zip(xs,pardists,means)]
+            means = [np.trapezoid(pardist*x,x) for x,pardist in zip(xs,pardists)]
+            std = [np.sqrt(np.trapezoid(pardist*(x-mean)**2,x)) for x,pardist,mean in zip(xs,pardists,means)]
             self.mean = means
             self.median = self.percentile(50)
             self.std = std
@@ -342,8 +342,8 @@ def pardist(self,n=0):
 
         # Ensure normalization of the probability density function (if not a Dirac delta function)
         if not isdelta:
-            if np.trapz(pdf, x)!=0:
-                pdf = pdf/np.trapz(pdf, x)
+            if np.trapezoid(pdf, x)!=0:
+                pdf = pdf/np.trapezoid(pdf, x)
 
         return x, pdf
     #--------------------------------------------------------------------------------

deerlab/dd_models.py

@@ -50,16 +50,17 @@ def _dd_docstring(model,notes):
     string += '\n'
     string += '\n'
     table = []
-    table.append(['Name','Lower','Upper','Type','Frozen','Unit','Description'])  
+    table.append(['Name','Lower','Upper','par0','Type','Frozen','Unit','Description'])  
     for n,paramname in enumerate(model._parameter_list(order='vector')): 
         param_str = f'``{paramname}``'
         lb_str = f'{np.atleast_1d(getattr(model,paramname).lb)[0]:5.3g}'
         ub_str = f'{np.atleast_1d(getattr(model,paramname).ub)[0]:5.3g}'
+        par0_str = f'{np.atleast_1d(getattr(model,paramname).par0)[0]:5.3g}'
         linear_str = "linear" if np.all(getattr(model,paramname).linear) else "nonlin"
         frozen_str = "Yes" if np.all(getattr(model,paramname).frozen) else "No"
         unit_str = str(getattr(model,paramname).unit)
         desc_str = str(getattr(model,paramname).description)
-        table.append([param_str,lb_str,ub_str,linear_str,frozen_str,unit_str,desc_str])
+        table.append([param_str,lb_str,ub_str,par0_str,linear_str,frozen_str,unit_str,desc_str])
     string += formatted_table(table)
     string += f'\n{notes}'
 
@@ -129,7 +130,7 @@ def _multirice3dfun(r,nu,sig):
     P[P<0] = 0
 
     # Normalization
-    P = np.squeeze(P)/np.sum([np.trapz(c,np.squeeze(r)) for c in P.T])
+    P = np.squeeze(P)/np.sum([np.trapezoid(c,np.squeeze(r)) for c in P.T])
     return P
 # =================================================================
 
@@ -197,7 +198,7 @@ def _gauss2(r,mean1,std1,mean2,std2):
 #=======================================================================================
 #                                     dd_gauss3
 #=======================================================================================
-ntoes = r"""
+notes = r"""
 **Model**
 
 :math:`P(r) = a_1\frac{1}{\sigma_1\sqrt{2\pi}}\exp\left(-\frac{(r-\left<r_1\right>)^2}{2\sigma_1^2}\right) + a_2\frac{1}{\sigma_2\sqrt{2\pi}}\exp\left(-\frac{(r-\left<r_2\right>)^2}{2\sigma_2^2}\right) + a_3\frac{1}{\sigma_3\sqrt{2\pi}}\exp\left(-\frac{(r-\left<r_3\right>)^2}{2\sigma_3^2}\right)`

deerlab/dipolarmodel.py

@@ -642,8 +642,8 @@ def dipolarpenalty(Pmodel, r, type, selection=None):
         def compactness_penalty(*args): 
             P = Pmodel(*[r]*Nconstants,*args)
             if not np.all(P==0):
-                P = P/np.trapz(P,r)
-            return np.sqrt(P*(r - np.trapz(P*r,r))**2*np.mean(np.diff(r)))
+                P = P/np.trapezoid(P,r)
+            return np.sqrt(P*(r - np.trapezoid(P*r,r))**2*np.mean(np.diff(r)))
         # Add the penalty to the Pmodel
         penalty = Penalty(compactness_penalty,selection,
                     signature = Pmodel._parameter_list(),

deerlab/diststats.py

@@ -103,7 +103,7 @@ def analyze_rmode(V):
 
     # Auxiliary functions
     # -------------------
-    int = np.trapz(P,r) if not np.all(P==0) else 1
+    int = np.trapezoid(P,r) if not np.all(P==0) else 1
     def normalize(P): 
         return P/int
     # Percentile function
@@ -114,7 +114,7 @@ def pctile(r,P,p):
         return rpctile
     # Expectation operator function
     def E(x,P,r):
-        return np.trapz(x*normalize(P),r)
+        return np.trapezoid(x*normalize(P),r)
 
     # Location estimators
     # -------------------
@@ -125,7 +125,7 @@ def E(x,P,r):
     # Interquartile mean
     def iqmfcn(P):
         IQrange = (r>pctile(r,P,25)) & (r<pctile(r,P,75))
-        return E(r[IQrange],P[IQrange]/np.trapz(normalize(P)[IQrange],r[IQrange]),r[IQrange]) 
+        return E(r[IQrange],P[IQrange]/np.trapezoid(normalize(P)[IQrange],r[IQrange]),r[IQrange]) 
     # Mode
     modefcn = lambda P: r[np.argmax(P)]
     # Modes
@@ -156,7 +156,7 @@ def iqmfcn(P):
     exkurtosisfcn = lambda P: 3 - E(((r - meanfcn(P))/stdfcn(P))**4,P,r)
 
     # 0th moment  - Integral 
-    intfcn = lambda P: np.trapz(P,r)
+    intfcn = lambda P: np.trapezoid(P,r)
 
     # Calculate distribution estimators
     estimators = {

deerlab/solvers.py

@@ -794,7 +794,7 @@ def uq_subset(uq_full,subset,subset_lb,subset_ub):
         # Jacobian (non-linear part)
         Jnonlin = Jacobian(_ResidualsFcn,nonlinfit,lb,ub)
         # Jacobian (linear part)
-        scale = np.trapz(linfit,np.arange(Nlin))
+        scale = np.trapezoid(linfit,np.arange(Nlin))
         Jlin = (weights[:,np.newaxis]*Amodel(nonlinfit))[mask,:]
         if includeExtrapenalty:
             for penalty in extrapenalty:

docsrc/source/changelog.rst

@@ -24,6 +24,15 @@ 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.2`` - December 2024
+------------------------------------------
+- |feature| : `modelUncert` is now returned in the `FitResult` object
+- |fix| : Fixes issues with bootrstrap uncertainties
+- |api| : N bootstrap samples no longer produces N+1 samples
+- |enhancement| : Support Python 3.13
+- |fix| : All gaussian models now normalise to 1
+
 Release ``v1.1.4`` - September 2024
 ------------------------------------------
 - |enhancement| : Expanded sophgrid to allow for closed phi integral. (:pr:`482`)

examples/advanced/ex_forcefield_fit.py

@@ -38,7 +38,7 @@ def forcefield_P(c0,c1,c2,c3):
     # Boltzmann distribution
     Pr = np.exp(-energy/thermal)
     # Ensure a probability density distribution
-    Pr /= np.trapz(Pr,r)
+    Pr /= np.trapezoid(Pr,r)
     return Pr
 
 # File location

examples/advanced/ex_global_twostates_parametric.py

@@ -56,7 +56,7 @@ def Ptwostates(meanA, meanB, stdA, stdB, fracA):
     PA = fracA * dl.dd_gauss(r, meanA, stdA)
     PB = (1 - fracA) * dl.dd_gauss(r, meanB, stdB)
     P = PA + PB
-    P /= np.trapz(P)
+    P /= np.trapezoid(P)
     return P
 
 

examples/advanced/ex_multigauss_fitting_4pdeer.py

@@ -81,8 +81,8 @@
     # Calculate the 95%-confidence intervals
     Pci = Puq.ci(95)
     # Normalize the probability density functions
-    Pci /= np.trapz(Pfit, r)
-    Pfit /= np.trapz(Pfit, r)
+    Pci /= np.trapezoid(Pfit, r)
+    Pfit /= np.trapezoid(Pfit, r)
     # Plot the optimal fit with a thicker line
     if n == np.argmin(aic):
         lw = 4

examples/advanced/ex_pseudotitration_parameter_free.py

@@ -119,7 +119,7 @@ def chemicalequilibrium(Kdis,L):
 for n,(xA,xB) in enumerate(zip(xAfit,xBfit)): 
 
     Pfit = Pmodel(P_1=results.P_1,P_2=results.P_2,weight_1=xA,weight_2=xB)
-    Pfit /= np.trapz(Pfit,r)
+    Pfit /= np.trapezoid(Pfit,r)
     if n>1: label=None
     plt.plot(r,2*n + Pfit,'k',label='Total contribution' if n<1 else None)
     plt.fill(r,2*n + xA*results.P_1,color=green,alpha=0.5,label='State A (natural)' if n<1 else None)

examples/basic/ex_fitting_4pdeer_gauss.py

@@ -52,7 +52,7 @@
 
 # Extract fitted distance distribution
 Pfit = results.evaluate(Pmodel,r)
-scale = np.trapz(Pfit,r)
+scale = np.trapezoid(Pfit,r)
 Puncert = results.propagate(Pmodel,r,lb=np.zeros_like(r))
 Pfit = Pfit/scale
 Pci95 = Puncert.ci(95)/scale

examples/basic/ex_restraints_4pdeer.py

@@ -70,9 +70,9 @@
 
 # Plot distribution and confidence bands
 violet = '#4550e6'
-Pci95 = fit.PUncert.ci(95)/np.trapz(fit.P,r)
-Pci50 = fit.PUncert.ci(50)/np.trapz(fit.P,r)
-plt.plot(r,fit.P/np.trapz(fit.P,r),linewidth=2,color=violet,label='Distance distribution fit')
+Pci95 = fit.PUncert.ci(95)/np.trapezoid(fit.P,r)
+Pci50 = fit.PUncert.ci(50)/np.trapezoid(fit.P,r)
+plt.plot(r,fit.P/np.trapezoid(fit.P,r),linewidth=2,color=violet,label='Distance distribution fit')
 plt.fill_between(r,Pci95[:,0],Pci95[:,1],color=violet,alpha=0.3)
 plt.fill_between(r,Pci50[:,0],Pci50[:,1],color=violet,alpha=0.4)
 

test/test_ddmodels.py

@@ -30,7 +30,7 @@ def assert_ddmodel(model):
     assert all(P1 >= 0)
     assert all(P1 >= 0) and all(P2 >= 0)
     assert all(~np.isnan(P1)) and all(~np.isnan(P2)) and all(~np.isnan(P3)) and all(~np.isnan(P4))
-    assert np.round(np.trapz(P5,rnus),2) == 1
+    assert np.round(np.trapezoid(P5,rnus),2) == 1
     assert len(lower)==nParam
     assert len(upper)==nParam
     assert len(meta['names'])==nParam

test/test_dipolarmodel.py

@@ -320,7 +320,7 @@ def test_fit_Pnonparametric_normalization(V1path):
 
     result = fit(Vmodel,V1path,ftol=1e-5)
 
-    assert np.isclose(np.trapz(result.P,r),1)
+    assert np.isclose(np.trapezoid(result.P,r),1)
 # ======================================================================
 
 # Fixtures

test/test_model_penalty.py

@@ -24,8 +24,8 @@ def mock_data():
 def penalty_fcn(): 
     def _penalty_fcn(mean,std): 
         P = dd_gauss(x,mean,std)
-        P = P/np.trapz(P,x)
-        return np.sqrt(P*(x - np.trapz(P*x,x))**2*np.mean(np.diff(x)))
+        P = P/np.trapezoid(P,x)
+        return np.sqrt(P*(x - np.trapezoid(P*x,x))**2*np.mean(np.diff(x)))
     return _penalty_fcn
 # -----------------------------------------------------------------------
 

test/test_snlls.py

@@ -207,7 +207,7 @@ def test_SNLLS_cost_value(mock_data,mock_Amodel):
 def test_SNLLS_fit_with_extra_penalty(mock_data,mock_Amodel):
     "Check that an additional penalty can be passed correctly to the SNLLS functional"
     beta = 0.05
-    compactness_penalty = lambda _,plin: beta*np.sqrt(plin*(r - np.trapz(plin*r,r))**2*dr)
+    compactness_penalty = lambda _,plin: beta*np.sqrt(plin*(r - np.trapezoid(plin*r,r))**2*dr)
     fit = snlls(mock_data,mock_Amodel,nlpar0,lb,ub,lbl,extrapenalty=compactness_penalty)
     assert np.all(abs(lin_param - fit.lin) < 1e-1) and np.all(abs(nonlin_param - fit.nonlin[0]) < 1e-1)
 # ======================================================================
@@ -216,7 +216,7 @@ def test_SNLLS_fit_with_multiple_extra_penalties(mock_data,mock_Amodel):
 # ======================================================================
     "Check that multiple additional penaltyies can be passed correctly to the SNLLS functional"
     beta = 0.05
-    compactness_penalty = lambda _,plin: beta*np.sqrt(plin*(r - np.trapz(plin*r,r))**2*dr)
+    compactness_penalty = lambda _,plin: beta*np.sqrt(plin*(r - np.trapezoid(plin*r,r))**2*dr)
     fit = snlls(mock_data,mock_Amodel,nlpar0,lb,ub,lbl,extrapenalty=[compactness_penalty])
     R = 0.5
     radial_penalty = lambda pnonlin,_: 1/R**2*(np.linalg.norm((pnonlin-nonlin_param)/nonlin_param-R))**2