diff --git a/README.md b/README.md
index 2757e06..6d0ec02 100644
--- a/README.md
+++ b/README.md
@@ -185,8 +185,8 @@ tolerated_error: float
 step_size: float
     if given, sets a constant step size rather than adapting it.
 
-parallel: int
-    how many parallel simulations to run. default : 1
+parallel: bool
+    whether to run simulations in parallel with the available ressources. default : false
 
 repeat: int
     how many simulations to run per parameter set. default : 1
diff --git a/developement_help.md b/developement_help.md
index 480d9b0..9adb5f2 100644
--- a/developement_help.md
+++ b/developement_help.md
@@ -4,4 +4,3 @@
 - add the necessary logic in the appropriate ```initialize.ensure_consistency``` subfunction
 - optional : add a default value
 - optional : add to valid varying
-- optional : add/update dependency map
\ No newline at end of file
diff --git a/scgenerator.log b/scgenerator.log
new file mode 100644
index 0000000..e69de29
diff --git a/src/scgenerator/const.py b/src/scgenerator/const.py
index 4401be5..847548b 100644
--- a/src/scgenerator/const.py
+++ b/src/scgenerator/const.py
@@ -132,7 +132,7 @@ valid_param_types = dict(
     ),
     simulation=dict(
         behaviors=behaviors,
-        parallel=integer,
+        parallel=boolean,
         raman_type=string(["measured", "agrawal", "stolen"]),
         ideal_gas=boolean,
         repeat=integer,
diff --git a/src/scgenerator/defaults.py b/src/scgenerator/defaults.py
index fc3fbdf..586d120 100644
--- a/src/scgenerator/defaults.py
+++ b/src/scgenerator/defaults.py
@@ -20,7 +20,7 @@ default_parameters = dict(
     frep=80e6,
     behaviors=["spm", "ss"],
     raman_type="agrawal",
-    parallel=1,
+    parallel=False,
     repeat=1,
     tolerated_error=1e-11,
     lower_wavelength_interp_limit=0,
diff --git a/src/scgenerator/initialize.py b/src/scgenerator/initialize.py
index 095eabc..033fd9c 100644
--- a/src/scgenerator/initialize.py
+++ b/src/scgenerator/initialize.py
@@ -14,8 +14,7 @@ from .utils import varying_iterator, count_variations
 
 class ParamSequence(Mapping):
     def __init__(self, config):
-        validate_types(config)
-        self.config = ensure_consistency(config)
+        self.config = validate(config)
         self.name = self.config["name"]
 
         self.num_sim, self.num_varying = count_variations(self.config)
@@ -24,8 +23,10 @@ class ParamSequence(Mapping):
     def __iter__(self) -> Iterator[Tuple[list, dict]]:
         """iterates through all possible parameters, yielding a config as welle as a flattened
         computed parameters set each time"""
-        for only_varying, full_config in varying_iterator(self.config):
-            yield only_varying, compute_init_parameters(full_config)
+        for varying_only, full_config in varying_iterator(self.config):
+            for i in range(self["simulation", "repeat"]):
+                varying = varying_only + [("num", i)]
+                yield varying, compute_init_parameters(full_config)
 
     def __len__(self):
         return self.num_sim
@@ -48,19 +49,24 @@ class RecoveryParamSequence(ParamSequence):
 
     def __iter__(self) -> Iterator[Tuple[list, dict]]:
         for varying_only, full_config in varying_iterator(self.config):
-            sub_folder = os.path.join(
-                io.get_data_folder(self.id), utils.format_varying_list(varying_only)
-            )
+            for i in range(self["simulation", "repeat"]):
+                varying = varying_only + [("num", i)]
+                print("varying ", varying_only, i)
+                sub_folder = os.path.join(
+                    io.get_data_folder(self.id), utils.format_varying_list(varying)
+                )
 
-            print(f"{io.propagation_initiated(sub_folder)=}, {sub_folder=}")
-            continue
+                if not io.propagation_initiated(sub_folder):
+                    yield varying, compute_init_parameters(full_config)
+                elif not io.propagation_completed(sub_folder, self.config["simulation"]["z_num"]):
+                    yield varying, recover_params(full_config, varying, self.id)
+                else:
+                    continue
 
-            if not io.propagation_initiated(vary_str):
-                yield varying_only, compute_init_parameters(full_config)
-            elif not io.propagation_completed(vary_str):
-                yield varying_only, recover_params(full_config, varying_only, self.id)
-            else:
-                continue
+
+def validate(config: dict) -> dict:
+    _validate_types(config)
+    return _ensure_consistency(config)
 
 
 def wspace(t, t_num=0):
@@ -142,7 +148,7 @@ def validate_single_parameter(parent, key, value):
     return
 
 
-def validate_types(config):
+def _validate_types(config):
     """validates the data types in the initial config dictionary
 
     Parameters
@@ -331,7 +337,7 @@ def _ensure_consistency_simulation(simulation):
     return simulation
 
 
-def ensure_consistency(config):
+def _ensure_consistency(config):
     """ensure the config dictionary is consistent and that certain parameters are set,
     either by filling in defaults or by raising an error. This is not where new values are calculated.
 
@@ -346,7 +352,7 @@ def ensure_consistency(config):
         the consistent config dict
     """
 
-    validate_types(config)
+    _validate_types(config)
 
     # ensure parameters are not specified multiple times
     for sub_dict in valid_param_types.values():
diff --git a/src/scgenerator/io.py b/src/scgenerator/io.py
index 8192c5d..a155c1b 100644
--- a/src/scgenerator/io.py
+++ b/src/scgenerator/io.py
@@ -16,7 +16,7 @@ from .errors import IncompleteDataFolderError
 from .logger import get_logger
 
 
-def load_toml(path):
+def load_toml(path: str):
     """returns a dictionary parsed from the specified toml file"""
     if not path.lower().endswith(".toml"):
         path += ".toml"
@@ -314,7 +314,6 @@ def check_data_integrity(sub_folders: List[str], init_z_num: int):
 
 
 def propagation_initiated(sub_folder) -> bool:
-    print(f"{sub_folder=}")
     if os.path.isdir(sub_folder):
         return find_last_spectrum_file(sub_folder) > 0
     return False
@@ -342,7 +341,7 @@ def propagation_completed(sub_folder: str, init_z_num: int):
     """
     params = load_toml(os.path.join(sub_folder, "params.toml"))
     z_num = params["z_num"]
-    num_spectra = find_last_spectrum_file(sub_folder)
+    num_spectra = find_last_spectrum_file(sub_folder) + 1  # because of zero-indexing
 
     if z_num != init_z_num:
         raise IncompleteDataFolderError(
diff --git a/src/scgenerator/physics/simulate.py b/src/scgenerator/physics/simulate.py
index d216984..ea9872a 100644
--- a/src/scgenerator/physics/simulate.py
+++ b/src/scgenerator/physics/simulate.py
@@ -1,6 +1,6 @@
 import os
 from datetime import datetime
-from typing import List, Tuple
+from typing import List, Tuple, Type
 
 import numpy as np
 from numpy.fft import fft, ifft
@@ -8,6 +8,7 @@ from numpy.fft import fft, ifft
 from .. import initialize, io, utils
 from ..logger import get_logger
 from . import pulse
+from ..errors import IncompleteDataFolderError
 from .fiber import create_non_linear_op, fast_dispersion_op
 
 using_ray = False
@@ -21,6 +22,51 @@ except ModuleNotFoundError:
 
 class RK4IP:
     def __init__(self, sim_params, save_data=False, job_identifier="", task_id=0, n_percent=10):
+        """A 1D solver using 4th order Runge-Kutta in the interaction picture
+
+        Parameters
+        ----------
+        sim_params : dict
+        a flattened parameter dictionary containing :
+            w_c : numpy.ndarray
+                angular frequencies centered around 0 generated with scgenerator.initialize.wspace
+            w0 : float
+                central angular frequency of the pulse
+            w_power_fact : numpy.ndarray
+                precomputed factorial/power operations on w_c (scgenerator.math.power_fact)
+            spec_0 : numpy.ndarray
+                initial spectral envelope as function of w_c
+            z_targets : list
+                target distances
+            length : float
+                length of the fiber
+            beta : numpy.ndarray or Callable[[float], numpy.ndarray]
+                beta coeficients (Taylor expansion of beta(w))
+            gamma : float or Callable[[float], float]
+                non-linear parameter
+            t : numpy.ndarray
+                time
+            dt : float
+                time resolution
+            behaviors : list(str {'ss', 'raman', 'spm'})
+                behaviors to include in the simulation given as a list of strings
+            raman_type : str, optional
+                type of raman modelisation if raman effect is present
+            f_r, hr_w : (opt) arguments of delayed_raman_t (see there for infos)
+            adapt_step_size : bool, optional
+                if True (default), adapts the step size with conserved quantity methode
+            error_ok : float
+            tolerated relative error for the adaptive step size if adaptive
+            step size is turned on, otherwise length of fixed steps in m
+        save_data : bool, optional
+            save calculated spectra to disk, by default False
+        job_identifier : str, optional
+            string  identifying the parameter set, by default ""
+        task_id : int, optional
+            unique identifier of the session, by default 0
+        n_percent : int, optional
+            print/log progress update every n_percent, by default 10
+        """
 
         self.job_identifier = job_identifier
         self.id = task_id
@@ -31,7 +77,7 @@ class RK4IP:
         self.save_data = save_data
         self._extract_params(sim_params)
         self._setup_functions()
-        self.starting_num = sim_params.get("recovery_last_store", 1) - 1
+        self.starting_num = sim_params.get("recovery_last_stored", 0)
         self._setup_sim_parameters()
 
     def _extract_params(self, params):
@@ -79,7 +125,7 @@ class RK4IP:
         # Initial setup of simulation parameters
         self.d_w = self.w_c[1] - self.w_c[0]  # resolution of the frequency grid
         self.z = self.z_targets.pop(0)
-        self.z_stored = [self.z]  # position of each stored spectrum (for display)
+        self.z_stored = list(self.z_targets.copy()[0 : self.starting_num + 1])
 
         self.progress_tracker = utils.ProgressTracker(
             self.z_final, percent_incr=self.n_percent, logger=self.logger
@@ -151,6 +197,7 @@ class RK4IP:
     #         self.initial_h = self.error_ok
 
     def run(self):
+
         # Print introduction
         self.logger.info(
             "Computing {} new spectra, first one at {}m".format(self.store_num, self.z_targets[0])
@@ -312,38 +359,50 @@ class Simulations:
         self.logger = io.get_logger(__name__)
         self.id = int(task_id)
 
-        self.param_seq = param_seq
-        self.name = param_seq.name
+        self.update(param_seq)
+
+        self.name = self.param_seq.name
         self.data_folder = io.get_data_folder(self.id, name_if_new=self.name)
         io.save_toml(os.path.join(self.data_folder, "initial_config.toml"), self.param_seq.config)
 
-        self.using_ray = False
-        self.sim_jobs = 1
+        self.sim_jobs_per_node = 1
 
         self.propagator = RK4IP
 
+    @property
+    def finished_and_complete(self):
+        try:
+            io.check_data_integrity(
+                io.get_data_subfolders(self.data_folder), self.param_seq["simulation", "z_num"]
+            )
+            return True
+        except IncompleteDataFolderError:
+            return False
+
+    def update(self, param_seq):
+        self.param_seq = param_seq
         self.progress_tracker = utils.ProgressTracker(
             len(self.param_seq), percent_incr=1, logger=self.logger
         )
 
     def run(self):
-        for varying_params, params in self.param_seq:
-            for i in range(self.param_seq["simulation", "repeat"]):
-                varying = varying_params + [("num", i)]
-                io.save_parameters(
-                    params,
-                    io.generate_file_path(
-                        "params.toml", self.id, utils.format_varying_list(varying)
-                    ),
-                )
-                self.new_sim(varying, params.copy())
+        self._run_available()
+        self.ensure_finised_and_complete()
 
-        self.finish()
         self.logger.info(f"Merging data...")
-
         self.merge_data()
+
         self.logger.info(f"Finished simulations from config {self.name} !")
 
+    def _run_available(self):
+        for varying, params in self.param_seq:
+            io.save_parameters(
+                params,
+                io.generate_file_path("params.toml", self.id, utils.format_varying_list(varying)),
+            )
+            self.new_sim(varying, params)
+        self.finish()
+
     def new_sim(self, varying_list: List[tuple], params: dict):
         """responsible to launch a new simulation
 
@@ -361,6 +420,12 @@ class Simulations:
         """called once all the simulations are launched."""
         raise NotImplementedError()
 
+    def ensure_finised_and_complete(self):
+        while not self.finished_and_complete:
+            self.logger.warning(f"Something wrong happened, running again to finish simulation")
+            self.update(initialize.RecoveryParamSequence(self.param_seq.config, self.id))
+            self._run_available()
+
     def stop(self):
         raise NotImplementedError()
 
@@ -380,10 +445,10 @@ class SequencialSimulations(Simulations, available=True, priority=0):
         ).run()
         self.progress_tracker.update()
 
-    def finish(self):
+    def stop(self):
         pass
 
-    def stop(self):
+    def finish(self):
         pass
 
 
@@ -396,22 +461,28 @@ class RaySimulations(Simulations, available=using_ray, priority=1):
 
     def _init_ray(self):
         nodes = ray.nodes()
-        nodes_num = len(nodes)
         self.logger.info(
-            f"{nodes_num} node{'s' if nodes_num > 1 else ''} in the Ray cluster : "
-            + str([node.get("NodeManagerHostname", "unknown") for node in nodes])
+            f"{len(nodes)} node{'s' if len(nodes) > 1 else ''} in the Ray cluster : "
+            + str(
+                [
+                    (node.get("NodeManagerHostname", "unknown"), node.get("Resources", {}))
+                    for node in nodes
+                ]
+            )
         )
 
-        self.sim_jobs = min(self.param_seq.num_sim, self.param_seq["simulation", "parallel"])
         self.propagator = ray.remote(self.propagator).options(
             override_environment_variables=io.get_all_environ()
         )
-
+        self.sim_jobs_per_node = min(
+            self.param_seq.num_sim, self.param_seq["simulation", "parallel"]
+        )
+        self.update_cluster_frequency = 5
         self.jobs = []
         self.actors = {}
 
     def new_sim(self, varying_list: List[tuple], params: dict):
-        while len(self.jobs) >= self.sim_jobs:
+        while len(self.jobs) >= self.sim_jobs_total:
             self._collect_1_job()
 
         v_list_str = utils.format_varying_list(varying_list)
@@ -431,267 +502,62 @@ class RaySimulations(Simulations, available=using_ray, priority=1):
             self._collect_1_job()
 
     def _collect_1_job(self):
-        ready, self.jobs = ray.wait(self.jobs)
-        ray.get(ready)
+        ready, self.jobs = ray.wait(self.jobs, timeout=self.update_cluster_frequency)
+
+        if len(ready) == 0:
+            return
+
+        try:
+            ray.get(ready)
+            self.progress_tracker.update()
+        except Exception as e:
+            self.logger.warning("A problem occured with 1 or more worker :")
+            self.logger.warning(e)
+            ray.kill(self.actors[ready[0].task_id()])
+
         del self.actors[ready[0].task_id()]
-        self.progress_tracker.update()
 
     def stop(self):
         ray.shutdown()
 
+    @property
+    def sim_jobs_total(self):
+        tot_cpus = sum([node.get("Resources", {}).get("CPU", 0) for node in ray.nodes()])
+        return min(self.param_seq.num_sim, tot_cpus)
 
-def new_simulations(config_file: str, task_id: int, data_folder="scgenerator/"):
+
+def new_simulations(
+    config_file: str, task_id: int, data_folder="scgenerator/", Method: Type[Simulations] = None
+):
 
     config = io.load_toml(config_file)
     param_seq = initialize.ParamSequence(config)
 
-    return _new_simulations(param_seq, task_id, data_folder)
+    return _new_simulations(param_seq, task_id, data_folder, Method)
 
 
-def resume_simulations(data_folder: str, task_id: int = 0):
+def resume_simulations(
+    data_folder: str, task_id: int = 0, Method: Type[Simulations] = None
+) -> Simulations:
 
     config = io.load_toml(os.path.join(data_folder, "initial_config.toml"))
     io.set_data_folder(task_id, data_folder)
     param_seq = initialize.RecoveryParamSequence(config, task_id)
 
-    return _new_simulations(param_seq, task_id, data_folder)
+    return _new_simulations(param_seq, task_id, data_folder, Method)
 
 
-def _new_simulations(param_seq: initialize.ParamSequence, task_id, data_folder):
-    if param_seq.num_sim > 1 and param_seq["simulation", "parallel"] > 1 and using_ray:
+def _new_simulations(
+    param_seq: initialize.ParamSequence, task_id, data_folder, Method: Type[Simulations]
+):
+    if Method is not None:
+        return Method(param_seq, task_id, data_folder=data_folder)
+    elif param_seq.num_sim > 1 and param_seq["simulation", "parallel"] and using_ray:
         return Simulations.get_best_method()(param_seq, task_id, data_folder=data_folder)
     else:
         return SequencialSimulations(param_seq, task_id, data_folder=data_folder)
 
 
-def RK4IP_func(sim_params, save_data=False, job_identifier="", task_id=0, n_percent=10):
-    """Computes the spectrum of a pulse as it propagates through a PCF
-
-    Parameters
-    ----------
-        sim_params : a dictionary containing the following :
-            w_c : array
-                angular frequencies centered around 0 generated with scgenerator.initialize.wspace
-            w0 : float
-                central angular frequency of the pulse
-            t : array
-                time
-            dt : float
-                time resolution
-            spec_0 : array
-                initial spectral envelope as function of w_c
-            z_targets : list
-                target distances
-            beta : array
-                beta coeficients (Taylor expansion of beta(w))
-            gamma : float
-                non-linear parameter
-            behaviors : list(str {'ss', 'raman', 'spm'})
-                behaviors to include in the simulation given as a list of strings
-            raman_type : str, optional
-                type of raman modelisation if raman effect is present
-            f_r, hr_w : (opt) arguments of delayed_raman_t (see there for infos)
-            adapt_step_size : bool, optional
-                if True (default), adapts the step size with conserved quantity methode
-            error_ok : float
-                tolerated relative error for the adaptive step size if adaptive
-                step size is turned on, otherwise length of fixed steps in m
-        save_data : bool
-            False : return the spectra (recommended, save manually later if necessary)
-            True : save in a temporary folder and return the folder name
-                   to be used for merging later
-        job_id : int
-            id of this particular simulation
-        param_id : int
-            id corresponding to the set of paramters. Files created with the same param_id will be
-            merged if an indexer is passed (this feature is mainly used for automated parallel simulations
-            using the parallel_simulations function).
-        task_id : int
-            id of the whole program (useful when many python instances run at once). None if not running in parallel
-        n_percent : int, float
-            log message every n_percent of the simulation done
-        pt : scgenerator.progresstracker.ProgressTracker object
-        indexer : indexer object
-        debug_return : bool
-            if True and save_data False, will return photon number and step sizes as well as the spectra.
-    Returns
-    ----------
-        stored_spectra : (z_num, nt) array
-            spectrum aligned on w_c array
-        h_stored : 1D array
-            length of each valid step
-        cons_qty : 1D array
-            conserved quantity at each valid step
-
-    """
-    # DEBUG
-    debug = False
-
-    w_c = sim_params.pop("w_c")
-    w0 = sim_params.pop("w0")
-    w_power_fact = sim_params.pop("w_power_fact")
-    spec_0 = sim_params.pop("spec_0")
-    z_targets = sim_params.pop("z_targets")
-    z_final = sim_params.pop("length")
-    beta = sim_params.pop("beta_func", sim_params.pop("beta"))
-    gamma = sim_params.pop("gamma_func", sim_params.pop("gamma"))
-    behaviors = sim_params.pop("behaviors")
-    raman_type = sim_params.pop("raman_type", "stolen")
-    f_r = sim_params.pop("f_r", 0)
-    hr_w = sim_params.pop("hr_w", None)
-    adapt_step_size = sim_params.pop("adapt_step_size", True)
-    error_ok = sim_params.pop("error_ok")
-    dynamic_dispersion = sim_params.pop("dynamic_dispersion", False)
-    del sim_params
-
-    logger = get_logger(job_identifier)
-
-    # Initial setup of both non linear and linear operators
-    N_func = create_non_linear_op(behaviors, w_c, w0, gamma, raman_type, f_r, hr_w)
-    if dynamic_dispersion:
-        disp = lambda r: fast_dispersion_op(w_c, beta(r), w_power_fact)
-    else:
-        disp = lambda r: fast_dispersion_op(w_c, beta, w_power_fact)
-
-    # Set up which quantity is conserved for adaptive step size
-    if adapt_step_size:
-        if "raman" in behaviors:
-            conserved_quantity_func = pulse.photon_number
-        else:
-            print("energy conserved")
-            conserved_quantity_func = pulse.pulse_energy
-    else:
-        conserved_quantity_func = lambda a, b, c, d: 0
-
-    # making sure to keep only the z that we want
-    z_targets = list(z_targets.copy())
-    z_targets.sort()
-    store_num = len(z_targets)
-
-    # Initial setup of simulation parameters
-    d_w = w_c[1] - w_c[0]  # resolution of the frequency grid
-    z = z_targets.pop(0)
-    z_stored = [z]  # position of each stored spectrum (for display)
-
-    pt = utils.ProgressTracker(z_final, percent_incr=n_percent, logger=logger)
-
-    # Setup initial values for every physical quantity that we want to track
-    current_spectrum = spec_0.copy()
-    stored_spectra = [current_spectrum.copy()]
-    stored_field = [ifft(current_spectrum.copy())]
-    cons_qty = [conserved_quantity_func(current_spectrum, w_c + w0, d_w, gamma), 0]
-    size_fac = 2 ** (1 / 5)
-
-    if save_data:
-        _save_current_spectrum(current_spectrum, cons_qty, 0, task_id, job_identifier)
-
-    # Initial step size
-    if adapt_step_size:
-        h = (z_targets[0] - z) / 2
-    else:
-        h = error_ok
-    newh = h
-
-    # Print introduction
-    logger.info("Computing {} new spectra, first one at {}m".format(store_num, z_targets[0]))
-    pt.set(z)
-
-    # Start of the integration
-    step = 1
-    keep = True  # keep a step
-    store = False  # store a spectrum
-    time_start = datetime.today()
-
-    while z < z_final:
-        h = newh
-        z_ratio = z / z_final
-
-        # Store Exp(h/2 * disp) to be used several times
-        expD = np.exp(h / 2 * disp(z_ratio))
-
-        # RK4 algorithm
-        A_I = expD * current_spectrum
-        k1 = expD * (h * N_func(current_spectrum, z_ratio))
-        k2 = h * N_func(A_I + k1 / 2, z_ratio)
-        k3 = h * N_func(A_I + k2 / 2, z_ratio)
-        k4 = h * N_func(expD * (A_I + k3), z_ratio)
-
-        end_spectrum = expD * (A_I + k1 / 6 + k2 / 3 + k3 / 3) + k4 / 6
-
-        # Check relative error and adjust next step size
-        if adapt_step_size:
-            cons_qty[step] = conserved_quantity_func(end_spectrum, w_c + w0, d_w, gamma)
-            curr_p_change = np.abs(cons_qty[step - 1] - cons_qty[step])
-            cons_qty_change_ok = error_ok * cons_qty[step - 1]
-
-            if curr_p_change > 2 * cons_qty_change_ok:
-                keep = False
-                newh = h / 2
-            elif cons_qty_change_ok < curr_p_change <= 2 * cons_qty_change_ok:
-                keep = True
-                newh = h / size_fac
-            elif curr_p_change < 0.1 * cons_qty_change_ok:
-                keep = True
-                newh = h * size_fac
-            else:
-                keep = True
-                newh = h
-
-        # consider storing anythin only if the step was valid
-        if keep:
-
-            # If step is accepted, z becomes the current position
-            z += h
-            step += 1
-            cons_qty.append(0)
-
-            current_spectrum = end_spectrum.copy()
-
-            # Whether the current spectrum has to be stored depends on previous step
-            if store:
-                pt.suffix = " ({} steps). z = {:.4f}, h = {:.5g}".format(step, z, h)
-                pt.set(z)
-
-                stored_spectra.append(end_spectrum)
-                stored_field.append(ifft(end_spectrum))
-                if save_data:
-                    _save_current_spectrum(
-                        end_spectrum, cons_qty, len(stored_spectra) - 1, task_id, job_identifier
-                    )
-
-                z_stored.append(z)
-                del z_targets[0]
-
-                # No more spectrum to store
-                if len(z_targets) == 0:
-                    break
-                store = False
-
-                # reset the constant step size after a spectrum is stored
-                if not adapt_step_size:
-                    newh = error_ok
-
-            # if the next step goes over a position at which we want to store
-            # a spectrum, we shorten the step to reach this position exactly
-            if z + newh >= z_targets[0]:
-                store = True
-                newh = z_targets[0] - z
-        else:
-            progress_str = f"step {step} rejected with h = {h:.4e}, doing over"
-            logger.info(progress_str)
-
-    logger.info(
-        "propagation finished in {} steps ({} seconds)".format(
-            step, (datetime.today() - time_start).total_seconds()
-        )
-    )
-
-    if save_data:
-        io.save_data(z_stored, "z.npy", task_id, job_identifier)
-
-    return stored_spectra
-
-
 def _save_current_spectrum(
     spectrum: np.ndarray, cons_qty: np.ndarray, num: int, task_id: int, job_identifier: str
 ):
diff --git a/src/scgenerator/utils.py b/src/scgenerator/utils.py
index 82cfe98..34fa414 100644
--- a/src/scgenerator/utils.py
+++ b/src/scgenerator/utils.py
@@ -12,6 +12,7 @@ from typing import Any, Callable, List, Tuple, Union
 
 import numpy as np
 import ray
+from copy import deepcopy
 
 from .const import PARAM_SEPARATOR, valid_varying
 from .logger import get_logger
@@ -98,17 +99,18 @@ class ProgressTracker:
 
 
 def count_variations(config: dict) -> Tuple[int, int]:
-    """returns True if the config specified by the config dict requires only on simulation run"""
-    num = 1
-    varying_params = 0
+    """returns (sim_num, varying_params_num) where sim_num is the total number of simulations required and
+    varying_params_num is the number of distinct parameters that will vary."""
+    sim_num = 1
+    varying_params_num = 0
 
     for section_name in valid_varying:
         for array in config.get(section_name, {}).get("varying", {}).values():
-            num *= len(array)
-            varying_params += 1
+            sim_num *= len(array)
+            varying_params_num += 1
 
-    num *= config["simulation"].get("repeat", 1)
-    return num, varying_params
+    sim_num *= config["simulation"].get("repeat", 1)
+    return sim_num, varying_params_num
 
 
 def format_varying_list(l: List[tuple]):
@@ -121,13 +123,13 @@ def format_varying_list(l: List[tuple]):
     return joints[0].join(str_list)
 
 
-def varying_list_from_path(s: str) -> List[tuple]:
-    s = s.replace("/", "")
-    str_list = s.split(PARAM_SEPARATOR)
-    out = []
-    for i in range(0, len(str_list) // 2 * 2, 2):
-        out.append((str_list[i], get_value(str_list[i + 1])))
-    return out
+# def varying_list_from_path(s: str) -> List[tuple]:
+#     s = s.replace("/", "")
+#     str_list = s.split(PARAM_SEPARATOR)
+#     out = []
+#     for i in range(0, len(str_list) // 2 * 2, 2):
+#         out.append((str_list[i], get_value(str_list[i + 1])))
+#     return out
 
 
 def format_value(value):
@@ -161,9 +163,9 @@ def get_value(s: str):
 
 
 def varying_iterator(config):
+    out = deepcopy(config)
     varying_dict = {
-        section_name: config.get(section_name, {}).pop("varying", {})
-        for section_name in valid_varying
+        section_name: out.get(section_name, {}).pop("varying", {}) for section_name in valid_varying
     }
 
     possible_keys = []
@@ -171,15 +173,13 @@ def varying_iterator(config):
 
     for section_name, section in varying_dict.items():
         for key in section:
-            arr = np.atleast_1d(varying_dict[section_name][key])
-            varying_dict[section_name][key] = arr
+            arr = varying_dict[section_name][key]
             possible_keys.append((section_name, key))
             possible_ranges.append(range(len(arr)))
 
     combinations = itertools.product(*possible_ranges)
 
     for combination in combinations:
-        out = config.copy()
         only_varying = []
         for i, key in enumerate(possible_keys):
             parameter_value = varying_dict[key[0]][key[1]][combination[i]]
diff --git a/testing/configs/compute_init_parameters/good.toml b/testing/configs/compute_init_parameters/good.toml
index 784e4b0..ebe0e03 100644
--- a/testing/configs/compute_init_parameters/good.toml
+++ b/testing/configs/compute_init_parameters/good.toml
@@ -17,7 +17,7 @@ width = 50e-15
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = false
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/single_sim/true1.toml b/testing/configs/count_variations/1sim_0vary.toml
similarity index 89%
rename from testing/configs/single_sim/true1.toml
rename to testing/configs/count_variations/1sim_0vary.toml
index 226f60b..86fc327 100644
--- a/testing/configs/single_sim/true1.toml
+++ b/testing/configs/count_variations/1sim_0vary.toml
@@ -7,12 +7,13 @@ model = "marcatili"
 gas_name = "air"
 
 [pulse]
+power = 100e3
 wavelength = 800e-9
 width = 250e-15
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 repeat = 1
 t_num = 16384
 time_window = 37e-12
diff --git a/testing/configs/count_variations/1sim_1vary.toml b/testing/configs/count_variations/1sim_1vary.toml
new file mode 100644
index 0000000..1d80f6e
--- /dev/null
+++ b/testing/configs/count_variations/1sim_1vary.toml
@@ -0,0 +1,23 @@
+[fiber]
+core_radius = 50e-6
+length = 50e-2
+model = "marcatili"
+
+[gas]
+gas_name = "air"
+
+[pulse]
+power = 100e3
+wavelength = 800e-9
+
+[pulse.varying]
+width = [250e-15]
+
+[simulation]
+behaviors = ["spm", "raman", "ss"]
+parallel = true
+repeat = 1
+t_num = 16384
+time_window = 37e-12
+tolerated_error = 1e-11
+z_num = 128
diff --git a/testing/configs/single_sim/false1.toml b/testing/configs/count_variations/2sim_0vary.toml
similarity index 89%
rename from testing/configs/single_sim/false1.toml
rename to testing/configs/count_variations/2sim_0vary.toml
index db79c2b..41f7229 100644
--- a/testing/configs/single_sim/false1.toml
+++ b/testing/configs/count_variations/2sim_0vary.toml
@@ -7,12 +7,13 @@ model = "marcatili"
 gas_name = "air"
 
 [pulse]
+power = 100e3
 wavelength = 800e-9
 width = 250e-15
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 repeat = 2
 t_num = 16384
 time_window = 37e-12
diff --git a/testing/configs/single_sim/false2.toml b/testing/configs/count_variations/2sim_1vary.toml
similarity index 90%
rename from testing/configs/single_sim/false2.toml
rename to testing/configs/count_variations/2sim_1vary.toml
index fd74e8a..ed5b9ca 100644
--- a/testing/configs/single_sim/false2.toml
+++ b/testing/configs/count_variations/2sim_1vary.toml
@@ -7,6 +7,7 @@ model = "marcatili"
 gas_name = "air"
 
 [pulse]
+soliton_num = 5
 wavelength = 800e-9
 width = 250e-15
 
@@ -15,7 +16,7 @@ shape = ["gaussian", "sech"]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 repeat = 1
 t_num = 16384
 time_window = 37e-12
diff --git a/testing/configs/ensure_consistency/bad1.toml b/testing/configs/ensure_consistency/bad1.toml
index 00ce34e..9aee427 100644
--- a/testing/configs/ensure_consistency/bad1.toml
+++ b/testing/configs/ensure_consistency/bad1.toml
@@ -20,7 +20,7 @@ intensity_noise = [0.05e-2, 0.1e-2]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/bad2.toml b/testing/configs/ensure_consistency/bad2.toml
index 443a71b..25c7b43 100644
--- a/testing/configs/ensure_consistency/bad2.toml
+++ b/testing/configs/ensure_consistency/bad2.toml
@@ -20,7 +20,7 @@ soliton_num = [1, 2, 3, 4]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/bad3.toml b/testing/configs/ensure_consistency/bad3.toml
index 46d61e6..82c6f77 100644
--- a/testing/configs/ensure_consistency/bad3.toml
+++ b/testing/configs/ensure_consistency/bad3.toml
@@ -21,7 +21,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/bad4.toml b/testing/configs/ensure_consistency/bad4.toml
index 25bda29..026c8b2 100644
--- a/testing/configs/ensure_consistency/bad4.toml
+++ b/testing/configs/ensure_consistency/bad4.toml
@@ -22,7 +22,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/bad5.toml b/testing/configs/ensure_consistency/bad5.toml
index 769f304..991abe0 100644
--- a/testing/configs/ensure_consistency/bad5.toml
+++ b/testing/configs/ensure_consistency/bad5.toml
@@ -23,7 +23,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/bad6.toml b/testing/configs/ensure_consistency/bad6.toml
index 1f1ca7b..be5bf23 100644
--- a/testing/configs/ensure_consistency/bad6.toml
+++ b/testing/configs/ensure_consistency/bad6.toml
@@ -22,7 +22,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/good1.toml b/testing/configs/ensure_consistency/good1.toml
index cc9beba..de6a009 100644
--- a/testing/configs/ensure_consistency/good1.toml
+++ b/testing/configs/ensure_consistency/good1.toml
@@ -20,7 +20,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/good2.toml b/testing/configs/ensure_consistency/good2.toml
index f660214..c9f93d6 100644
--- a/testing/configs/ensure_consistency/good2.toml
+++ b/testing/configs/ensure_consistency/good2.toml
@@ -21,7 +21,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 repeat = 4
 t_num = 16384
 time_window = 37e-12
diff --git a/testing/configs/ensure_consistency/good3.toml b/testing/configs/ensure_consistency/good3.toml
index 7189cc3..7e8a56e 100644
--- a/testing/configs/ensure_consistency/good3.toml
+++ b/testing/configs/ensure_consistency/good3.toml
@@ -19,7 +19,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/good4.toml b/testing/configs/ensure_consistency/good4.toml
index 3606e02..b357f2c 100644
--- a/testing/configs/ensure_consistency/good4.toml
+++ b/testing/configs/ensure_consistency/good4.toml
@@ -29,7 +29,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/good5.toml b/testing/configs/ensure_consistency/good5.toml
index 98670fb..48588b3 100644
--- a/testing/configs/ensure_consistency/good5.toml
+++ b/testing/configs/ensure_consistency/good5.toml
@@ -20,7 +20,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/ensure_consistency/good6.toml b/testing/configs/ensure_consistency/good6.toml
index e841e22..0778cf1 100644
--- a/testing/configs/ensure_consistency/good6.toml
+++ b/testing/configs/ensure_consistency/good6.toml
@@ -22,7 +22,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/run_simulations/full_anomalous.toml b/testing/configs/run_simulations/full_anomalous.toml
index fca1a2f..af4b1bb 100644
--- a/testing/configs/run_simulations/full_anomalous.toml
+++ b/testing/configs/run_simulations/full_anomalous.toml
@@ -24,7 +24,7 @@ wavelength = [835e-9, 830e-9]
 
 [simulation]
 dt = 1e-15
-parallel = 3
+parallel = true
 raman_type = "measured"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/validate_types/bad1.toml b/testing/configs/validate_types/bad1.toml
index 80eb7ee..6f8c7ca 100644
--- a/testing/configs/validate_types/bad1.toml
+++ b/testing/configs/validate_types/bad1.toml
@@ -21,7 +21,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/validate_types/bad2.toml b/testing/configs/validate_types/bad2.toml
index 1eb052a..30ff28b 100644
--- a/testing/configs/validate_types/bad2.toml
+++ b/testing/configs/validate_types/bad2.toml
@@ -21,7 +21,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss", "q_noise"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/validate_types/bad3.toml b/testing/configs/validate_types/bad3.toml
index fe8055b..1e1a291 100644
--- a/testing/configs/validate_types/bad3.toml
+++ b/testing/configs/validate_types/bad3.toml
@@ -21,7 +21,7 @@ width = ["gaussian", "sech"]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/validate_types/bad4.toml b/testing/configs/validate_types/bad4.toml
index 2024dad..81dcc69 100644
--- a/testing/configs/validate_types/bad4.toml
+++ b/testing/configs/validate_types/bad4.toml
@@ -29,4 +29,4 @@ tolerated_error = 1e-11
 z_num = 1
 
 [simulation.varying]
-parallel = [2, 4]
+parallel = [true, false]
diff --git a/testing/configs/validate_types/bad5.toml b/testing/configs/validate_types/bad5.toml
index 975a088..114515f 100644
--- a/testing/configs/validate_types/bad5.toml
+++ b/testing/configs/validate_types/bad5.toml
@@ -21,7 +21,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/validate_types/bad6.toml b/testing/configs/validate_types/bad6.toml
index f6925c8..6904f23 100644
--- a/testing/configs/validate_types/bad6.toml
+++ b/testing/configs/validate_types/bad6.toml
@@ -21,7 +21,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 0
 t_num = 16384
diff --git a/testing/configs/validate_types/bad7.toml b/testing/configs/validate_types/bad7.toml
index f05f548..ccf5466 100644
--- a/testing/configs/validate_types/bad7.toml
+++ b/testing/configs/validate_types/bad7.toml
@@ -19,7 +19,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/configs/validate_types/good.toml b/testing/configs/validate_types/good.toml
index 49abcb0..f8ec86c 100644
--- a/testing/configs/validate_types/good.toml
+++ b/testing/configs/validate_types/good.toml
@@ -19,7 +19,7 @@ width = [50e-15, 100e-15, 200e-15]
 
 [simulation]
 behaviors = ["spm", "raman", "ss"]
-parallel = 2
+parallel = true
 raman_type = "agrawal"
 repeat = 4
 t_num = 16384
diff --git a/testing/test_initialize.py b/testing/test_initialize.py
index 0e2fad7..a090232 100644
--- a/testing/test_initialize.py
+++ b/testing/test_initialize.py
@@ -22,33 +22,33 @@ class TestInitializeMethods(unittest.TestCase):
     def test_validate_types(self):
         conf = lambda s: load_conf("validate_types/" + s)
         with self.assertRaisesRegex(TypeError, "belong"):
-            init.validate_types(conf("bad1"))
+            init._validate_types(conf("bad1"))
 
         with self.assertRaisesRegex(TypeError, "valid list of behaviors"):
-            init.validate_types(conf("bad2"))
+            init._validate_types(conf("bad2"))
 
         with self.assertRaisesRegex(TypeError, "single, real, non-negative number"):
-            init.validate_types(conf("bad3"))
+            init._validate_types(conf("bad3"))
 
         with self.assertRaisesRegex(TypeError, "'parallel' is not a valid variable parameter"):
-            init.validate_types(conf("bad4"))
+            init._validate_types(conf("bad4"))
 
         with self.assertRaisesRegex(TypeError, "Varying parameters should be specified in a list"):
-            init.validate_types(conf("bad5"))
+            init._validate_types(conf("bad5"))
 
         with self.assertRaisesRegex(
             TypeError,
             "value '0' of type <class 'int'> for key 'repeat' is not valid, must be a strictly positive integer",
         ):
-            init.validate_types(conf("bad6"))
+            init._validate_types(conf("bad6"))
 
         with self.assertRaisesRegex(
             ValueError,
             r"Varying parameters lists should contain at least 1 element",
         ):
-            init.ensure_consistency(conf("bad7"))
+            init._ensure_consistency(conf("bad7"))
 
-        self.assertIsNone(init.validate_types(conf("good")))
+        self.assertIsNone(init._validate_types(conf("good")))
 
     def test_ensure_consistency(self):
         conf = lambda s: load_conf("ensure_consistency/" + s)
@@ -56,68 +56,68 @@ class TestInitializeMethods(unittest.TestCase):
             MissingParameterError,
             r"1 of '\['t0', 'width'\]' is required and no defaults have been set",
         ):
-            init.ensure_consistency(conf("bad1"))
+            init._ensure_consistency(conf("bad1"))
 
         with self.assertRaisesRegex(
             MissingParameterError,
             r"1 of '\['power', 'energy', 'width', 't0'\]' is required when 'soliton_num' is specified and no defaults have been set",
         ):
-            init.ensure_consistency(conf("bad2"))
+            init._ensure_consistency(conf("bad2"))
 
         with self.assertRaisesRegex(
             MissingParameterError,
             r"2 of '\['dt', 't_num', 'time_window'\]' are required and no defaults have been set",
         ):
-            init.ensure_consistency(conf("bad3"))
+            init._ensure_consistency(conf("bad3"))
 
         with self.assertRaisesRegex(
             DuplicateParameterError,
             r"got multiple values for parameter 'width'",
         ):
-            init.ensure_consistency(conf("bad4"))
+            init._ensure_consistency(conf("bad4"))
 
         with self.assertRaisesRegex(
             MissingParameterError,
             r"'capillary_thickness' is a required parameter for fiber model 'hasan' and no defaults have been set",
         ):
-            init.ensure_consistency(conf("bad5"))
+            init._ensure_consistency(conf("bad5"))
 
         with self.assertRaisesRegex(
             MissingParameterError,
             r"1 of '\['capillary_spacing', 'capillary_outer_d'\]' is required for fiber model 'hasan' and no defaults have been set",
         ):
-            init.ensure_consistency(conf("bad6"))
+            init._ensure_consistency(conf("bad6"))
 
         self.assertLessEqual(
-            {"model": "pcf"}.items(), init.ensure_consistency(conf("good1"))["fiber"].items()
+            {"model": "pcf"}.items(), init._ensure_consistency(conf("good1"))["fiber"].items()
         )
 
-        self.assertNotIn("gas", init.ensure_consistency(conf("good1")))
+        self.assertNotIn("gas", init._ensure_consistency(conf("good1")))
 
-        self.assertNotIn("gamma", init.ensure_consistency(conf("good4"))["fiber"])
+        self.assertNotIn("gamma", init._ensure_consistency(conf("good4"))["fiber"])
 
         self.assertLessEqual(
             {"raman_type": "agrawal"}.items(),
-            init.ensure_consistency(conf("good2"))["simulation"].items(),
+            init._ensure_consistency(conf("good2"))["simulation"].items(),
         )
 
         self.assertLessEqual(
-            {"name": "no name"}.items(), init.ensure_consistency(conf("good3")).items()
+            {"name": "no name"}.items(), init._ensure_consistency(conf("good3")).items()
         )
 
         self.assertLessEqual(
             {"capillary_nested": 0, "capillary_resonance_strengths": []}.items(),
-            init.ensure_consistency(conf("good4"))["fiber"].items(),
+            init._ensure_consistency(conf("good4"))["fiber"].items(),
         )
 
         self.assertLessEqual(
             dict(he_mode=(1, 1)).items(),
-            init.ensure_consistency(conf("good5"))["fiber"].items(),
+            init._ensure_consistency(conf("good5"))["fiber"].items(),
         )
 
         self.assertLessEqual(
             dict(temperature=300, pressure=1e5, gas_name="vacuum", plasma_density=0).items(),
-            init.ensure_consistency(conf("good5"))["gas"].items(),
+            init._ensure_consistency(conf("good5"))["gas"].items(),
         )
 
         self.assertLessEqual(
@@ -127,29 +127,14 @@ class TestInitializeMethods(unittest.TestCase):
                 lower_wavelength_interp_limit=0,
                 upper_wavelength_interp_limit=1900e-9,
             ).items(),
-            init.ensure_consistency(conf("good6"))["simulation"].items(),
+            init._ensure_consistency(conf("good6"))["simulation"].items(),
         )
 
-    def test_single_sim(self):
-        conf = conf_maker("single_sim")
-
-        self.assertTrue(init.single_sim(conf("true1")))
-
-        self.assertFalse(init.single_sim(conf("false1")))
-
-        self.assertFalse(init.single_sim(conf("false2")))
-
     # def test_compute_init_parameters(self):
     #     conf = lambda s: load_conf("compute_init_parameters/" + s)
 
 
 if __name__ == "__main__":
     conf = conf_maker("validate_types")
-    config = conf("good")
-    pprint(config)
-    config = init.ensure_consistency(config)
-    pprint(config)
-    params = init.compute_init_parameters(config)
-    pprint(params)
 
     unittest.main()
diff --git a/testing/test_utils.py b/testing/test_utils.py
new file mode 100644
index 0000000..036bb46
--- /dev/null
+++ b/testing/test_utils.py
@@ -0,0 +1,30 @@
+import unittest
+from scgenerator import utils, initialize
+import toml
+
+
+def load_conf(name):
+    with open("testing/configs/" + name + ".toml") as file:
+        conf = toml.load(file)
+    return conf
+
+
+def conf_maker(folder):
+    def conf(name):
+        return initialize.validate(load_conf(folder + "/" + name))
+
+    return conf
+
+
+class TestUtilsMethods(unittest.TestCase):
+    def test_count_variations(self):
+        conf = conf_maker("count_variations")
+
+        self.assertEqual((1, 0), utils.count_variations(conf("1sim_0vary")))
+        self.assertEqual((1, 1), utils.count_variations(conf("1sim_1vary")))
+        self.assertEqual((2, 1), utils.count_variations(conf("2sim_1vary")))
+        self.assertEqual((2, 0), utils.count_variations(conf("2sim_0vary")))
+
+
+if __name__ == "__main__":
+    unittest.main()
\ No newline at end of file