diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 00000000..f8dcb6ab
Binary files /dev/null and b/.DS_Store differ
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 76eba228..0b6b644f 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,7 @@
 Changelog
 =========
 
+- Enable `numpy >= 2.0` 
 - Add RobustGPyRegression model that can handle non-finite output values
 - Update surrogate model initialisation to use all initial evidence
 - Update BOLFI and BOLFIRE to use a shared sample class that returns individual chains in the arviz inference data
diff --git a/docs/.DS_Store b/docs/.DS_Store
new file mode 100644
index 00000000..f1ff1ed9
Binary files /dev/null and b/docs/.DS_Store differ
diff --git a/elfi/examples/bdm.py b/elfi/examples/bdm.py
index 15d8948d..d303a414 100644
--- a/elfi/examples/bdm.py
+++ b/elfi/examples/bdm.py
@@ -27,7 +27,7 @@ def prepare_inputs(*inputs, **kwinputs):
     meta = kwinputs['meta']
     # Organize the parameters to an array. The broadcasting works nicely with constant
     # arguments.
-    param_array = np.row_stack(
+    param_array = np.vstack(
         [(a, d, t, N) for (a, d, t, N) in np.broadcast(alpha, delta, tau, N)]
         )
 
diff --git a/elfi/examples/bignk.py b/elfi/examples/bignk.py
index 4e67f38e..e169b7ca 100644
--- a/elfi/examples/bignk.py
+++ b/elfi/examples/bignk.py
@@ -83,7 +83,7 @@ def BiGNK(A1, A2, B1, B2, g1, g2, k1, k2, rho, c=.8, n_obs=150, batch_size=1, ra
     z_batches = []
     for i in range(batch_size):
         # Initialising a separate covariance matrix for each batch.
-        matrix_cov = np.array([[1, float(rho[i])], [float(rho[i]), 1]])
+        matrix_cov = np.array([[1, rho[i, 0]], [rho[i, 0], 1]])
 
         z_batch = ss.multivariate_normal.rvs(cov=matrix_cov, size=n_obs, random_state=random_state)
         z_batches.append(z_batch)
diff --git a/elfi/methods/bo/gpy_regression.py b/elfi/methods/bo/gpy_regression.py
index d921f811..e42dbf66 100644
--- a/elfi/methods/bo/gpy_regression.py
+++ b/elfi/methods/bo/gpy_regression.py
@@ -149,8 +149,8 @@ def predict(self, x, noiseless=False):
     # TODO: find a more general solution
     # cache some RBF-kernel-specific values for faster sampling
     def _cache_RBF_kernel(self):
-        self._rbf_var = float(self._gp.kern.rbf.variance)
-        self._rbf_factor = -0.5 / float(self._gp.kern.rbf.lengthscale)**2
+        self._rbf_var = self._gp.kern.rbf.variance.item()
+        self._rbf_factor = -0.5 / self._gp.kern.rbf.lengthscale.item()**2
         self._rbf_bias = float(self._gp.kern.bias.K(self._gp.X)[0, 0])
         self._rbf_noisevar = float(self._gp.likelihood.variance[0])
         self._rbf_woodbury = self._gp.posterior.woodbury_vector
@@ -319,7 +319,7 @@ def optimize(self):
         logger.debug("Optimizing GP hyperparameters")
         try:
             self._gp.optimize(self.optimizer, max_iters=self.max_opt_iters)
-        except np.linalg.linalg.LinAlgError:
+        except np.linalg.LinAlgError:
             logger.warning("Numerical error in GP optimization. Stopping optimization")
 
     @property
diff --git a/elfi/methods/bo/utils.py b/elfi/methods/bo/utils.py
index 3e94043e..29b0f7bc 100644
--- a/elfi/methods/bo/utils.py
+++ b/elfi/methods/bo/utils.py
@@ -94,11 +94,15 @@ def minimize(fun,
     # Run the optimisation from each initialization point.
     locs = []
     vals = np.empty(n_start_points)
+    options = {'maxiter': maxiter}
+    if method == 'SLSQP':
+        options['ftol'] = 1e-10
+
     for i in range(n_start_points):
         result = scipy.optimize.minimize(fun, start_points[i, :],
                                          method=method, jac=grad,
                                          bounds=bounds, constraints=constraints,
-                                         options={'maxiter': maxiter})
+                                         options=options)
         locs.append(result['x'])
         vals[i] = result['fun']
 
diff --git a/elfi/methods/bsl/pdf_methods.py b/elfi/methods/bsl/pdf_methods.py
index 680733cc..5c1008c4 100644
--- a/elfi/methods/bsl/pdf_methods.py
+++ b/elfi/methods/bsl/pdf_methods.py
@@ -224,9 +224,9 @@ def semi_param_kernel_estimate(ssx, ssy, shrinkage=None, penalty=None, whitening
         # NOTE: bw_method - "silverman" is being used here is slightly
         #       different than "nrd0" - silverman's rule of thumb in R.
         kernel = ss.gaussian_kde(ssx_j, bw_method="silverman")
-        logpdf_y[j] = kernel.logpdf(y)
+        logpdf_y[j] = np.asarray(kernel.logpdf(y)).item()
 
-        y_u[j] = kernel.integrate_box_1d(np.NINF, y)
+        y_u[j] = kernel.integrate_box_1d(-np.inf, y)
         y_u[j] = min(1, y_u[j])  # fix numerical errors, CDF values cannot exceed 1
 
         if whitening is not None:
diff --git a/elfi/methods/bsl/pre_sample_methods.py b/elfi/methods/bsl/pre_sample_methods.py
index fbf37a0a..09e43976 100644
--- a/elfi/methods/bsl/pre_sample_methods.py
+++ b/elfi/methods/bsl/pre_sample_methods.py
@@ -300,8 +300,9 @@ def select_penalty(model, n_sim, theta, feature_names, likelihood=None,
                                             whitening=whitening)
                     except FloatingPointError as err:
                         logger.warning('Floating point error: {}'.format(err))
-                        loglik = np.NINF
-                    logliks[m_iteration, n_iteration, lmda_iteration] = loglik
+                        loglik = -np.inf
+                    logliks[m_iteration, n_iteration, lmda_iteration] = \
+                        np.asarray(loglik).item()
 
     # choose the lambda with the empirical s.d. of the log SL estimates
     # closest to sigma
diff --git a/elfi/methods/density_ratio_estimation.py b/elfi/methods/density_ratio_estimation.py
index 479a03c9..8c0dbef4 100644
--- a/elfi/methods/density_ratio_estimation.py
+++ b/elfi/methods/density_ratio_estimation.py
@@ -162,7 +162,7 @@ def _KLIEP_lcv(self, x, y, sigma):
         non_null = np.any(A > 1e-64, axis=1)
         non_null_length = sum(non_null)
         if non_null_length == 0:
-            return np.Inf
+            return np.inf
 
         A_full = A[non_null, :]
         x_full = x[non_null, :]
diff --git a/elfi/methods/inference/bsl.py b/elfi/methods/inference/bsl.py
index d8514053..a69e1f7e 100644
--- a/elfi/methods/inference/bsl.py
+++ b/elfi/methods/inference/bsl.py
@@ -155,7 +155,7 @@ def _init_state(self, n_samples, params0=None, gamma0=None):
         self.state['params'] = np.zeros((n_samples, len(self.parameter_names)))
         self.state['params'][0] = params0
         self.state['logprior'] = np.zeros((n_samples))
-        self.state['logprior'][0] = self.prior.logpdf(params0)
+        self.state['logprior'][0] = np.asarray(self.prior.logpdf(params0)).item()
         self.state['logposterior'] = np.zeros((n_samples))
         if self.is_misspec:
             self.state['gamma'] = np.zeros((n_samples, self.observed.size))
@@ -214,7 +214,7 @@ def _init_round(self):
                 self.state['logposterior'][n-1] = ll + self.state['logprior'][n-1]
             # sample candidate parameter values
             prop = self._propagate_state()
-            logprior = self.prior.logpdf(prop)
+            logprior = np.asarray(self.prior.logpdf(prop)).item()
             if np.isfinite(logprior):
                 # start data collection with the proposed parameter values
                 self.state['logprior'][n] = logprior
@@ -240,7 +240,7 @@ def _process_simulated(self):
         """
         # estimate synthetic likelihood
         if not np.all(np.isfinite(self.simulated)):
-            loglikelihood = np.NINF
+            loglikelihood = -np.inf
         else:
             if self.is_misspec:
                 gamma = self.gamma_sampler_state['gamma']
@@ -248,6 +248,7 @@ def _process_simulated(self):
             else:
                 loglikelihood = self.likelihood(self.simulated, self.observed)
 
+        loglikelihood = np.asarray(loglikelihood).item()
         n = self.state['n_samples']
         if not np.isfinite(loglikelihood):
             if n == 0:
diff --git a/elfi/methods/inference/romc.py b/elfi/methods/inference/romc.py
index 5cd21503..1418d650 100644
--- a/elfi/methods/inference/romc.py
+++ b/elfi/methods/inference/romc.py
@@ -571,7 +571,7 @@ def _det_generator(self, theta, seed):
         param_dict = flat_array_to_dict(model.parameter_names, theta)
         dict_outputs = model.generate(
             batch_size=1, outputs=[output_node], with_values=param_dict, seed=int(seed))
-        return float(dict_outputs[output_node]) ** 2
+        return np.asarray(dict_outputs[output_node]).item() ** 2
 
     def _freeze_seed(self, seed):
         """Freeze the model.generate with a specific seed.
@@ -1565,7 +1565,7 @@ def fit_local_surrogate(self, **kwargs) -> None:
         def local_surrogate(theta, model_scikit):
             assert theta.ndim == 1
             theta = np.expand_dims(theta, 0)
-            return float(model_scikit.predict(theta))
+            return np.asarray(model_scikit.predict(theta)).item()
 
         def create_local_surrogate(model):
             return partial(local_surrogate, model_scikit=model)
@@ -1899,8 +1899,12 @@ def _find_rotation_vector(self, hess_appr: np.ndarray) -> np.ndarray:
             hess_appr = np.eye(dim)
         eig_val, eig_vec = np.linalg.eig(hess_appr)
 
+        if np.iscomplexobj(eig_vec) and np.allclose(eig_vec.imag, 0):
+            eig_vec = eig_vec.real
+
         # if extreme values appear, return the I matrix
-        if np.isnan(np.sum(eig_vec)) or np.isinf(np.sum(eig_vec)) or (eig_vec.dtype == complex):
+        if (np.isnan(np.sum(eig_vec)) or np.isinf(np.sum(eig_vec))
+                or np.iscomplexobj(eig_vec)):
             logger.info("Eye matrix return as rotation.")
             eig_vec = np.eye(dim)
         if np.linalg.matrix_rank(eig_vec) < dim:
diff --git a/elfi/methods/inference/samplers.py b/elfi/methods/inference/samplers.py
index e09ba267..17e6e6d4 100644
--- a/elfi/methods/inference/samplers.py
+++ b/elfi/methods/inference/samplers.py
@@ -116,7 +116,7 @@ def set_objective(self, n_samples, threshold=None, quantile=None, n_sim=None):
         """
         if quantile is None and threshold is None and n_sim is None:
             quantile = .01
-        self.state = dict(samples=None, threshold=np.Inf,
+        self.state = dict(samples=None, threshold=np.inf,
                           n_sim=0, accept_rate=1, n_batches=0)
 
         if quantile:
diff --git a/elfi/methods/mcmc.py b/elfi/methods/mcmc.py
index f6912c24..5c8a30c3 100644
--- a/elfi/methods/mcmc.py
+++ b/elfi/methods/mcmc.py
@@ -326,7 +326,7 @@ def _build_tree_nuts(params, momentum, log_slicevar, step, depth, log_joint0, ta
         momentum1 = momentum1 + 0.5 * step * grad_target(params1)
 
         log_joint = target(params1) - 0.5 * np.inner(momentum1, momentum1)
-        n_ok = float(log_slicevar <= log_joint)
+        n_ok = float(np.asarray(log_slicevar <= log_joint).item())
         sub_ok = log_slicevar < (1000. + log_joint)  # check for diverging error
         is_out = False
         if not sub_ok:
@@ -338,7 +338,7 @@ def _build_tree_nuts(params, momentum, log_slicevar, step, depth, log_joint0, ta
                     "momentum1={}.".format(log_joint, params, params1, momentum, momentum1))
             mh_ratio = 0.  # reject
         else:
-            mh_ratio = min(1., np.exp(log_joint - log_joint0))
+            mh_ratio = np.minimum(1., np.exp(log_joint - log_joint0)).item()
 
         return params1, momentum1, params1, momentum1, params1, n_ok, sub_ok, mh_ratio, 1., \
             not sub_ok, is_out
diff --git a/elfi/methods/posteriors.py b/elfi/methods/posteriors.py
index 9de41ba4..698086b0 100644
--- a/elfi/methods/posteriors.py
+++ b/elfi/methods/posteriors.py
@@ -485,7 +485,7 @@ def _pdf_unnorm_single_point(self, theta: np.ndarray) -> float:
         assert theta.ndim == 1
 
         prior = self.prior
-        pr = float(prior.pdf(np.expand_dims(theta, 0)))
+        pr = np.asarray(prior.pdf(np.expand_dims(theta, 0))).item()
 
         if self.surrogate_used:
             indicator_sum = self._sum_over_regions_indicators(theta)
@@ -570,7 +570,7 @@ def _worker_compute_weight(self, args):
             q = region.pdf(cur_theta)
             if q == 0.0:
                 logger.warning("Zero q")
-            pr = float(prior.pdf(np.expand_dims(cur_theta, 0)))
+            pr = np.asarray(prior.pdf(np.expand_dims(cur_theta, 0))).item()
             dist = func(cur_theta)
             distances.append(dist)
             ind = dist < eps
@@ -742,7 +742,7 @@ def sample(self, n2: int, seed=None) -> (np.ndarray, np.ndarray):
                     if q == 0.0:
                         logger.warning("Zero q")
                     # (ii) p
-                    pr = float(prior.pdf(np.expand_dims(cur_theta, 0)))
+                    pr = np.asarray(prior.pdf(np.expand_dims(cur_theta, 0))).item()
 
                     # (iii) indicator
                     dist = funcs[i](cur_theta)
diff --git a/tests/unit/test_failure_robust_model.py b/tests/unit/test_failure_robust_model.py
index c0b4b47c..222b0d66 100644
--- a/tests/unit/test_failure_robust_model.py
+++ b/tests/unit/test_failure_robust_model.py
@@ -116,7 +116,7 @@ def test_success_proba(model_2):
 def test_success_proba_default(model_1):
     target_model = model_1
     prob = target_model.success_proba([3, 3])
-    assert float(prob) == 1
+    assert np.array_equal(prob, np.ones((1, 1)))
 
 
 def test_success_proba_gradients(model_2):
diff --git a/tests/unit/test_testbench.py b/tests/unit/test_testbench.py
index adf191d6..49d2f2e9 100644
--- a/tests/unit/test_testbench.py
+++ b/tests/unit/test_testbench.py
@@ -1,5 +1,4 @@
 import numpy as np
-from numpy.lib.function_base import quantile
 import pytest
 
 import elfi
diff --git a/tests/unit/test_utils.py b/tests/unit/test_utils.py
index 565bf8bf..375fd09d 100644
--- a/tests/unit/test_utils.py
+++ b/tests/unit/test_utils.py
@@ -203,8 +203,8 @@ def test_sample_object_to_dict():
 
 
 def test_numpy_to_python_type():
-    data = dict(a=np.array([1, 2, 3, 4]), b=np.uint(5), c=np.float_(10),
-                d=dict(a=np.array([0, 9, 8, 7]), b=np.uint(15), c=np.float_(12)))
+    data = dict(a=np.array([1, 2, 3, 4]), b=np.uint(5), c=np.float64(10),
+                d=dict(a=np.array([0, 9, 8, 7]), b=np.uint(15), c=np.float64(12)))
     numpy_to_python_type(data)
 
     # checking that our objects are jsonable is enough to be sure that numpy_to_python_type