Compatability update

ares/__init__.py

@@ -1,13 +1,14 @@
 import os as _os
-import imp as _imp
+import importlib as _imp
+# updates to imp: https://docs.python.org/3/whatsnew/3.12.html#whatsnew312-removed-imp
 
 _HOME = _os.environ.get('HOME')
 
 # Load custom defaults    
-if _os.path.exists('{!s}/.ares/defaults.py'.format(_HOME)):
+if _os.path.exists('~/.ares/defaults.py'.format(_HOME)):
     (_f, _filename, _data) =\
-        _imp.find_module('defaults', ['{!s}/.ares/'.format(_HOME)])
-    rcParams = _imp.load_module('defaults.py', _f, _filename, _data).pf
+        _imp.util.find_spec('defaults', ['~/.ares/'.format(_HOME)])
+    rcParams = _imp.import_module('defaults.py', _f, _filename, _data).pf
 else:
     rcParams = {}
 

ares/analysis/MetaGalacticBackground.py

@@ -14,7 +14,7 @@
 from ..util import labels
 from ..util.Pickling import read_pickle_file
 import matplotlib.pyplot as pl
-from scipy.integrate import trapz
+from scipy.integrate import   trapezoid
 from ..util.ReadData import flatten_energies
 from ..physics.Constants import erg_per_ev, J21_num, h_P, c, E_LL, E_LyA, \
     sqdeg_per_std

ares/analysis/MultiPhaseMedium.py

@@ -15,9 +15,8 @@
 import matplotlib.pyplot as pl
 from ..util.Stats import get_nu
 from ..util.Pickling import read_pickle_file
-from scipy.misc import derivative
 from ..physics.Constants import *
-from scipy.integrate import cumtrapz
+from scipy.integrate import cumulative_trapezoid
 from scipy.interpolate import interp1d
 from ..physics import Cosmology, Hydrogen
 from ..util.SetDefaultParameterValues import *
@@ -341,7 +340,7 @@ def tau_CMB(self, include_He=True, z_HeII_EoR=3.):
 
         integrand *= sigma_T * dldz
 
-        tau = cumtrapz(integrand, self.history_asc['z'], initial=0)
+        tau = cumulative_trapezoid(integrand, self.history_asc['z'], initial=0)
 
         tau[self.history_asc['z'] > 100] = 0.0
 
@@ -803,7 +802,7 @@ def tau_post_EoR(self, include_He=True, z_HeII_EoR=3.):
 
         integrand *= sigma_T * dldz
 
-        tau = cumtrapz(integrand, ztmp, initial=0)
+        tau = cumulative_trapezoid(integrand, ztmp, initial=0)
 
         return ztmp, tau
 

ares/analysis/TurningPoints.py

@@ -11,8 +11,8 @@
 """
 
 import numpy as np
+import numdifftools as nd
 from ..util import ParameterFile
-from scipy.misc import derivative
 from scipy.optimize import minimize
 from ..physics.Constants import nu_0_mhz
 from ..util.Math import central_difference
@@ -202,9 +202,10 @@ def is_stopping_point(self, z, dTb):
                         continue
 
                 else:
+                    raise NotImplemented('havent revisited since scipy.misc.derivative Deprecation')
                     # Compute curvature at turning point (mK**2 / MHz**2)
                     nuTP = nu_0_mhz / (1. + zTP)
-                    d2 = float(derivative(lambda zz: splev(zz, Bspl_fit1),
+                    d2 = float(nd.Derivative(lambda zz: splev(zz, Bspl_fit1),
                         x0=float(zTP), n=2, dx=1e-4, order=5) * nu_0_mhz**2 / nuTP**4)
 
                     self.turning_points[TP] = (zTP, TTP, d2)

ares/phenom/Tanh21cm.py

@@ -12,8 +12,8 @@
 
 import time
 import numpy as np
+import numdifftools as nd
 from ..util import ParameterFile
-from scipy.misc import derivative
 from ..physics import Hydrogen, Cosmology
 from ..physics.Constants import k_B, J21_num, nu_0_mhz
 from ..physics.RateCoefficients import RateCoefficients
@@ -52,7 +52,7 @@ def __init__(self, **kwargs):
         self.hydr = Hydrogen(cosm=self.cosm, **kwargs)
 
     def dTgas_dz(self, z):
-        return derivative(self.cosm.Tgas, x0=z)
+        return nd.Derivative(self.cosm.Tgas)(z)
 
     def electron_density(self, z):
         return np.interp(z, self.cosm.thermal_history['z'],

ares/physics/Cosmology.py

@@ -12,7 +12,7 @@
 import os
 import numpy as np
 from ..data import ARES
-from scipy.misc import derivative
+import numdifftools as nd
 from scipy.optimize import fsolve
 from scipy.integrate import quad, ode
 from ..util.Math import interp1d
@@ -462,7 +462,7 @@ def cooling_rate(self, z, T=None):
             ##s
             #func = lambda zz: np.interp(zz, self.inits['z'], self.inits['Tk'])
 
-            dTdz = derivative(self._Tgas_CosmoRec, z, dx=1e-2)
+            dTdz = nd.Derivative(self._Tgas_CosmoRec)(z)
 
             xe = np.interp(z, self.inits['z'], self.inits['xe'])
 
@@ -483,7 +483,7 @@ def cooling_rate(self, z, T=None):
             return dTdz + xe_cool * mult
 
         else:
-            return derivative(self.Tgas, z)
+            return nd.Derivative(self.Tgas)(z)
 
     def log_cooling_rate(self, z):
         if self.pf['approx_thermal_history'] == 'exp':

ares/physics/DustEmission.py

@@ -12,7 +12,7 @@
 """
 
 import numpy as np
-from scipy.integrate import simps
+from scipy.integrate import simpson
 from ares.physics.Constants import c, h, k_B, g_per_msun, cm_per_kpc, Lsun
 
 # T_dust parameters
@@ -376,8 +376,8 @@ def __T_dust(self, z, L_nu, tau_nu, R_dust, T_cmb):
             tmp_cmb = 8 * np.pi * h / c**2 * cmb_kappa_nu * (cmb_freqs[None, :, None])**3 \
                 / (np.exp(h * cmb_freqs[None,:,None] / k_B / T_cmb[:,None,:]) - 1)
 
-            tmp_power = simps(tmp_stellar, self.frequencies, axis = 1)
-            tmp_power += simps(tmp_cmb, cmb_freqs, axis = 1)
+            tmp_power = simpson(tmp_stellar, self.frequencies, axis = 1)
+            tmp_power += simpson(tmp_cmb, cmb_freqs, axis = 1)
 
             if self.pf.get('pop_dust_experimental'):
                 print("power =", tmp_power)
@@ -519,7 +519,7 @@ def Luminosity(self, z, wave = 3e5, band=None, idnum=None, window=1,
                 fmax = c / (8 * 1e-4)
                 fmin = c / (1000 * 1e-4)
                 freqs, luminosities = self.dust_sed(fmin, fmax, 1000)
-                luminosities = simps(luminosities[:,:,index], freqs, axis = 1)
+                luminosities = simpson(luminosities[:,:,index], freqs, axis = 1)
 
         # is not cached, we calculate everything for the given z and wave
         else:
@@ -557,7 +557,7 @@ def Luminosity(self, z, wave = 3e5, band=None, idnum=None, window=1,
                 kappa_nu = 0.1 * (nu / 1e12)**2
                 luminosities = 8 * np.pi * h / c**2 * nu[None,:]**3 * kappa_nu[None,:] \
                     / (np.exp(h * nu[None,:] / k_B / T_dust[:,None]) - 1) * (M_dust[:,None] * g_per_msun)
-                luminosities = simps(luminosities, nu, axis = 1)
+                luminosities = simpson(luminosities, nu, axis = 1)
 
 
         if idnum is not None:

ares/physics/ExcursionSet.py

@@ -15,9 +15,8 @@
 from .Cosmology import Cosmology
 from ..util.Math import central_difference
 from ..util.ParameterFile import ParameterFile
-from scipy.integrate import simps, quad
+from scipy.integrate import simpson, quad
 from scipy.interpolate import interp1d
-from scipy.misc import derivative
 
 two_pi = 2. * np.pi
 four_pi = 4. * np.pi

ares/physics/HaloMassFunction.py

@@ -19,10 +19,9 @@
 from ..data import ARES
 from types import FunctionType
 from ..util import ParameterFile
-from scipy.misc import derivative
 from scipy.optimize import fsolve
 from ..util.Warnings import no_hmf
-from scipy.integrate import cumtrapz, simps
+from scipy.integrate import cumulative_trapezoid, simpson
 from ..util.PrintInfo import print_hmf
 from ..util.ProgressBar import ProgressBar
 from ..util.ParameterFile import ParameterFile
@@ -437,14 +436,14 @@ def _get_ngtm_mgtm_from_dndm(self):
                 mf_func = InterpolatedUnivariateSpline(np.log(m), np.log(dndlnm), k=1)
                 mf = mf_func(m_upper)
 
-                int_upper_n = simps(np.exp(mf), dx=m_upper[2] - m_upper[1], even='first')
-                int_upper_m = simps(np.exp(m_upper + mf), dx=m_upper[2] - m_upper[1], even='first')
+                int_upper_n = simpson(np.exp(mf), dx=m_upper[2] - m_upper[1], even='first')
+                int_upper_m = simpson(np.exp(m_upper + mf), dx=m_upper[2] - m_upper[1], even='first')
             else:
                 int_upper_n = 0
                 int_upper_m = 0
 
-            ngtm_ = np.concatenate((cumtrapz(dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
-            mgtm_ = np.concatenate((cumtrapz(m[::-1] * dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
+            ngtm_ = np.concatenate((cumulative_trapezoid(dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
+            mgtm_ = np.concatenate((cumulative_trapezoid(m[::-1] * dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
 
             ngtm.append(ngtm_ + int_upper_n)
             mgtm.append(mgtm_ + int_upper_m)
@@ -567,10 +566,10 @@ def _load_hmf(self):
                 mgtm_0 = np.trapz(self.tab_dndm[i] * self.tab_M**2,
                     x=np.log(self.tab_M))
                 self.tab_ngtm[i,:] = ngtm_0 \
-                    - cumtrapz(self.tab_dndm[i] * self.tab_M,
+                    - cumulative_trapezoid(self.tab_dndm[i] * self.tab_M,
                         x=np.log(self.tab_M), initial=0.0)
                 self.tab_mgtm[i,:] = mgtm_0 \
-                    - cumtrapz(self.tab_dndm[i] * self.tab_M**2,
+                    - cumulative_trapezoid(self.tab_dndm[i] * self.tab_M**2,
                         x=np.log(self.tab_M), initial=0.0)
 
             # Keep it positive please.

ares/physics/RateCoefficients.py

@@ -12,7 +12,7 @@
 """
 
 import numpy as np
-from scipy.misc import derivative
+import numdifftools as nd
 from ..util.Math import interp1d
 from ..util.Math import central_difference
 
@@ -71,7 +71,7 @@ def _dCollisionalIonizationRate(self):
         if not hasattr(self, '_dCollisionalIonizationRate_'):
             self._dCollisionalIonizationRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.CollisionalIonizationRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.CollisionalIonizationRate(i, T))(self.Tarr)
                 self._dCollisionalIonizationRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
 
@@ -132,7 +132,7 @@ def _dRadiativeRecombinationRate(self):
         if not hasattr(self, '_dRadiativeRecombinationRate_'):
             self._dRadiativeRecombinationRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.RadiativeRecombinationRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.RadiativeRecombinationRate(i, T))(self.Tarr)
 
                 self._dRadiativeRecombinationRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
@@ -161,7 +161,7 @@ def DielectricRecombinationRate(self, T):
     def _dDielectricRecombinationRate(self):
         if not hasattr(self, '_dDielectricRecombinationRate_'):
             self._dDielectricRecombinationRate_ = {}
-            tmp = derivative(lambda T: self.DielectricRecombinationRate(T), self.Tarr)
+            tmp = nd.Derivative(lambda T: self.DielectricRecombinationRate(T))(self.Tarr)
             self._dDielectricRecombinationRate_ = interp1d(self.Tarr, tmp, 
                 kind=self.interp_rc)
 
@@ -197,7 +197,7 @@ def _dCollisionalIonizationCoolingRate(self):
         if not hasattr(self, '_dCollisionalIonizationCoolingRate_'):
             self._dCollisionalIonizationCoolingRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T))(self.Tarr)
                 self._dCollisionalIonizationCoolingRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
 
@@ -233,7 +233,7 @@ def _dCollisionalExcitationCoolingRate(self):
         if not hasattr(self, '_dCollisionalExcitationCoolingRate_'):
             self._dCollisionalExcitationCoolingRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T))(self.Tarr)
                 self._dCollisionalExcitationCoolingRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
 
@@ -272,7 +272,7 @@ def _dRecombinationCoolingRate(self):
         if not hasattr(self, '_dRecombinationCoolingRate_'):
             self._dRecombinationCoolingRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.RecombinationCoolingRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.RecombinationCoolingRate(i, T))(self.Tarr)
                 self._dRecombinationCoolingRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
 
@@ -300,7 +300,7 @@ def DielectricRecombinationCoolingRate(self, T):
     @property 
     def _dDielectricRecombinationCoolingRate(self):
         if not hasattr(self, '_dDielectricRecombinationCoolingRate_'):
-            tmp = derivative(lambda T: self.DielectricRecombinationCoolingRate(T), self.Tarr)
+            tmp = nd.Derivative(lambda T: self.DielectricRecombinationCoolingRate(T))(self.Tarr)
             self._dDielectricRecombinationCoolingRate_ = interp1d(self.Tarr, tmp, 
                 kind=self.interp_rc)
 

ares/populations/GalaxyAggregate.py

@@ -18,7 +18,7 @@
 from .Halo import HaloPopulation
 from collections import namedtuple
 from ..util.Math import interp1d
-from scipy.integrate import quad, simps
+from scipy.integrate import quad, simpson
 from ..util.Warnings import negative_SFRD
 from ..util.ParameterFile import get_pq_pars, pop_id_num
 from scipy.interpolate import interp1d as interp1d_scipy