Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion src/mechanalyzer/calculator/bf.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,12 @@
import numpy as np
import pandas as pd
from mechanalyzer import calculator
# AVC: Backward compatibility for numpy < 2.0, where the function is called `trapz`
try:
from numpy import trapezoid
except ImportError:
from numpy import trapz as trapezoid


######################## wrapper functions #########################################################

Expand Down Expand Up @@ -160,7 +166,7 @@ def bf_tp_df_full(ped_df, hotbf_df):
kind='cubic', fill_value=(hoten_spc.values[0], hoten_spc.values[-1]))
hoten_vect = f_hoten(ene_vect)
# recompute in an appropriate range
bf_series[spc] = np.trapz(ped_vect*hoten_vect, x=ene_vect)
bf_series[spc] = trapezoid(ped_vect*hoten_vect, x=ene_vect)
# renormalize for all species and put in dataframe
if any(bf_series.values < 0):
print('Warning: found negative BFs at {:1.0f} K and {:1.1e} atm'
Expand Down
23 changes: 14 additions & 9 deletions src/mechanalyzer/calculator/ene_partition.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,11 @@
from phydat import phycon
from mechanalyzer import calculator
from mechanalyzer.calculator.spinfo_frommess import get_dof_info_fromspcdct
# AVC: Backward compatibility for numpy < 2.0, where the function is called `trapz`
try:
from numpy import trapezoid
except ImportError:
from numpy import trapz as trapezoid


#################### wrapper function that calls the class ################################
Expand Down Expand Up @@ -144,7 +149,7 @@ def ped_df_rescale(starthot_df, ped_df_fromhot, save=False, name=''):
kernelsize)/kernelsize, mode='same')
cycles -= 1
# renormalize and put in dataframe
prob_vect /= np.trapz(prob_vect, x=ene_vect)
prob_vect /= trapezoid(prob_vect, x=ene_vect)
ped_df.at[temp, pressure] = pd.Series(prob_vect, index=ene_vect)
if ped_df[pressure][temp].empty:
T_del.append(temp)
Expand Down Expand Up @@ -280,7 +285,7 @@ def equip_simple(self):
for pressure in self.ped_df.columns:
for temp in self.ped_df.sort_index().index:
idx_new = self.ped_df[pressure][temp].index * beta_prod
norm_factor = np.trapz(
norm_factor = trapezoid(
self.ped_df[pressure][temp].values, x=idx_new)
vals = self.ped_df[pressure][temp].values/norm_factor
ped_df_prod.at[temp, pressure] = pd.Series(vals, index=idx_new)
Expand Down Expand Up @@ -348,7 +353,7 @@ def init_dos(pressure, temp):
rho_rovib_prod2[idx_ene_int] *
rho_trasl[idx_ene_minus_ene_int]
)
rho_non1.append(np.trapz(rho_non1_integrand,
rho_non1.append(trapezoid(rho_non1_integrand,
x=self.ene1_vect[idx_ene_int]))

rho_non1 = np.array(rho_non1)
Expand Down Expand Up @@ -378,7 +383,7 @@ def dos(idx_ene1, idx_ene_new):
# rhonon1(ene-ene1) with ene1<ene (fixed ene)
rho_non1_array = self.rho_non1[idx_ene_minus_ene1_array]
num = rho1_ene1 * rho_non1
den = np.trapz(rho1_ene1_array*rho_non1_array,
den = trapezoid(rho1_ene1_array*rho_non1_array,
x=self.ene1_vect[idx_ene1_array])
# if for some reason you get den=0: append 0 as value
if den == 0:
Expand Down Expand Up @@ -408,13 +413,13 @@ def dostherm_rhovibtrasl(idx_ene_vect, temp):
rho1_ene1_array = self.rho_rovib_prod1[idx_ene1_array]
rho_non1_array = self.rho_non1[idx_ene_minus_ene1_array]

rhotot_E1 = np.trapz(rho1_ene1_array*rho_non1_array,
rhotot_E1 = trapezoid(rho1_ene1_array*rho_non1_array,
x=self.ene1_vect[idx_ene1_array])

ene = self.ene1_vect[idx_ene1]
f_Etot_num.append(rhotot_E1*np.exp(-ene/phycon.RC_KCAL/temp))

den = np.trapz(f_Etot_num, x=self.ene1_vect[idx_ene_vect])
den = trapezoid(f_Etot_num, x=self.ene1_vect[idx_ene_vect])
f_Etot = pd.Series(
f_Etot_num/den, index=self.ene1_vect[idx_ene_vect])

Expand All @@ -425,7 +430,7 @@ def dostherm_rhovib(temp):
# rovib only
rho1_ene1_array = self.rho_rovib_prod1 * \
np.exp(-self.ene1_vect/phycon.RC_KCAL/temp)
den = np.trapz(rho1_ene1_array, x=self.ene1_vect)
den = trapezoid(rho1_ene1_array, x=self.ene1_vect)
f_Etot = pd.Series(rho1_ene1_array/den, index=self.ene1_vect)

return f_Etot
Expand Down Expand Up @@ -489,11 +494,11 @@ def dostherm_rhovib(temp):
ped_series.values[idx_ene_vect >= idx_ene1]
)

prob_ene1 = np.trapz(
prob_ene1 = trapezoid(
prob_ene1ene_tot_pressure_ped, ene_new)
prob_ene1_vect.append(prob_ene1)

norm_factor_prob_ene1 = np.trapz(
norm_factor_prob_ene1 = trapezoid(
prob_ene1_vect, x=self.ene1_vect)
prob_ene1_norm = prob_ene1_vect/norm_factor_prob_ene1

Expand Down
33 changes: 19 additions & 14 deletions src/mechanalyzer/tests/test__prompt.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,11 @@
from mechanalyzer.calculator import spinfo_frommess
from mechanalyzer.calculator import ene_partition
from mechanalyzer.calculator import bf
# AVC: Backward compatibility for numpy < 2.0, where the function is called `trapz`
try:
from numpy import trapezoid
except ImportError:
from numpy import trapz as trapezoid

PATH = os.path.dirname(os.path.realpath(__file__))
INP_PATH = os.path.join(PATH, 'data', 'prompt', 'C3H8_OH')
Expand Down Expand Up @@ -70,8 +75,8 @@ def test_equip_simple():

assert np.isclose((ped_600.iloc[100]), 0.1043, atol=1e-4, rtol=1e-4)
assert np.isclose((ped_1200.iloc[100]), 0.02105, atol=1e-4, rtol=1e-4)
assert np.isclose(np.trapz(ped_600.values, x=ped_600.index), 1)
assert np.isclose(np.trapz(ped_1200.values, x=ped_1200.index), 1)
assert np.isclose(trapezoid(ped_600.values, x=ped_600.index), 1)
assert np.isclose(trapezoid(ped_1200.values, x=ped_1200.index), 1)


def test_equip_phi():
Expand All @@ -89,8 +94,8 @@ def test_equip_phi():

assert np.isclose((ped_600.iloc[166]), 0.045, atol=1e-3, rtol=1e-2)
assert np.isclose((ped_1200.iloc[166]), 0.0376, atol=1e-3, rtol=1e-2)
assert np.isclose(np.trapz(ped_600.values, x=ped_600.index), 1)
assert np.isclose(np.trapz(ped_1200.values, x=ped_1200.index), 1)
assert np.isclose(trapezoid(ped_600.values, x=ped_600.index), 1)
assert np.isclose(trapezoid(ped_1200.values, x=ped_1200.index), 1)


def test_beta_phi1a():
Expand All @@ -108,8 +113,8 @@ def test_beta_phi1a():

assert np.isclose((ped_400.iloc[166]), 0.001, atol=1e-3, rtol=1e-2)
assert np.isclose((ped_800.iloc[166]), 0.0589, atol=1e-3, rtol=1e-2)
assert np.isclose(np.trapz(ped_400.values, x=ped_400.index), 1)
assert np.isclose(np.trapz(ped_800.values, x=ped_800.index), 1)
assert np.isclose(trapezoid(ped_400.values, x=ped_400.index), 1)
assert np.isclose(trapezoid(ped_800.values, x=ped_800.index), 1)


def test_beta_phi2a():
Expand All @@ -127,8 +132,8 @@ def test_beta_phi2a():

assert np.isclose((ped_400.iloc[166]), 0.00045, atol=1e-4, rtol=1e-2)
assert np.isclose((ped_800.iloc[166]), 0.0522, atol=1e-3, rtol=1e-2)
assert np.isclose(np.trapz(ped_400.values, x=ped_400.index), 1)
assert np.isclose(np.trapz(ped_800.values, x=ped_800.index), 1)
assert np.isclose(trapezoid(ped_400.values, x=ped_400.index), 1)
assert np.isclose(trapezoid(ped_800.values, x=ped_800.index), 1)


def test_beta_phi3a():
Expand All @@ -146,8 +151,8 @@ def test_beta_phi3a():

assert np.isclose((ped_400.iloc[166]), 0.001, atol=1e-3, rtol=1e-2)
assert np.isclose((ped_800.iloc[166]), 0.0589, atol=1e-3, rtol=1e-2)
assert np.isclose(np.trapz(ped_400.values, x=ped_400.index), 1)
assert np.isclose(np.trapz(ped_800.values, x=ped_800.index), 1)
assert np.isclose(trapezoid(ped_400.values, x=ped_400.index), 1)
assert np.isclose(trapezoid(ped_800.values, x=ped_800.index), 1)


def test_rovib_dos():
Expand All @@ -167,8 +172,8 @@ def test_rovib_dos():

assert np.isclose((ped_1800.iloc[166]), 0.01985, atol=1e-3, rtol=1e-2)
assert np.isclose((ped_2000.iloc[166]), 0.01724, atol=1e-3, rtol=1e-2)
assert np.isclose(np.trapz(ped_1800.values, x=ped_1800.index), 1)
assert np.isclose(np.trapz(ped_2000.values, x=ped_2000.index), 1)
assert np.isclose(trapezoid(ped_1800.values, x=ped_1800.index), 1)
assert np.isclose(trapezoid(ped_2000.values, x=ped_2000.index), 1)


def test_thermal():
Expand All @@ -187,8 +192,8 @@ def test_thermal():

assert np.isclose((ped_1800.iloc[166]), 0.01924792, atol=1e-5, rtol=1e-5)
assert np.isclose((ped_2000.iloc[166]), 0.01132143, atol=1e-5, rtol=1e-5)
assert np.isclose(np.trapz(ped_1800.values, x=ped_1800.index), 1)
assert np.isclose(np.trapz(ped_2000.values, x=ped_2000.index), 1)
assert np.isclose(trapezoid(ped_1800.values, x=ped_1800.index), 1)
assert np.isclose(trapezoid(ped_2000.values, x=ped_2000.index), 1)


def test_bf_from_phi1a():
Expand Down
Loading