diff --git a/tests/test_active_learning.py b/tests/test_active_learning.py
index 00f26b8555fcdd9c700dc7cc2fec2ba6ef40f51f..e045be5e0c629a647e9e308baf3150edd82265ff 100644
--- a/tests/test_active_learning.py
+++ b/tests/test_active_learning.py
@@ -1,21 +1,19 @@
+import os
 import pytest
 import numpy as np
+import matplotlib.pyplot as plt
 from scipy import optimize
 from scipy.stats import norm, multivariate_normal, uniform
-from psimpy.inference.active_learning import ActiveLearning
-from psimpy.simulator.run_simulator import RunSimulator
-from psimpy.simulator.mass_point_model import MassPointModel
-from psimpy.sampler.latin import LHS
-from psimpy.sampler.saltelli import Saltelli
-from psimpy.emulator.robustgasp import ScalarGaSP, PPGaSP
-from psimpy.inference.bayes_inference import GridEstimation
-from psimpy.inference.bayes_inference import MetropolisHastingsEstimation
-from psimpy.sampler.metropolis_hastings import MetropolisHastings
-from scipy.stats import uniform, norm
-from scipy import optimize
 from beartype.roar import BeartypeCallHintParamViolation
-import matplotlib.pyplot as plt
-import os
+from psimpy.inference import ActiveLearning
+from psimpy.simulator import RunSimulator
+from psimpy.simulator import MassPointModel
+from psimpy.sampler import LHS
+from psimpy.sampler import Saltelli
+from psimpy.emulator import ScalarGaSP, PPGaSP
+from psimpy.inference import GridEstimation
+from psimpy.inference import MetropolisHastingsEstimation
+from psimpy.sampler import MetropolisHastings
 
 dir_test = os.path.abspath(os.path.join(__file__, '../'))
 
diff --git a/tests/test_bayes_inference.py b/tests/test_bayes_inference.py
index dacc2c3d9b4ca2114165d765426a98a91d5adbee..b1fdd35cc5cee22d50220df86e88955cfd1c5fb7 100644
--- a/tests/test_bayes_inference.py
+++ b/tests/test_bayes_inference.py
@@ -1,12 +1,12 @@
+import os
 import pytest
+import itertools
 import numpy as np
-from scipy.stats import norm, multivariate_normal, uniform
-from psimpy.inference.bayes_inference import GridEstimation
-from psimpy.inference.bayes_inference import MetropolisHastingsEstimation
-from psimpy.sampler.metropolis_hastings import MetropolisHastings
 import matplotlib.pyplot as plt
-import os
-import itertools
+from scipy.stats import norm, multivariate_normal, uniform
+from psimpy.inference import GridEstimation
+from psimpy.inference import MetropolisHastingsEstimation
+from psimpy.sampler import MetropolisHastings
 
 dir_test = os.path.abspath(os.path.join(__file__, '../'))
 
diff --git a/tests/test_latin.py b/tests/test_latin.py
index 7dd3f283473e70fc7af9fd234b877428521dab91..ca35184ed98e2015d3f17fb3f3a9e9382e68bd7b 100644
--- a/tests/test_latin.py
+++ b/tests/test_latin.py
@@ -1,6 +1,6 @@
 import pytest
 import numpy as np
-from psimpy.sampler.latin import LHS
+from psimpy.sampler import LHS
 from beartype.roar import BeartypeCallHintParamViolation
 
 @pytest.mark.parametrize(
diff --git a/tests/test_mass_point_model.py b/tests/test_mass_point_model.py
index 9381d3ce424a0f65b3ba97d89b2608225861ec40..30fd3ee6afc18ad86e0cf6d690fd286ef610ac16 100644
--- a/tests/test_mass_point_model.py
+++ b/tests/test_mass_point_model.py
@@ -1,7 +1,7 @@
-from psimpy.simulator.mass_point_model import MassPointModel
-import numpy as np
 import os
+import numpy as np
 import pytest
+from psimpy.simulator import MassPointModel
 
 @pytest.mark.parametrize(
     "elevation, x0, y0",
diff --git a/tests/test_metropolis_hastings.py b/tests/test_metropolis_hastings.py
index f8bc3a4124f305e81f4a8e3a73c4b25d64e72e9b..996f4ca884af653a726ee03aaf3de6cd0a507f82 100644
--- a/tests/test_metropolis_hastings.py
+++ b/tests/test_metropolis_hastings.py
@@ -1,9 +1,9 @@
+import os
 import pytest
 import numpy as np
-from scipy.stats import norm, multivariate_normal, uniform
-from psimpy.sampler.metropolis_hastings import MetropolisHastings
 import matplotlib.pyplot as plt
-import os
+from scipy.stats import norm, multivariate_normal, uniform
+from psimpy.sampler import MetropolisHastings
 
 @pytest.mark.parametrize(
     "ndim, init_state, f_sample, target, ln_target, bounds, f_density, symmetric",
diff --git a/tests/test_ravaflow24.py b/tests/test_ravaflow24.py
index fa7c5961ccaf95de0ec0672db014a0bd8d9271ec..cc04fe8283d575c0298a20c0a6fa4aecf1167838 100644
--- a/tests/test_ravaflow24.py
+++ b/tests/test_ravaflow24.py
@@ -1,7 +1,7 @@
-from psimpy.simulator.ravaflow24 import Ravaflow24Mixture
-import numpy as np
 import os
 import pytest
+import numpy as np
+from psimpy.simulator import Ravaflow24Mixture
 
 @pytest.mark.parametrize(
     "dir_sim, conversion_control, curvature_control, surface_control, \
diff --git a/tests/test_robustgasp.py b/tests/test_robustgasp.py
index 6e8af65352018dbb3559550a31fd73892fe9bcca..1b0113a36b5468325fa5e892f8710ea93db2958f 100644
--- a/tests/test_robustgasp.py
+++ b/tests/test_robustgasp.py
@@ -1,9 +1,10 @@
+import os
 import pytest
 import numpy as np
-from psimpy.emulator.robustgasp import RobustGaSPBase, ScalarGaSP, PPGaSP
-from beartype.roar import BeartypeCallHintParamViolation
 from rpy2.rinterface import NA
-import os
+from beartype.roar import BeartypeCallHintParamViolation
+from psimpy.emulator import ScalarGaSP, PPGaSP
+from psimpy.emulator.robustgasp import RobustGaSPBase
 
 @pytest.mark.parametrize(
     "ndim, zero_mean, nugget, nugget_est, range_par, method, prior_choice, \
diff --git a/tests/test_run_mass_point_model.py b/tests/test_run_mass_point_model.py
index 310fcce2de46a77b5f619e7c313f302b781d4ceb..26d5a79580d3c96dc8b48fad37f605ff13957ac7 100644
--- a/tests/test_run_mass_point_model.py
+++ b/tests/test_run_mass_point_model.py
@@ -1,9 +1,9 @@
-from psimpy.simulator.run_simulator import RunSimulator
-from psimpy.simulator.mass_point_model import MassPointModel
 import os
-import numpy as np
-import itertools
 import time
+import itertools
+import numpy as np
+from psimpy.simulator import RunSimulator
+from psimpy.simulator import MassPointModel
 
 def test_run_mass_point_model():
     mpm = MassPointModel()
diff --git a/tests/test_run_ravaflow24.py b/tests/test_run_ravaflow24.py
index d2b34660909cf4e318c7e7f8da7c317261c29fae..940630cc01760fa79aa8502da82f89c4446953cd 100644
--- a/tests/test_run_ravaflow24.py
+++ b/tests/test_run_ravaflow24.py
@@ -1,9 +1,9 @@
-from psimpy.simulator.run_simulator import RunSimulator
-from psimpy.simulator.ravaflow24 import Ravaflow24Mixture
-import numpy as np
-import itertools
-import time
 import os
+import time
+import itertools
+import numpy as np
+from psimpy.simulator import RunSimulator
+from psimpy.simulator import Ravaflow24Mixture
 
 dir_test = os.path.abspath(os.path.join(__file__, '../'))
 
diff --git a/tests/test_run_simulator.py b/tests/test_run_simulator.py
index 79ab5ba2ac14b77449a9bdfa8487c725d03752e3..7096e9782238d65df7a1cd012b673c34a49d2c0a 100644
--- a/tests/test_run_simulator.py
+++ b/tests/test_run_simulator.py
@@ -1,8 +1,8 @@
-from psimpy.simulator.run_simulator import RunSimulator
+import os
 import pytest
 import numpy as np
-import os
 from beartype.roar import BeartypeCallHintParamViolation
+from psimpy.simulator import RunSimulator
 
 def add(a, b, c , d=100, save=False, filename=None):
     if save is True:
diff --git a/tests/test_saltelli.py b/tests/test_saltelli.py
index d9331c53372e66c63a0d07bc4e000bb07808bcec..fffcadc986c1c314e849925582a9c4a6ed2a0d4a 100644
--- a/tests/test_saltelli.py
+++ b/tests/test_saltelli.py
@@ -1,7 +1,7 @@
 import pytest
 import numpy as np
 from SALib.sample.saltelli import sample
-from psimpy.sampler.saltelli import Saltelli
+from psimpy.sampler import Saltelli
 from beartype.roar import BeartypeCallHintParamViolation
 
 
diff --git a/tests/test_sobol.py b/tests/test_sobol.py
index f4c651006422d0acbafd0ed6577d0dd3a5bee95e..0db41bdc9b278e0705468166caea1f9a6e360aed 100644
--- a/tests/test_sobol.py
+++ b/tests/test_sobol.py
@@ -1,10 +1,9 @@
 import pytest
 import numpy as np
-from psimpy.sampler.saltelli import Saltelli
-from psimpy.sampler.latin import LHS
-from psimpy.emulator.robustgasp import ScalarGaSP
-from psimpy.sensitivity.sobol import SobolAnalyze
-
+from psimpy.sampler import Saltelli
+from psimpy.sampler import LHS
+from psimpy.emulator import ScalarGaSP
+from psimpy.sensitivity import SobolAnalyze
 
 def f(x1,x2,x3):
     return np.sin(x1) + 7*np.sin(x2)**2 + 0.1*x3**4*np.sin(x1)
@@ -42,7 +41,6 @@ def test_SobolAnalyze(calc_second_order, skip_values, nbase, seed, mode,
     print('analytical S1: \n', [0.314, 0.442, 0])
 
 
-
 @pytest.mark.parametrize(
     "calc_second_order, skip_values, nbase, seed, mode, max_workers",
     [
diff --git a/tests/test_util_funcs.py b/tests/test_util_funcs.py
index c2699960b458ebf7e07112ef81620fb9051b2d55..d8141a1d72138af7f6fa4aed7ecb3d06b467f9d2 100644
--- a/tests/test_util_funcs.py
+++ b/tests/test_util_funcs.py
@@ -1,7 +1,8 @@
-import numpy as np
 import pytest
+import numpy as np
+from beartype.roar import BeartypeCallHintParamViolation
 from psimpy.utility.util_funcs import scale_samples
-import beartype
+
 
 @pytest.mark.parametrize(
     'samples, bounds',
@@ -11,7 +12,7 @@ import beartype
      ]
 )
 def test_scale_samples_TypeError(samples, bounds):
-    with pytest.raises(beartype.roar.BeartypeCallHintParamViolation):
+    with pytest.raises(BeartypeCallHintParamViolation):
         scale_samples(samples, bounds)
 
 @pytest.mark.parametrize(