Numpy Fixes

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.1.5

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

@@ -129,7 +129,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
 # =================================================================
 

deerlab/dipolarkernel.py

@@ -403,7 +403,7 @@ def K0_3spin(tdip):
         if Nspin==2:
 
             # Determine the modulated dipolar interaction
-            q = int(np.where([np.any(δq!=0) for δq in δ])[0])
+            q = int(np.where([np.any(δq!=0) for δq in δ])[0][0])
             # Construc the effective dipolar evolution time for the modulated dipolar interactions
             tdip = np.sum(np.array([δ_qd*(t_d-tref_qd) for t_d,δ_qd,tref_qd in zip(t,δ[q],tref[q])],dtype=object), axis=0).astype(float)
             # Compute and accumulate the two-spin dipolar pathway contribution to the dipolar kernel

deerlab/dipolarmodel.py

@@ -240,7 +240,7 @@ def dipolarpathways(*param):
         for idx in range(npathways):
             # Construct all the permutations of the two-spin interaction pathways (without repetitions)
             for perm in set(set(permutations([0]+[None]*(Q-1)))):
-                q = int(np.where(np.array(perm)==0)[0])
+                q = int(np.where(np.array(perm)==0)[0][0])
                 # Compute amplitude of the two-spin interaction pathway 
                 λp = λs[q][idx]
                 λ2k = factorial(Q-1)*λp*λu**(Q-1)
@@ -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/model.py

@@ -1426,7 +1426,7 @@ def _linkparameter(model,parameters,newname):
                 else: 
                     model.Nnonlin -= Nremoved
                     # Update the parameter vector map
-                    mapping[unlinked_nonlinear_idx[n]] = link_indices
+                    mapping[unlinked_nonlinear_idx[n]] = link_indices[0]
                 # Delete the linked parameter from the model
                 delattr(model,param)
 

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,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.5`` - January 2025
+- |fix|: Moves to numpy 2.0 as a mimum requirement, and removes all `np.trapz` calls to `np.trapezoid`.
+
 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