Fix for numpy 2.x

Author: AnkitBarik

python/magic/checker.py

@@ -2,9 +2,13 @@
 from magic import MagicTs, avgField, MagicSpectrum, MagicRadial, scanDir, MagicSetup
 import numpy as np
 from scipy.signal import argrelextrema
-import matplotlib.pyplot as plt
 import os
 
+try:
+    from numpy import trapz as trapz
+except:
+    from numpy import trapezoid as trapz
+
 class bcolors:
     HEADER = '\033[95m'
     OKBLUE = '\033[34m'
@@ -76,7 +80,7 @@ def MagicCheck(tstart=None):
         print(bcolors.MODERATE + 'Sudden variations detected in power balance!' + bcolors.ENDC)
         ones = np.ones_like(ts.time[mask])
         ones[~mask_spike] = 0.
-        ttot_spikes = np.trapz(ones, ts.time[mask])
+        ttot_spikes = trapz(ones, ts.time[mask])
         ttot = ts.time[-1]-ts.time[0]
         ratio = ttot_spikes/ttot
         print('    -Time fraction with spikes: {:.3f} %'.format(100.*ratio))

python/magic/coeff.py

@@ -2,14 +2,17 @@
 from magic import npfile, scanDir, MagicSetup, hammer2cart, symmetrize, progressbar
 from scipy.interpolate import interp1d
 from scipy import signal
-from scipy.version import version
 import os, re
 import numpy as np
 import matplotlib.pyplot as plt
 from .spectralTransforms import SpectralTransforms
 from magic.setup import labTex
 import copy
 
+try:
+    from numpy import trapz as trapz
+except:
+    from numpy import trapezoid as trapz
 
 def deriv(x, y, axis=0):
     """
@@ -224,8 +227,8 @@ def __init__(self, tag=None, datadir='.', ratio_cmb_surface=1, scale_b=1,
         # Time-averaged Gauss coefficient
         if not ave:
             facT = 1./(self.time[-1]-self.time[0])
-            self.glmM = facT * np.trapz(self.glm, self.time, axis=0)
-            self.hlmM = facT * np.trapz(self.hlm, self.time, axis=0)
+            self.glmM = facT * trapz(self.glm, self.time, axis=0)
+            self.hlmM = facT * trapz(self.hlm, self.time, axis=0)
 
             if len(self.time) > 3:
                 self.dglmdt = deriv(self.time, self.glm.T, axis=1)
@@ -257,11 +260,11 @@ def __init__(self, tag=None, datadir='.', ratio_cmb_surface=1, scale_b=1,
 
         if not ave:
             # Time-averaged energy
-            self.ElM = facT * np.trapz(self.El, self.time, axis=0)
-            self.EmM = facT * np.trapz(self.Em, self.time, axis=0)
+            self.ElM = facT * trapz(self.El, self.time, axis=0)
+            self.EmM = facT * trapz(self.Em, self.time, axis=0)
 
             # Secular variation
-            self.ESVlM = facT * np.trapz(self.ESVl, self.time, axis=0)
+            self.ESVlM = facT * trapz(self.ESVl, self.time, axis=0)
             if abs(self.ESVlM[1:]).min() > 0.:
                 self.taul = np.sqrt(self.ElM[1:]/self.ESVlM[1:])
             else:
@@ -838,10 +841,10 @@ def __init__(self, tag, datadir='.', ratio_cmb_surface=1, scale_b=1, iplot=True,
             # Time-averaged energy
             facT = 1./(self.time[-1]-self.time[0])
 
-            self.e_pol_lM = facT * np.trapz(self.e_pol_l, self.time, axis=0)
-            self.e_tor_lM = facT * np.trapz(self.e_tor_l, self.time, axis=0)
-            self.e_pol_axi_lM = facT * np.trapz(self.e_pol_axi_l, self.time, axis=0)
-            self.e_tor_axi_lM = facT * np.trapz(self.e_tor_axi_l, self.time, axis=0)
+            self.e_pol_lM = facT * trapz(self.e_pol_l, self.time, axis=0)
+            self.e_tor_lM = facT * trapz(self.e_tor_l, self.time, axis=0)
+            self.e_pol_axi_lM = facT * trapz(self.e_pol_axi_l, self.time, axis=0)
+            self.e_tor_axi_lM = facT * trapz(self.e_tor_axi_l, self.time, axis=0)
 
     def truncate(self, lCut):
         """

python/magic/cyl.py

@@ -11,6 +11,11 @@
 import os
 import pickle
 
+try:
+    from numpy import trapz as trapz
+except:
+    from numpy import trapezoid as trapz
+
 if buildSo:
     try:
         from magic.cylavg import *
@@ -234,7 +239,7 @@ def zavg(input, radius, ns, minc, save=True, filename='vp.pickle',
                 Z, S, out2D = sph2cyl_plane([input[0][iphi, ...]], radius, ns)
                 S = S[:, 1:-1]
                 Z = Z[:, 1:-1]
-                output[iphi, :] = np.trapz(out2D[0][:, 1:-1], z, axis=0)
+                output[iphi, :] = trapz(out2D[0][:, 1:-1], z, axis=0)
                 if normed:
                     output[iphi, :] /= height
 
@@ -252,7 +257,7 @@ def zavg(input, radius, ns, minc, save=True, filename='vp.pickle',
             output = []
             outIntZ = np.zeros((ns-2), dtype=input[0].dtype)
             for k,out in enumerate(out2D):
-                outIntZ = np.trapz(out[:, 1:-1], z, axis=0)
+                outIntZ = trapz(out[:, 1:-1], z, axis=0)
                 if normed:
                     outIntZ /= height
                 output.append(outIntZ)

python/magic/libmagic.py

@@ -9,6 +9,10 @@
 except:
     from scipy.integrate import simpson as simps
 
+try:
+    from numpy import trapz as trapz
+except:
+    from numpy import trapezoid as trapz
 
 def selectField(obj, field, labTex=True, ic=False):
     """
@@ -313,9 +317,9 @@ def avgField(time, field, tstart=None, std=False, fix_missing_series=False,
                 stdField = 0
         else:
             fac = 1./(time[ind2]-time[ind1])
-            avgField = fac*np.trapz(field[ind1:ind2+1, ...], time[ind1:ind2+1], axis=0)
+            avgField = fac*trapz(field[ind1:ind2+1, ...], time[ind1:ind2+1], axis=0)
             if std:
-                stdField = np.sqrt(fac*np.trapz((field[ind1:ind2+1, ...]-avgField)**2,
+                stdField = np.sqrt(fac*trapz((field[ind1:ind2+1, ...]-avgField)**2,
                                    time[ind1:ind2+1], axis=0))
 
     if std:
@@ -344,9 +348,9 @@ def writeVpEq(par, tstart):
     mask = np.where(abs(par.time-tstart) == min(abs(par.time-tstart)), 1, 0)
     ind = np.nonzero(mask)[0][0]
     fac = 1./(par.time.max()-par.time[ind])
-    avgReEq = fac*np.trapz(par.reEquat[ind:], par.time[ind:])
+    avgReEq = fac*trapz(par.reEquat[ind:], par.time[ind:])
     roEq = avgReEq*par.ek*(1.-par.radratio)
-    avgRolC = fac*np.trapz(par.rolc[ind:], par.time[ind:])
+    avgRolC = fac*trapz(par.rolc[ind:], par.time[ind:])
     st = '{:10.3e}{:5.2f}{:6.2f}{:11.3e}{:11.3e}{:11.3e}'.format(par.ek,
         par.strat, par.pr, par.ra, roEq, avgRolC)