NumPy 2.0 support (#255)

espei/core_utils.py

@@ -32,7 +32,7 @@ def filter_configurations(desired_data: List[Dataset], configuration, symmetry)
                                      for sblconf in data['solver']['sublattice_configurations']])
         matching_configs = np.arange(len(data['solver']['sublattice_configurations']))[matching_configs]
         # Rewrite output values with filtered data
-        data['values'] = np.array(data['values'], dtype=np.float_)[..., matching_configs]
+        data['values'] = np.array(data['values'], dtype=np.float64)[..., matching_configs]
         data['solver']['sublattice_configurations'] = recursive_tuplify(np.array(data['solver']['sublattice_configurations'], dtype=np.object_)[matching_configs].tolist())
         if 'sublattice_occupancies' in data['solver']:
             data['solver']['sublattice_occupancies'] = np.array(data['solver']['sublattice_occupancies'], dtype=np.object_)[matching_configs].tolist()
@@ -62,7 +62,7 @@ def filter_temperatures(desired_data: List[Dataset]) -> List[Dataset]:
     for data in desired_data:
         temp_filter = np.atleast_1d(data['conditions']['T']) >= 298.15
         data['conditions']['T'] = np.atleast_1d(data['conditions']['T'])[temp_filter]
-        data['values'] = np.array(data['values'], dtype=np.float_)[..., temp_filter, :].tolist()
+        data['values'] = np.array(data['values'], dtype=np.float64)[..., temp_filter, :].tolist()
     return desired_data
 
 

espei/datasets.py

@@ -183,7 +183,7 @@ def check_dataset(dataset: Dataset):
                 if len(component_list) != len(mole_fraction_list):
                     raise DatasetError('The length of the components list and mole fractions list in tieline {} for the ZPF point {} should be the same.'.format(tieline, zpf))
                 # check that all mole fractions are less than one
-                mf_sum = np.nansum(np.array(mole_fraction_list, dtype=np.float_))
+                mf_sum = np.nansum(np.array(mole_fraction_list, dtype=np.float64))
                 if any([mf is not None for mf in mole_fraction_list]) and mf_sum > 1.0:
                     raise DatasetError('Mole fractions for tieline {} for the ZPF point {} sum to greater than one.'.format(tieline, zpf))
 

espei/error_functions/equilibrium_thermochemical_error.py

@@ -211,8 +211,8 @@ def calc_prop_differences(eqpropdata: EqPropData,
     update_phase_record_parameters(phase_records, parameters)
     params_dict = OrderedDict(zip(map(str, eqpropdata.params_keys), parameters))
     output = eqpropdata.output
-    weights = np.array(eqpropdata.weight, dtype=np.float_)
-    samples = np.array(eqpropdata.samples, dtype=np.float_)
+    weights = np.array(eqpropdata.weight, dtype=np.float64)
+    samples = np.array(eqpropdata.samples, dtype=np.float64)
 
     calculated_data = []
     for comp_conds in eqpropdata.composition_conds:
@@ -230,7 +230,7 @@ def calc_prop_differences(eqpropdata: EqPropData,
         vals = getattr(propdata, output).flatten().tolist()
         calculated_data.extend(vals)
 
-    calculated_data = np.array(calculated_data, dtype=np.float_)
+    calculated_data = np.array(calculated_data, dtype=np.float64)
 
     assert calculated_data.shape == samples.shape, f"Calculated data shape {calculated_data.shape} does not match samples shape {samples.shape}"
     assert calculated_data.shape == weights.shape, f"Calculated data shape {calculated_data.shape} does not match weights shape {weights.shape}"

espei/error_functions/non_equilibrium_thermochemical_error.py

@@ -49,7 +49,7 @@ def filter_sublattice_configurations(desired_data: List[Dataset], subl_model) ->
         matching_configs = np.asarray(matching_configs, dtype=np.bool_)
 
         # Rewrite output values with filtered data
-        data['values'] = np.array(data['values'], dtype=np.float_)[..., matching_configs]
+        data['values'] = np.array(data['values'], dtype=np.float64)[..., matching_configs]
         data['solver']['sublattice_configurations'] = np.array(data['solver']['sublattice_configurations'], dtype=np.object_)[matching_configs].tolist()
         if 'sublattice_occupancies' in data['solver']:
             data['solver']['sublattice_occupancies'] = np.array(data['solver']['sublattice_occupancies'], dtype=np.object_)[matching_configs].tolist()

espei/error_functions/zpf_error.py

@@ -104,7 +104,7 @@ def _phase_is_stoichiometric(mod):
 def _compute_vertex_composition(comps: Sequence[str], comp_conds: Dict[str, float]):
     """Compute the overall composition in a vertex assuming an N=1 normalization condition"""
     pure_elements = sorted(c for c in comps if c != 'VA')
-    vertex_composition = np.empty(len(pure_elements), dtype=np.float_)
+    vertex_composition = np.empty(len(pure_elements), dtype=np.float64)
     unknown_indices = []
     for idx, el in enumerate(pure_elements):
         amt = comp_conds.get(v.X(el), None)
@@ -271,7 +271,7 @@ def estimate_hyperplane(phase_region: PhaseRegion, parameters: np.ndarray, appro
                 target_hyperplane_chempots.append(np.full_like(MU_values, np.nan))
             else:
                 target_hyperplane_chempots.append(MU_values)
-    target_hyperplane_mean_chempots = np.nanmean(target_hyperplane_chempots, axis=0, dtype=np.float_)
+    target_hyperplane_mean_chempots = np.nanmean(target_hyperplane_chempots, axis=0, dtype=np.float64)
     return target_hyperplane_mean_chempots
 
 
@@ -301,7 +301,7 @@ def driving_force_to_hyperplane(target_hyperplane_chempots: np.ndarray,
         # Compute residual driving force
         # TODO: Check that it actually makes sense to declare this phase 'disordered'
         num_dof = sum([len(subl) for subl in models[current_phase].constituents])
-        desired_sitefracs = np.ones(num_dof, dtype=np.float_)
+        desired_sitefracs = np.ones(num_dof, dtype=np.float64)
         dof_idx = 0
         for subl in models[current_phase].constituents:
             dof = sorted(subl, key=str)
@@ -313,7 +313,7 @@ def driving_force_to_hyperplane(target_hyperplane_chempots: np.ndarray,
                 else:
                     sitefracs_to_add = [1.0]
             else:
-                sitefracs_to_add = np.array([cond_dict.get(v.X(d)) for d in dof], dtype=np.float_)
+                sitefracs_to_add = np.array([cond_dict.get(v.X(d)) for d in dof], dtype=np.float64)
                 # Fix composition of dependent component
                 sitefracs_to_add[np.isnan(sitefracs_to_add)] = 1 - np.nansum(sitefracs_to_add)
             desired_sitefracs[dof_idx:dof_idx + num_subl_dof] = sitefracs_to_add

espei/optimizers/opt_mcmc.py

@@ -272,7 +272,7 @@ def predict(params, **ctx):
 
         # Important to coerce to floats here because the values _must_ be floats if
         # they are used to update PhaseRecords directly
-        params = np.asarray(params, dtype=np.float_)
+        params = np.asarray(params, dtype=np.float64)
 
         # lnprior
         prior_rvs = ctx.get('prior_rvs', [rv_zero() for _ in range(params.size)])

espei/parameter_selection/fitting_steps.py

@@ -152,7 +152,7 @@ def get_response_vector(cls, fixed_model: Model, fixed_portions: [symengine.Basi
         site_fractions = []
         for ds in data:
             for _ in ds['conditions']['T']:
-                sf = build_sitefractions(fixed_model.phase_name, ds['solver']['sublattice_configurations'], ds['solver'].get('sublattice_occupancies', np.ones((len(ds['solver']['sublattice_configurations']), len(ds['solver']['sublattice_configurations'][0])), dtype=np.float_)))
+                sf = build_sitefractions(fixed_model.phase_name, ds['solver']['sublattice_configurations'], ds['solver'].get('sublattice_occupancies', np.ones((len(ds['solver']['sublattice_configurations']), len(ds['solver']['sublattice_configurations'][0])), dtype=np.float64)))
                 site_fractions.append(sf)
         site_fractions = list(itertools.chain(*site_fractions))
         # If any site fractions show up in our rhs that aren't in these
@@ -166,7 +166,7 @@ def get_response_vector(cls, fixed_model: Model, fixed_portions: [symengine.Basi
         # also replace with database symbols in case we did higher order fitting
         rhs = [fixed_model.symbol_replace(symengine.S(value_with_symbols).xreplace(sf), fixed_model._symbols).evalf() for value_with_symbols, sf in zip(rhs, site_fractions)]
         # cast to float, confirming that these are concrete values with no sybolics
-        rhs = np.asarray(rhs, dtype=np.float_)
+        rhs = np.asarray(rhs, dtype=np.float64)
         return rhs
 
 
@@ -248,7 +248,7 @@ def get_response_vector(cls, fixed_model: Model, fixed_portions: [symengine.Basi
         site_fractions = []
         for ds in data:
             for _ in ds['conditions']['T']:
-                sf = build_sitefractions(phase_name, ds['solver']['sublattice_configurations'], ds['solver'].get('sublattice_occupancies', np.ones((len(ds['solver']['sublattice_configurations']), len(ds['solver']['sublattice_configurations'][0])), dtype=np.float_)))
+                sf = build_sitefractions(phase_name, ds['solver']['sublattice_configurations'], ds['solver'].get('sublattice_occupancies', np.ones((len(ds['solver']['sublattice_configurations']), len(ds['solver']['sublattice_configurations'][0])), dtype=np.float64)))
                 site_fractions.append(sf)
         site_fractions = list(itertools.chain(*site_fractions))
 
@@ -267,7 +267,7 @@ def get_response_vector(cls, fixed_model: Model, fixed_portions: [symengine.Basi
             sf.update(cond_dict)
         # also replace with database symbols in case we did higher order fitting
         data_qtys = [fixed_model.symbol_replace(symengine.S(i).xreplace(sf), fixed_model._symbols).evalf() for i, sf in zip(data_qtys, site_fractions)]
-        data_qtys = np.asarray(data_qtys, dtype=np.float_)
+        data_qtys = np.asarray(data_qtys, dtype=np.float64)
         return data_qtys
 
 

espei/parameter_selection/selection.py

@@ -78,7 +78,7 @@ def score_model(feature_matrix, data_quantities, model_coefficients, feature_lis
     """
     p = aicc_factor if aicc_factor is not None else 1.0
     k = len(feature_list)
-    rss = np.square((np.dot(feature_matrix, model_coefficients) - data_quantities.astype(np.float_))*np.array(weights)).sum()
+    rss = np.square((np.dot(feature_matrix, model_coefficients) - data_quantities.astype(np.float64))*np.array(weights)).sum()
     if np.abs(rss) < rss_numerical_limit:
         rss = rss_numerical_limit
     # Compute the corrected Akaike Information Criterion

espei/plot.py

@@ -597,7 +597,7 @@ def plot_interaction(dbf, comps, phase_name, configuration, output, datasets=Non
     symbol_map = bib_marker_map(bib_reference_keys)
     for data in desired_data:
         indep_var_data = None
-        response_data = np.zeros_like(data['values'], dtype=np.float_)
+        response_data = np.zeros_like(data['values'], dtype=np.float64)
         if disordered_config:
             # Take the second element of the first interacting sublattice as the coordinate
             # Because it's disordered all sublattices should be equivalent
@@ -625,15 +625,15 @@ def plot_interaction(dbf, comps, phase_name, configuration, output, datasets=Non
                 points[point_idx].update({key: 0 for key in missing_variables})
                 # Change entry to a sorted array of site fractions
                 points[point_idx] = list(OrderedDict(sorted(points[point_idx].items(), key=str)).values())
-            points = np.array(points, dtype=np.float_)
+            points = np.array(points, dtype=np.float64)
             # TODO: Real temperature support
             points = points[None, None]
             stability = calculate(dbf, comps, [phase_name], output=data['output'][:-5],
                                     T=temps, P=pressures, points=points,
                                     model=mod_srf)
             response_data -= stability[data['output'][:-5]].values.squeeze()
 
-        response_data += np.array(data['values'], dtype=np.float_)
+        response_data += np.array(data['values'], dtype=np.float64)
         response_data = response_data.flatten()
         ref = data.get('reference', '')
         dataplot_kwargs.setdefault('markersize', 8)
@@ -734,10 +734,10 @@ def plot_endmember(dbf, comps, phase_name, configuration, output, datasets=None,
     symbol_map = bib_marker_map(bib_reference_keys)
     for data in desired_data:
         indep_var_data = None
-        response_data = np.zeros_like(data['values'], dtype=np.float_)
-        indep_var_data = np.array(data['conditions'][x], dtype=np.float_).flatten()
+        response_data = np.zeros_like(data['values'], dtype=np.float64)
+        indep_var_data = np.array(data['conditions'][x], dtype=np.float64).flatten()
 
-        response_data += np.array(data['values'], dtype=np.float_)
+        response_data += np.array(data['values'], dtype=np.float64)
         response_data = response_data.flatten()
         ref = data.get('reference', '')
         dataplot_kwargs.setdefault('markersize', 8)
@@ -875,9 +875,9 @@ def _compare_data_to_parameters(dbf, comps, phase_name, desired_data, mod, confi
 
     for data in desired_data:
         indep_var_data = None
-        response_data = np.zeros_like(data['values'], dtype=np.float_)
+        response_data = np.zeros_like(data['values'], dtype=np.float64)
         if x == 'T' or x == 'P':
-            indep_var_data = np.array(data['conditions'][x], dtype=np.float_).flatten()
+            indep_var_data = np.array(data['conditions'][x], dtype=np.float64).flatten()
         elif x == 'Z':
             if disordered_config:
                 # Take the second element of the first interacting sublattice as the coordinate
@@ -909,15 +909,15 @@ def _compare_data_to_parameters(dbf, comps, phase_name, desired_data, mod, confi
                     points[point_idx].update({key: 0 for key in missing_variables})
                     # Change entry to a sorted array of site fractions
                     points[point_idx] = list(OrderedDict(sorted(points[point_idx].items(), key=str)).values())
-                points = np.array(points, dtype=np.float_)
+                points = np.array(points, dtype=np.float64)
                 # TODO: Real temperature support
                 points = points[None, None]
                 stability = calculate(dbf, comps, [phase_name], output=data['output'][:-5],
                                       T=temps, P=pressures, points=points,
                                       model=mod_latticeonly, mode='numpy')
                 response_data -= stability[data['output'][:-5]].values.squeeze()
 
-        response_data += np.array(data['values'], dtype=np.float_)
+        response_data += np.array(data['values'], dtype=np.float64)
         response_data = response_data.flatten()
         if not bar_chart:
             extra_kwargs = {}

espei/priors.py

@@ -18,11 +18,11 @@ class rv_zero(object):
     --------
     >>> import numpy as np
     >>> rv = rv_zero()
-    >>> np.isclose(rv.logpdf(-np.inf), 0.0)
+    >>> bool(np.isclose(rv.logpdf(-np.inf), 0.0))
     True
-    >>> np.isclose(rv.logpdf(1.0), 0.0)
+    >>> bool(np.isclose(rv.logpdf(1.0), 0.0))
     True
-    >>> np.isclose(rv.logpdf(0.0), 0.0)
+    >>> bool(np.isclose(rv.logpdf(0.0), 0.0))
     True
 
     """
@@ -41,16 +41,16 @@ class DistributionParameter(object):
     Examples
     --------
     >>> dp = DistributionParameter(5.0, 'absolute')  # always get back 5
-    >>> dp.value(1.0) == 5.0
+    >>> bool(dp.value(1.0) == 5.0)
     True
     >>> dp = DistributionParameter(-2.0, 'relative')  # multiply by -2
-    >>> dp.value(2.0) == -4.0
+    >>> bool(dp.value(2.0) == -4.0)
     True
     >>> dp = DistributionParameter(-1.0, 'shift_absolute')  # subtract 1
-    >>> dp.value(2.0) == 1.0
+    >>> bool(dp.value(2.0) == 1.0)
     True
     >>> dp = DistributionParameter(-0.5, 'shift_relative')  # subtract 1/2 value
-    >>> dp.value(2.0) == 1.0
+    >>> bool(dp.value(2.0) == 1.0)
     True
 
     """
@@ -155,9 +155,9 @@ def get_prior(self, value):
         >>> import numpy as np
         >>> from espei.priors import PriorSpec
         >>> tri_spec = {'name': 'triangular', 'loc_shift_relative': -0.5, 'scale_shift_relative': 0.5, 'c': 0.5}
-        >>> np.isneginf(PriorSpec(**tri_spec).get_prior(10).logpdf(5.1))
+        >>> bool(np.isneginf(PriorSpec(**tri_spec).get_prior(10).logpdf(5.1)))
         False
-        >>> np.isneginf(PriorSpec(**tri_spec).get_prior(10).logpdf(4.9))
+        >>> bool(np.isneginf(PriorSpec(**tri_spec).get_prior(10).logpdf(4.9)))
         True
 
         """

espei/shadow_functions.py

@@ -25,7 +25,7 @@ def update_phase_record_parameters(phase_records: Dict[str, PhaseRecord], parame
             # very important that these are floats, otherwise parameters can end up
             # with garbage data. `np.asarray` does not create a copy if the type is
             # correct
-            phase_record.parameters[:] = np.asarray(parameters, dtype=np.float_)
+            phase_record.parameters[:] = np.asarray(parameters, dtype=np.float64)
 
 def _single_phase_start_point(conditions, state_variables, phase_records, grid):
     """Return a single CompositionSet object to use in a point calculation

espei/utils.py

@@ -100,7 +100,7 @@ def optimal_parameters(trace_array, lnprob_array, kth=0):
     >>> from espei.utils import optimal_parameters
     >>> trace = np.array([[[1, 0], [2, 0], [3, 0], [0, 0]], [[0, 2], [0, 4], [0, 6], [0, 0]]])  # 3 iterations of 4 allocated
     >>> lnprob = np.array([[-6, -4, -2, 0], [-3, -1, -2, 0]])
-    >>> np.all(np.isclose(optimal_parameters(trace, lnprob), np.array([0, 4])))
+    >>> bool(np.all(np.isclose(optimal_parameters(trace, lnprob), np.array([0, 4]))))
     True
 
     """

tests/test_error_functions.py

@@ -476,7 +476,7 @@ def test_equilibrium_thermochemical_error_computes_correct_probability(datasets_
     errors, weights = calc_prop_differences(eqdata[0], np.array([-31626.6]))
     assert np.all(np.isclose(errors, [-31626.6*0.5*0.5]))
     # change to -40000
-    errors, weights = calc_prop_differences(eqdata[0], np.array([-40000], np.float_))
+    errors, weights = calc_prop_differences(eqdata[0], np.array([-40000], np.float64))
     assert np.all(np.isclose(errors, [-40000*0.5*0.5]))
 
 

tests/test_mcmc.py

@@ -164,7 +164,7 @@ def test_equilibrium_thermochemical_correct_probability(datasets_db):
     assert np.isclose(prob, expected_prob)
 
     # change to -40000
-    prob = opt.predict(np.array([-40000], dtype=np.float_), **ctx)
+    prob = opt.predict(np.array([-40000], dtype=np.float64), **ctx)
     expected_prob = norm(loc=0, scale=500).logpdf([-40000*0.5*0.5]).sum()
     assert np.isclose(prob, expected_prob)