diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index e1d8b90..40401d5 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -9,7 +9,7 @@ jobs:
     strategy:
       max-parallel: 3
       matrix:
-        python-version: [3.9, '3.10']
+        python-version: [3.9, '3.12']
 
     steps:
     - uses: actions/checkout@v1
diff --git a/examples/example_pineappl.py b/examples/example_pineappl.py
index cf43544..344bc17 100644
--- a/examples/example_pineappl.py
+++ b/examples/example_pineappl.py
@@ -61,7 +61,7 @@ def hadronic_pspgen(xarr, mmin, mmax):
     x2 = tau / y
     s = tau * smax
 
-    jacobian = np.math.log(tau0) * np.math.log(tau0) * tau * r1
+    jacobian = np.log(tau0) * np.log(tau0) * tau * r1
 
     # theta integration (in the CMS)
     cos_theta = 2.0 * r3 - 1.0
@@ -71,7 +71,7 @@ def hadronic_pspgen(xarr, mmin, mmax):
     u = -0.5 * s * (1.0 + cos_theta)
 
     # phi integration
-    jacobian *= 2.0 * np.math.acos(-1.0)
+    jacobian *= 2.0 * np.acos(-1.0)
 
     return s, t, u, x1, x2, jacobian
 
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..3caa792
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,58 @@
+[build-system]
+requires = ["setuptools >= 61.0"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "vegasflow"
+version = "1.4.0"
+description="Hardware agnostic Monte Carlo integration"
+authors = [
+	{name = "S. Carrazza", email = "stefano.carraza@cern.ch"},
+	{name = "J.Cruz-Martinez", email = "juacrumar@lairen.eu"}
+]
+classifiers=[
+	'Operating System :: Unix',
+	'Programming Language :: Python',
+	'Programming Language :: Python :: 3',
+	'Topic :: Scientific/Engineering',
+	'Topic :: Scientific/Engineering :: Physics',
+]
+dependencies = [
+	"joblib",
+	"numpy",
+	"tensorflow>2.2"
+]
+
+[project.optional-dependencies] 
+docs = [
+'sphinx_rtd_theme',
+'recommonmark',
+'sphinxcontrib-bibtex',
+]
+examples = [
+'cffi',
+'pineappl',
+'pdfflow',
+'scipy'
+]
+benchmark = [
+'vegas', # Lepage's Vegas for benchmarking
+]
+distribute = [
+'dask',
+'distributed',
+'dask-jobqueue',
+]
+
+[tool.black]
+line-length = 100
+skip_magic_trailing_comma = true
+
+[tool.isort]
+atomic = true
+line_length = 120
+profile = "black" # https://black.readthedocs.io/en/stable/guides/using_black_with_other_tools.html#custom-configuration
+skip_gitignore = true
+known_first_party = ["bgtrees"]
+force_sort_within_sections = true
+
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 9e09021..0000000
--- a/setup.py
+++ /dev/null
@@ -1,67 +0,0 @@
-# Installation script for python
-from setuptools import setup, find_packages
-import os
-import re
-
-
-requirements = ['joblib', 'numpy', 'tensorflow']
-PACKAGE = 'vegasflow'
-
-def get_version():
-    """ Gets the version from the package's __init__ file
-    if there is some problem, let it happily fail """
-    VERSIONFILE = os.path.join('src', PACKAGE, '__init__.py')
-    initfile_lines = open(VERSIONFILE, 'rt').readlines()
-    VSRE = r"^__version__ = ['\"]([^'\"]*)['\"]"
-    for line in initfile_lines:
-        mo = re.search(VSRE, line, re.M)
-        if mo:
-            return mo.group(1)
-
-this_directory = os.path.abspath(os.path.dirname(__file__))
-with open(os.path.join(this_directory, "README.md"), encoding="utf-8") as f:
-    long_description = f.read()
-
-setup(name='vegasflow',
-      version=get_version(),
-      description='Monte Carlo integration with Tensorflow',
-      author = 'S.Carrazza, J.Cruz-Martinez',
-      author_email='stefano.carrazza@cern.ch, juan.cruz@mi.infn.it',
-      url='https://github.com/N3PDF/VegasFlow',
-      package_dir={'':'src'},
-      packages=find_packages('src'),
-      zip_safe=False,
-      classifiers=[
-          'Operating System :: Unix',
-          'Programming Language :: Python',
-          'Programming Language :: Python :: 3',
-          'Topic :: Scientific/Engineering',
-          'Topic :: Scientific/Engineering :: Physics',
-      ],
-      install_requires=requirements,
-      extras_require={
-          'docs' : [
-            'sphinx_rtd_theme',
-            'recommonmark',
-            'sphinxcontrib-bibtex',
-            ],
-          'examples' : [
-            'cffi',
-            'pineappl',
-            'pdfflow',
-            'scipy'
-            ],
-          'benchmark' : [
-            'vegas', # Lepage's Vegas for benchmarking
-            ],
-          'distribute' : [
-            'dask',
-            'distributed',
-            'dask-jobqueue',
-            'tensorflow>2.2', # dask needs a more-pickable tensorflow
-            ],
-          },
-      python_requires='>=3.9',
-      long_description=long_description,
-      long_description_content_type="text/markdown",
-)
diff --git a/src/vegasflow/__init__.py b/src/vegasflow/__init__.py
index 5cad0fd..84d6129 100644
--- a/src/vegasflow/__init__.py
+++ b/src/vegasflow/__init__.py
@@ -1,10 +1,10 @@
 """Monte Carlo integration with Tensorflow"""
 
-from vegasflow.configflow import int_me, float_me, run_eager, DTYPE, DTYPEINT
+from vegasflow.configflow import DTYPE, DTYPEINT, float_me, int_me, run_eager
+from vegasflow.plain import PlainFlow, plain_sampler, plain_wrapper
 
 # Expose the main interfaces
-from vegasflow.vflow import VegasFlow, vegas_wrapper, vegas_sampler
-from vegasflow.plain import PlainFlow, plain_wrapper, plain_sampler
-from vegasflow.vflowplus import VegasFlowPlus, vegasflowplus_wrapper, vegasflowplus_sampler
+from vegasflow.vflow import VegasFlow, vegas_sampler, vegas_wrapper
+from vegasflow.vflowplus import VegasFlowPlus, vegasflowplus_sampler, vegasflowplus_wrapper
 
-__version__ = "1.3.0"
+__version__ = "1.4.0"
diff --git a/src/vegasflow/configflow.py b/src/vegasflow/configflow.py
index 4887d3f..06f2cf2 100644
--- a/src/vegasflow/configflow.py
+++ b/src/vegasflow/configflow.py
@@ -1,9 +1,12 @@
 """
 Define some constants, header style
 """
+
+import logging
+
 # Most of this can be moved to a yaml file without loss of generality
 import os
-import logging
+
 import numpy as np
 
 # Some global parameters
@@ -94,6 +97,8 @@ def run_eager(flag=True):
 
 
 FMAX = tf.constant(np.finfo(np.float64).max, dtype=DTYPE)
+
+
 # The wrappers below transform tensors and array to the correct type
 def int_me(i):
     """Cast the input to the `DTYPEINT` type"""
diff --git a/src/vegasflow/monte_carlo.py b/src/vegasflow/monte_carlo.py
index 029c8f3..482f9ba 100644
--- a/src/vegasflow/monte_carlo.py
+++ b/src/vegasflow/monte_carlo.py
@@ -32,28 +32,29 @@
     as it is the first one found idle.
 """
 
-import inspect
-import time
+from abc import ABC, abstractmethod
 import copy
-import threading
+import inspect
 import logging
-from abc import abstractmethod, ABC
+import threading
+import time
+
 import joblib
 import numpy as np
 import tensorflow as tf
+
 from vegasflow.configflow import (
-    MAX_EVENTS_LIMIT,
     DEFAULT_ACTIVE_DEVICES,
     DTYPE,
     DTYPEINT,
+    MAX_EVENTS_LIMIT,
     TECH_CUT,
     float_me,
-    int_me,
     fone,
     fzero,
+    int_me,
 )
 
-
 logger = logging.getLogger(__name__)
 
 
@@ -304,7 +305,8 @@ def _apply_integration_limits(self, rand):
 
     def _can_run_vectorial(self, expected_shape=None):
         """Accepting vectorial integrands depends on the algorithm,
-        if an algorithm can run on vectorial algorithms it should implement this method and return True"""
+        if an algorithm can run on vectorial algorithms it should implement this method and return True
+        """
         return self._CAN_RUN_VECTORIAL
 
     #### Integration management
@@ -561,7 +563,7 @@ def compile(self, integrand, compilable=True, signature=None, trace=False, check
                     signature = False
 
             compile_options = {
-                "experimental_autograph_options": tf.autograph.experimental.Feature.ALL,
+                "experimental_autograph_options": tf.autograph.experimental.Feature.ALL
             }
             if signature is None:
                 signature = autodiscover_signature
@@ -606,14 +608,14 @@ def batch_events(**kwargs):
                 )
             test_array = tf.random.uniform((event_size, self.n_dim), dtype=DTYPE)
             wgt = tf.random.uniform((event_size,), dtype=DTYPE)
-            res_tmp = new_integrand(test_array, weight=wgt).numpy()
+            res_tmp = new_integrand(test_array, weight=wgt)  # .numpy()
             res_shape = res_tmp.shape
 
             expected_shape = (event_size,)
 
             if len(res_shape) == 2:
                 self._vectorial = True
-                expected_shape = res_tmp.reshape(event_size, -1).shape
+                expected_shape = tf.reshape(res_tmp, (event_size, -1)).shape
                 if not self._can_run_vectorial(expected_shape):
                     raise NotImplementedError(
                         f"""The {self.__class__.__name__} algorithm does not support vectorial integrands
diff --git a/src/vegasflow/plain.py b/src/vegasflow/plain.py
index d9fd081..bc60e90 100644
--- a/src/vegasflow/plain.py
+++ b/src/vegasflow/plain.py
@@ -2,10 +2,11 @@
     Plain implementation of the plainest possible MonteCarlo
 """
 
-from vegasflow.configflow import fone, fzero
-from vegasflow.monte_carlo import MonteCarloFlow, wrapper, sampler
 import tensorflow as tf
 
+from vegasflow.configflow import fone, fzero
+from vegasflow.monte_carlo import MonteCarloFlow, sampler, wrapper
+
 
 class PlainFlow(MonteCarloFlow):
     """
diff --git a/src/vegasflow/tests/test_algs.py b/src/vegasflow/tests/test_algs.py
index 9d18531..a2bbb1d 100644
--- a/src/vegasflow/tests/test_algs.py
+++ b/src/vegasflow/tests/test_algs.py
@@ -7,14 +7,16 @@
 everything works """
 import json
 import tempfile
-import pytest
+
 import numpy as np
+import pytest
+import tensorflow as tf
+
+from vegasflow import plain_sampler, vegas_sampler
 from vegasflow.configflow import DTYPE, run_eager
-from vegasflow.vflow import VegasFlow
 from vegasflow.plain import PlainFlow
+from vegasflow.vflow import VegasFlow
 from vegasflow.vflowplus import VegasFlowPlus
-from vegasflow import plain_sampler, vegas_sampler
-import tensorflow as tf
 
 # Test setup
 dim = 2
diff --git a/src/vegasflow/tests/test_config.py b/src/vegasflow/tests/test_config.py
index dad1136..8142b07 100644
--- a/src/vegasflow/tests/test_config.py
+++ b/src/vegasflow/tests/test_config.py
@@ -2,11 +2,13 @@
     Test that the configuration is consistent
 """
 
+import importlib
 import os
+
 import numpy as np
-import importlib
+
 import vegasflow.configflow
-from vegasflow.configflow import DTYPE, DTYPEINT, int_me, float_me
+from vegasflow.configflow import DTYPE, DTYPEINT, float_me, int_me
 
 
 def test_int_me():
diff --git a/src/vegasflow/tests/test_gradients.py b/src/vegasflow/tests/test_gradients.py
index c03642f..42a0eaa 100644
--- a/src/vegasflow/tests/test_gradients.py
+++ b/src/vegasflow/tests/test_gradients.py
@@ -4,11 +4,10 @@
 
 import numpy as np
 from pytest import mark
-
-from vegasflow import float_me, run_eager
-from vegasflow import VegasFlow, VegasFlowPlus, PlainFlow
 import tensorflow as tf
 
+from vegasflow import PlainFlow, VegasFlow, VegasFlowPlus, float_me, run_eager
+
 
 def generate_integrand(variable):
     """Generate an integrand that depends on an input variable"""
diff --git a/src/vegasflow/tests/test_misc.py b/src/vegasflow/tests/test_misc.py
index 62d50ed..4fd46ac 100644
--- a/src/vegasflow/tests/test_misc.py
+++ b/src/vegasflow/tests/test_misc.py
@@ -1,13 +1,14 @@
 """
     Miscellaneous tests that don't really fit anywhere else
 """
-import pytest
-import numpy as np
 
-from vegasflow import VegasFlow, VegasFlowPlus, PlainFlow
+import numpy as np
+import pytest
 import tensorflow as tf
 
-from .test_algs import instance_and_compile, check_is_one
+from vegasflow import PlainFlow, VegasFlow, VegasFlowPlus
+
+from .test_algs import check_is_one, instance_and_compile
 
 
 def _vector_integrand(xarr, weight=None):
diff --git a/src/vegasflow/tests/test_utils.py b/src/vegasflow/tests/test_utils.py
index cae9a4e..c403344 100644
--- a/src/vegasflow/tests/test_utils.py
+++ b/src/vegasflow/tests/test_utils.py
@@ -1,13 +1,11 @@
 """ Test the utilities """
 
 import numpy as np
-import tensorflow as tf
 import pytest
+import tensorflow as tf
 
 from vegasflow.configflow import DTYPEINT, int_me
-from vegasflow.utils import consume_array_into_indices, py_consume_array_into_indices
-
-from vegasflow.utils import generate_condition_function
+from vegasflow.utils import consume_array_into_indices, generate_condition_function, py_consume_array_into_indices
 
 
 def test_consume_array_into_indices():
diff --git a/src/vegasflow/utils.py b/src/vegasflow/utils.py
index c748546..e1b13ab 100644
--- a/src/vegasflow/utils.py
+++ b/src/vegasflow/utils.py
@@ -2,9 +2,10 @@
     This module contains tensorflow_compiled utilities
 """
 
-from vegasflow.configflow import DTYPEINT, DTYPE, float_me, int_me, fzero
 import tensorflow as tf
 
+from vegasflow.configflow import DTYPE, DTYPEINT, float_me, fzero, int_me
+
 
 @tf.function(
     input_signature=[
diff --git a/src/vegasflow/vflow.py b/src/vegasflow/vflow.py
index 48323c8..d2f3ab6 100644
--- a/src/vegasflow/vflow.py
+++ b/src/vegasflow/vflow.py
@@ -4,26 +4,27 @@
     The main interfaces of this class are the class `VegasFlow` and the
     `vegas_wrapper`
 """
+
 import json
+import logging
+
 import numpy as np
 import tensorflow as tf
 
 from vegasflow.configflow import (
+    ALPHA,
+    BINS_MAX,
     DTYPE,
     DTYPEINT,
+    float_me,
     fone,
     fzero,
-    float_me,
-    ione,
     int_me,
-    BINS_MAX,
-    ALPHA,
+    ione,
 )
-from vegasflow.monte_carlo import MonteCarloFlow, wrapper, sampler
+from vegasflow.monte_carlo import MonteCarloFlow, sampler, wrapper
 from vegasflow.utils import consume_array_into_indices
 
-import logging
-
 logger = logging.getLogger(__name__)
 
 FBINS = float_me(BINS_MAX)
@@ -198,10 +199,7 @@ def while_body(bin_weight, n_bin, cur, prev):
     prev = fzero
     for _ in range(BINS_MAX - 1):
         bin_weight, n_bin, cur, prev = tf.while_loop(
-            while_check,
-            while_body,
-            (bin_weight, n_bin, cur, prev),
-            parallel_iterations=1,
+            while_check, while_body, (bin_weight, n_bin, cur, prev), parallel_iterations=1
         )
         bin_weight -= ave_t
         delta = (cur - prev) * bin_weight / wei_t[n_bin]
diff --git a/src/vegasflow/vflowplus.py b/src/vegasflow/vflowplus.py
index 57e228f..f329dec 100644
--- a/src/vegasflow/vflowplus.py
+++ b/src/vegasflow/vflowplus.py
@@ -6,27 +6,28 @@
 
     The main interface is the `VegasFlowPlus` class.
 """
-from itertools import product
+
 from functools import partial
+from itertools import product
+import logging
+
 import numpy as np
 import tensorflow as tf
 
 from vegasflow.configflow import (
+    BETA,
+    BINS_MAX,
     DTYPE,
     DTYPEINT,
+    MAX_NEVAL_HCUBE,
+    float_me,
     fone,
     fzero,
-    float_me,
     int_me,
-    BINS_MAX,
-    BETA,
-    MAX_NEVAL_HCUBE,
 )
-from vegasflow.monte_carlo import wrapper, sampler, MonteCarloFlow
-from vegasflow.vflow import VegasFlow, importance_sampling_digest
+from vegasflow.monte_carlo import MonteCarloFlow, sampler, wrapper
 from vegasflow.utils import consume_array_into_indices
-
-import logging
+from vegasflow.vflow import VegasFlow, importance_sampling_digest
 
 logger = logging.getLogger(__name__)
 
@@ -165,11 +166,7 @@ def _digest_random_generation(self, rnds, n_ev):
         """Generate a random array for a given number of events divided in hypercubes"""
         # Get random numbers from hypercubes
         x, w, ind, segm = generate_samples_in_hypercubes(
-            rnds,
-            self._n_strat,
-            n_ev,
-            self._hypercubes,
-            self.divisions,
+            rnds, self._n_strat, n_ev, self._hypercubes, self.divisions
         )
         return x, w, ind, segm