Repository: SALib/SALib

Issue:
Add parallel analysis
Users are currently expected to apply the desired analysis method to every column in the result array.

```python
# example taken from readme
sobol_indices = [sobol.analyze(problem, Y) for Y in y.T]
```

The hybrid OO/functional approach introduced in v1.4 handles this [automatically](https://github.com/SALib/SALib/blob/4012c1313ee360d0cc3be0fba084f4f86048e2ea/src/SALib/util/problem.py#L302) but does so sequentially as in the above example.

The approach could be extended to allow parallel analyses.

Discussion:

Patch:
diff --git a/examples/Problem/multi_output.py b/examples/Problem/multi_output.py
index a308bff..717aec0 100644
--- a/examples/Problem/multi_output.py
+++ b/examples/Problem/multi_output.py
@@ -18,7 +18,7 @@ sp = ProblemSpec({
     'outputs': ['max_P', 'Utility', 'Inertia', 'Reliability']
 })
 
-(sp.sample_saltelli(1000)
+(sp.sample_saltelli(1024)
     .evaluate(lake_problem.evaluate)
     .analyze_sobol())
 
diff --git a/examples/Problem/problem_spec.py b/examples/Problem/problem_spec.py
index 3c7a49a..e250a7c 100644
--- a/examples/Problem/problem_spec.py
+++ b/examples/Problem/problem_spec.py
@@ -1,15 +1,15 @@
 """Example showing how to use the ProblemSpec approach.
 
-Showcases method chaining, and parallel model runs using
-all available processors.
+Showcases method chaining, and parallel model runs using 2 processors.
 
 The following caveats apply:
 
-1. Functions passed into `sample`, `analyze`, `evaluate` and `evaluate_*` must 
+1. Functions passed into `sample`, `analyze` and `evaluate` must 
    accept a numpy array of `X` values as the first parameter, and return a 
    numpy array of results.
-2. Parallel evaluation is only beneficial for long-running models
-3. Currently, model results must fit in memory - no on-disk caching is provided.
+2. Parallel evaluation/analysis is only beneficial for long-running models 
+   or large datasets
+3. Currently, results must fit in memory - no on-disk caching is provided.
 """
 
 from SALib.analyze import sobol
@@ -33,29 +33,40 @@ if __name__ == '__main__':
 
     # Single core example
     start = time.perf_counter()
-    (sp.sample_saltelli(25000)
+    (sp.sample_saltelli(2**15)
         .evaluate(Ishigami.evaluate)
         .analyze_sobol(calc_second_order=True, conf_level=0.95))
     print("Time taken with 1 core:", time.perf_counter() - start, '\n')
 
+    # Same above example, but passing in specific functions
+    # (sp.sample(saltelli.sample, 25000, calc_second_order=True)
+    #     .evaluate(Ishigami.evaluate)
+    #     .analyze(sobol.analyze, calc_second_order=True, conf_level=0.95))
+
     # Samples, model results and analyses can be extracted:
     # print(sp.samples)
     # print(sp.results)
     # print(sp.analysis)
     # print(sp.to_df())
 
-    # Same above, but passing in specific functions
-    # (sp.sample(saltelli.sample, 25000, calc_second_order=True)
-    #     .evaluate(Ishigami.evaluate)
-    #     .analyze(sobol.analyze, calc_second_order=True, conf_level=0.95))
+    # Can set pre-existing samples/results as needed
+    # sp.samples = some_numpy_array
+    # sp.set_samples(some_numpy_array)
+    #
+    # Using method chaining...
+    # (sp.set_samples(some_numpy_array)
+    #    .set_results(some_result_array)
+    #    .analyze_sobol(calc_second_order=True, conf_level=0.95))
 
     # Parallel example
     start = time.perf_counter()
-    (sp.sample(saltelli.sample, 25000)
+    (sp.sample(saltelli.sample, 2**15)
          # can specify number of processors to use with `nprocs`
-        .evaluate_parallel(Ishigami.evaluate, nprocs=2)
-        .analyze(sobol.analyze, calc_second_order=True, conf_level=0.95))
-    print("Time taken with all available cores:", time.perf_counter() - start, '\n')
+         # this will be capped to the number of detected processors
+         # or, in the case of analysis, the number of outputs.
+        .evaluate(Ishigami.evaluate, nprocs=2)
+        .analyze_sobol(calc_second_order=True, conf_level=0.95, nprocs=2))
+    print("Time taken with 2 cores:", time.perf_counter() - start, '\n')
 
     print(sp)
     
@@ -66,7 +77,7 @@ if __name__ == '__main__':
                'localhost:55776')
 
     start = time.perf_counter()
-    (sp.sample(saltelli.sample, 25000)
+    (sp.sample(saltelli.sample, 2**15)
         .evaluate_distributed(Ishigami.evaluate, nprocs=2, servers=servers, verbose=True)
         .analyze(sobol.analyze, calc_second_order=True, conf_level=0.95))
     print("Time taken with distributed cores:", time.perf_counter() - start, '\n')
diff --git a/src/SALib/analyze/sobol.py b/src/SALib/analyze/sobol.py
index 002ab6b..289e407 100644
--- a/src/SALib/analyze/sobol.py
+++ b/src/SALib/analyze/sobol.py
@@ -65,7 +65,11 @@ def analyze(problem, Y, calc_second_order=True, num_resamples=100,
 
     """
     if seed:
-        np.random.seed(seed)
+        # Set seed to ensure CIs are the same
+        rng = np.random.default_rng(seed).integers
+    else:
+        rng = np.random.randint
+
     # determining if groups are defined and adjusting the number
     # of rows in the cross-sampled matrix accordingly
     groups = _check_groups(problem)
@@ -90,7 +94,7 @@ def analyze(problem, Y, calc_second_order=True, num_resamples=100,
     Y = (Y - Y.mean()) / Y.std()
 
     A, B, AB, BA = separate_output_values(Y, D, N, calc_second_order)
-    r = np.random.randint(N, size=(N, num_resamples))
+    r = rng(N, size=(N, num_resamples))
     Z = norm.ppf(0.5 + conf_level / 2)
 
     if not parallel:
diff --git a/src/SALib/util/problem.py b/src/SALib/util/problem.py
index 6a2dd78..177e941 100644
--- a/src/SALib/util/problem.py
+++ b/src/SALib/util/problem.py
@@ -5,6 +5,7 @@ from types import MethodType
 from multiprocess import Pool, cpu_count
 from pathos.pp import ParallelPythonPool as pp_Pool
 from functools import partial, wraps
+import itertools as it
 
 import numpy as np
 
@@ -61,7 +62,7 @@ class ProblemSpec(dict):
     @property
     def results(self):
         return self._results
-    
+
     @results.setter
     def results(self, vals):
         val_shape = vals.shape
@@ -85,7 +86,7 @@ class ProblemSpec(dict):
                 raise ValueError(msg)
 
         self._results = vals
-    
+
     @property
     def analysis(self):
         return self._analysis
@@ -99,10 +100,10 @@ class ProblemSpec(dict):
             Sampling method to use. The given function must accept the SALib
             problem specification as the first parameter and return a numpy
             array.
-        
+
         *args : list,
             Additional arguments to be passed to `func`
-        
+
         **kwargs : dict,
             Additional keyword arguments passed to `func`
 
@@ -118,17 +119,19 @@ class ProblemSpec(dict):
 
         return self
 
-    def set_samples(self, samples):
+    def set_samples(self, samples: np.ndarray):
         """Set previous samples used."""
         self.samples = samples
 
         return self
 
-    def set_results(self, results):
+    def set_results(self, results: np.ndarray):
         """Set previously available model results."""
+        if self.samples is not None:
+            assert self.samples.shape[0] == results.shape[0], \
+                "Provided result array does not match existing number of existing samples!"
+
         self.results = results
-        # if self.samples is not None:
-        #     warnings.warn('Existing samples found - make sure these results are for those samples!')
 
         return self
 
@@ -142,10 +145,13 @@ class ProblemSpec(dict):
             The provided function is required to accept a numpy array of
             inputs as its first parameter and must return a numpy array of
             results.
-        
+
         *args : list,
             Additional arguments to be passed to `func`
-        
+
+        nprocs : int,
+            If specified, attempts to parallelize model evaluations
+
         **kwargs : dict,
             Additional keyword arguments passed to `func`
 
@@ -153,10 +159,14 @@ class ProblemSpec(dict):
         ----------
         self : ProblemSpec object
         """
+        if 'nprocs' in kwargs:
+            nprocs = kwargs.pop('nprocs')
+            return self.evaluate_parallel(func, *args, nprocs=nprocs, **kwargs)
+
         self._results = func(self._samples, *args, **kwargs)
 
         return self
-    
+
     def evaluate_parallel(self, func, *args, nprocs=None, **kwargs):
         """Evaluate model locally in parallel.
 
@@ -166,15 +176,16 @@ class ProblemSpec(dict):
         ----------
         func : function,
             Model, or function that wraps a model, to be run in parallel.
-            The provided function needs to accept a numpy array of inputs as 
+            The provided function needs to accept a numpy array of inputs as
             its first parameter and must return a numpy array of results.
-        
+
         nprocs : int,
-            Number of processors to use. Uses all available if not specified.
-        
+            Number of processors to use.
+            Capped to the number of available processors.
+
         *args : list,
             Additional arguments to be passed to `func`
-        
+
         **kwargs : dict,
             Additional keyword arguments passed to `func`
 
@@ -186,9 +197,14 @@ class ProblemSpec(dict):
 
         if self._samples is None:
             raise RuntimeError("Sampling not yet conducted")
-        
+
+        max_procs = cpu_count()
         if nprocs is None:
-            nprocs = cpu_count()
+            nprocs = max_procs
+        else:
+            if nprocs > max_procs:
+                warnings.warn(f"{nprocs} processors requested but only {max_procs} found.")
+            nprocs = min(max_procs, nprocs)
 
         # Create wrapped partial function to allow passing of additional args
         tmp_f = self._wrap_func(func, *args, **kwargs)
@@ -211,7 +227,7 @@ class ProblemSpec(dict):
         """Distribute model evaluation across a cluster.
 
         Usage Conditions:
-        * The provided function needs to accept a numpy array of inputs as 
+        * The provided function needs to accept a numpy array of inputs as
           its first parameter
         * The provided function must return a numpy array of results
 
@@ -219,19 +235,19 @@ class ProblemSpec(dict):
         ----------
         func : function,
             Model, or function that wraps a model, to be run in parallel
-        
+
         nprocs : int,
             Number of processors to use for each node. Defaults to 1.
-        
+
         servers : list[str] or None,
             IP addresses or alias for each server/node to use.
 
         verbose : bool,
             Display job execution statistics. Defaults to False.
-        
+
         *args : list,
             Additional arguments to be passed to `func`
-        
+
         **kwargs : dict,
             Additional keyword arguments passed to `func`
 
@@ -266,17 +282,19 @@ class ProblemSpec(dict):
     def analyze(self, func, *args, **kwargs):
         """Analyze sampled results using given function.
 
-
         Parameters
         ----------
         func : function,
-            Analysis method to use. The provided function must accept the 
-            problem specification as the first parameter, X values if needed, 
+            Analysis method to use. The provided function must accept the
+            problem specification as the first parameter, X values if needed,
             Y values, and return a numpy array.
-        
+
         *args : list,
             Additional arguments to be passed to `func`
-        
+
+        nprocs : int,
+            If specified, attempts to parallelize model evaluations
+
         **kwargs : dict,
             Additional keyword arguments passed to `func`
 
@@ -284,12 +302,18 @@ class ProblemSpec(dict):
         ----------
         self : ProblemSpec object
         """
+        if 'nprocs' in kwargs:
+            # Call parallel method instead
+            return self.analyze_parallel(func, *args, **kwargs)
+
         if self._results is None:
             raise RuntimeError("Model not yet evaluated")
 
         if 'X' in func.__code__.co_varnames:
             # enforce passing of X if expected
             func = partial(func, *args, X=self._samples, **kwargs)
+        else:
+            func = partial(func, *args, **kwargs)
 
         out_cols = self.get('outputs', None)
         if out_cols is None:
@@ -302,9 +326,87 @@ class ProblemSpec(dict):
         if len(self['outputs']) > 1:
             self._analysis = {}
             for i, out in enumerate(self['outputs']):
-                self._analysis[out] = func(self, *args, Y=self._results[:, i], **kwargs)
+                self._analysis[out] = func(self, Y=self._results[:, i])
         else:
-            self._analysis = func(self, *args, Y=self._results, **kwargs)
+            self._analysis = func(self, Y=self._results)
+
+        return self
+
+    def analyze_parallel(self, func, *args, nprocs=None, **kwargs):
+        """Analyze sampled results using the given function in parallel.
+
+        Parameters
+        ----------
+        func : function,
+            Analysis method to use. The provided function must accept the
+            problem specification as the first parameter, X values if needed,
+            Y values, and return a numpy array.
+
+        *args : list,
+            Additional arguments to be passed to `func`
+
+        nprocs : int,
+            Number of processors to use.
+            Capped to the number of outputs or available processors.
+
+        **kwargs : dict,
+            Additional keyword arguments passed to `func`
+
+        Returns
+        ----------
+        self : ProblemSpec object
+        """
+        warnings.warn("This is an experimental feature and may not work.")
+
+        if self._results is None:
+            raise RuntimeError("Model not yet evaluated")
+
+        if 'X' in func.__code__.co_varnames:
+            # enforce passing of X if expected
+            func = partial(func, *args, X=self._samples, **kwargs)
+        else:
+            func = partial(func, *args, **kwargs)
+
+        out_cols = self.get('outputs', None)
+        if out_cols is None:
+            if len(self._results.shape) == 1:
+                self['outputs'] = ['Y']
+            else:
+                num_cols = self._results.shape[1]
+                self['outputs'] = [f'Y{i}' for i in range(1, num_cols+1)]
+
+        # Cap number of processors used
+        Yn = len(self['outputs'])
+        if Yn == 1:
+            # Only single output, cannot parallelize
+            warnings.warn(f"Analysis was not parallelized: {nprocs} processors requested for 1 output.")
+
+            res = func(self, Y=self._results)
+        else:
+            max_procs = cpu_count()
+            if nprocs is None:
+                nprocs = max_procs
+            else:
+                nprocs = min(Yn, nprocs, max_procs)
+
+            if ptqdm_available:
+                # Display progress bar if available
+                res = p_imap(lambda y: func(self, Y=y),
+                             [self._results[:, i] for i in range(Yn)],
+                             num_cpus=nprocs)
+            else:
+                with Pool(nprocs) as pool:
+                    res = list(pool.imap(lambda y: func(self, Y=y),
+                               [self._results[:, i] for i in range(Yn)]))
+
+        # Assign by output name if more than 1 output, otherwise
+        # attach directly
+        if Yn > 1:
+            self._analysis = {}
+            for out, Si in zip(self['outputs'], list(res)):
+                self._analysis[out] = Si
+        else:
+            self._analysis = res
 
         return self
 
@@ -316,9 +418,9 @@ class ProblemSpec(dict):
         elif isinstance(an_res, dict):
             # case where analysis result is a dict of ResultDicts
             return [an.to_df() for an in list(an_res.values())]
-        
+
         raise RuntimeError("Analysis not yet conducted")
-    
+
     def plot(self):
         """Plot results.
 
@@ -346,7 +448,7 @@ class ProblemSpec(dict):
 
         p_width = max(num_cols*3, 5)
         p_height = max(num_rows*3, 6)
-        _, axes = plt.subplots(num_rows, num_cols, sharey=True, 
+        _, axes = plt.subplots(num_rows, num_cols, sharey=True,
                                figsize=(p_width, p_height))
         for res, ax in zip(self._analysis, axes):
             self._analysis[res].plot(ax=ax)
@@ -366,9 +468,9 @@ class ProblemSpec(dict):
         tmp_f = func
         if (len(args) > 0) or (len(kwargs) > 0):
             tmp_f = partial(func, *args, **kwargs)
-        
+
         return tmp_f
-    
+
     def _setup_result_array(self):
         if len(self['outputs']) > 1:
             res_shape = (len(self._samples), len(self['outputs']))
@@ -388,24 +490,25 @@ class ProblemSpec(dict):
             r_len = len(r)
             final_res[i:i+r_len] = r
             i += r_len
-        
+
         return final_res
 
     def _method_creator(self, func, method):
+        """Generate convenience methods for specified `method`."""
         @wraps(func)
         def modfunc(self, *args, **kwargs):
             return getattr(self, method)(func, *args, **kwargs)
-        
+
         return modfunc
 
     def _add_samplers(self):
         """Dynamically add available SALib samplers as ProblemSpec methods."""
         for sampler in avail_approaches(samplers):
             func = getattr(importlib.import_module('SALib.sample.{}'.format(sampler)), 'sample')
-            method_name = "sample_{}".format(sampler.replace('_sampler', ''))            
+            method_name = "sample_{}".format(sampler.replace('_sampler', ''))
 
             self.__setattr__(method_name, MethodType(self._method_creator(func, 'sample'), self))
-    
+
     def _add_analyzers(self):
         """Dynamically add available SALib analyzers as ProblemSpec methods."""
         for analyzer in avail_approaches(analyzers):
@@ -414,13 +517,25 @@ class ProblemSpec(dict):
 
             self.__setattr__(method_name, MethodType(self._method_creator(func, 'analyze'), self))
 
-    def __repr__(self):
+    def __str__(self):
         if self._samples is not None:
-            print('Samples:', self._samples.shape, "\n")
+            arr_shape = self._samples.shape
+            if len(arr_shape) == 1:
+                arr_shape = (arr_shape[0], 1)
+            nr, nx = arr_shape
+            print('Samples:')
+            print(f'\t{nx} parameters:', self['names'])
+            print(f'\t{nr} evaluations', '\n')
         if self._results is not None:
-            print('Outputs:', self._results.shape, "\n")
+            arr_shape = self._results.shape
+            if len(arr_shape) == 1:
+                arr_shape = (arr_shape[0], 1)
+            nr, ny = arr_shape
+            print('Outputs:')
+            print(f"\t{ny} outputs:", self['outputs'])
+            print(f'\t{nr} evaluations', '\n')
         if self._analysis is not None:
-            print('Analysis:\n')
+            print('Analysis:')
             an_res = self._analysis
 
             allowed_types = (list, tuple)
@@ -448,7 +563,7 @@ def _check_spec_attributes(spec: ProblemSpec):
     assert 'bounds' in spec, "Bounds not defined"
     assert len(spec['bounds']) == len(spec['names']), \
         f"""Number of bounds do not match number of names
-        Number of names: 
+        Number of names:
         {len(spec['names'])} | {spec['names']}
         ----------------
         Number of bounds: {len(spec['bounds'])}


Test patch:
diff --git a/src/SALib/test_functions/lake_problem.py b/src/SALib/test_functions/lake_problem.py
index 88cf0e2..a0792d9 100644
--- a/src/SALib/test_functions/lake_problem.py
+++ b/src/SALib/test_functions/lake_problem.py
@@ -55,7 +55,7 @@ def lake_problem(X: FLOAT_OR_ARRAY, a: FLOAT_OR_ARRAY = 0.1, q: FLOAT_OR_ARRAY =
     return X_t1
 
 
-def evaluate_lake(values: np.ndarray) -> np.ndarray:
+def evaluate_lake(values: np.ndarray, seed=101) -> np.ndarray:
     """Evaluate the Lake Problem with an array of parameter values.
 
     .. [1] Hadka, D., Herman, J., Reed, P., Keller, K., (2015).
@@ -86,15 +86,17 @@ def evaluate_lake(values: np.ndarray) -> np.ndarray:
     -------
     np.ndarray, of Phosphorus pollution over time `t`
     """
+    rng = np.random.default_rng(seed)
+
     nvars = values.shape[0]
 
     a, q, b, mean, stdev = values.T
 
     sq_mean = mean**2
     sq_std = stdev**2
-    eps = np.random.lognormal(log(sq_mean / sqrt(sq_std + sq_mean)),
-                              sqrt(log(1.0 + sq_std / sq_mean)),
-                              size=nvars)
+    eps = rng.lognormal(log(sq_mean / sqrt(sq_std + sq_mean)),
+                        sqrt(log(1.0 + sq_std / sq_mean)),
+                        size=nvars)
 
     Y = np.zeros((nvars, nvars))
     for t in range(nvars):
@@ -104,7 +106,7 @@ def evaluate_lake(values: np.ndarray) -> np.ndarray:
     return Y
 
 
-def evaluate(values: np.ndarray, nvars: int = 100):
+def evaluate(values: np.ndarray, nvars: int = 100, seed=101):
     """Evaluate the Lake Problem with an array of parameter values.
 
     Parameters
@@ -126,7 +128,7 @@ def evaluate(values: np.ndarray, nvars: int = 100):
     nsamples = len(a)
     Y = np.empty((nsamples, 4))
     for i in range(nsamples):
-        tmp = evaluate_lake(values[i, :5])
+        tmp = evaluate_lake(values[i, :5], seed=seed)
 
         a_i = a[i]
         q_i = q[i]
diff --git a/tests/test_cli_analyze.py b/tests/test_cli_analyze.py
index 411a3fc..1c4f01b 100644
--- a/tests/test_cli_analyze.py
+++ b/tests/test_cli_analyze.py
@@ -243,7 +243,7 @@ def test_sobol():
     result = subprocess.check_output(analyze_cmd, universal_newlines=True)
     result = re.sub(r'[\n\t\s]*', '', result)
 
-    expected_output = 'STST_confx10.5579470.084460x20.4421890.044082x30.2414020.028068S1S1_confx10.3105760.060494x20.4436530.054648x3-0.0129620.054765S2S2_conf(x1,x2)-0.0143970.084384(x1,x3)0.2462310.103131(x2,x3)0.0005390.064658'
+    expected_output = "STST_confx10.5579470.085851x20.4421890.041396x30.2414020.028607S1S1_confx10.3105760.059615x20.4436530.053436x3-0.0129620.053891S2S2_conf(x1,x2)-0.0143970.083679(x1,x3)0.2462310.103117(x2,x3)0.0005390.064169"
     assert len(result) > 0 and result == expected_output, \
         "Results did not match expected values:\n\n Expected: \n{} \n\n Got: \n{}".format(
             expected_output, result)
diff --git a/tests/test_problem_spec.py b/tests/test_problem_spec.py
index c826646..fc8226b 100644
--- a/tests/test_problem_spec.py
+++ b/tests/test_problem_spec.py
@@ -1,4 +1,5 @@
 import pytest
+import copy
 
 import numpy as np
 
@@ -83,6 +84,56 @@ def test_sp_setters():
        .analyze_sobol(calc_second_order=True, conf_level=0.95))
 
 
-if __name__ == '__main__':
-    test_sp_setters()
+def test_parallel_single_output():
+    # Create the SALib Problem specification
+    sp = ProblemSpec({
+        'names': ['x1', 'x2', 'x3'],
+        'groups': None,
+        'bounds': [[-np.pi, np.pi]]*3,
+        'outputs': ['Y']
+    })
+
+    # Single core example
+    (sp.sample_saltelli(2**8)
+        .evaluate(Ishigami.evaluate)
+        .analyze_sobol(calc_second_order=True, conf_level=0.95, seed=101))
+
+    # Parallel example
+    psp = copy.deepcopy(sp)
+    (psp.sample_saltelli(2**8)
+        .evaluate_parallel(Ishigami.evaluate, nprocs=2)
+        .analyze_sobol(calc_second_order=True, conf_level=0.95, nprocs=2, seed=101))
+
+    assert np.testing.assert_equal(sp.results, psp.results) is None, "Model results not equal!"
+    assert np.testing.assert_equal(sp.analysis, psp.analysis) is None, "Analysis results not equal!"
+
+
+def test_parallel_multi_output():
+    from SALib.test_functions import lake_problem
+
+    # Create the SALib Problem specification
+    sp = ProblemSpec({
+        'names': ['a', 'q', 'b', 'mean', 'stdev', 'delta', 'alpha'],
+        'bounds': [[0.0, 0.1],
+                   [2.0, 4.5],
+                   [0.1, 0.45],
+                   [0.01, 0.05],
+                   [0.001, 0.005],
+                   [0.93, 0.99],
+                   [0.2, 0.5]],
+        'outputs': ['max_P', 'Utility', 'Inertia', 'Reliability']
+    })
+
+    # Single core example
+    (sp.sample_saltelli(2**8)
+        .evaluate(lake_problem.evaluate)
+        .analyze_sobol(calc_second_order=True, conf_level=0.95, seed=101))
+
+    # Parallel example
+    psp = copy.deepcopy(sp)
+    (psp.sample_saltelli(2**8)
+        .evaluate_parallel(lake_problem.evaluate, nprocs=2)
+        .analyze_sobol(calc_second_order=True, conf_level=0.95, nprocs=2, seed=101))
 
+    assert np.testing.assert_equal(sp.results, psp.results) is None, "Model results not equal!"
+    assert np.testing.assert_equal(sp.analysis, psp.analysis) is None, "Analysis results not equal!"