Feature synthetic data (#472)

dilpath · web-flow · commit fb152000d0a8 · 2020-10-13T21:13:43.000+02:00
* Base class that implements noise generation for simulated data, and an interface for synthetic data generation, optionally with noise
* Noise scaling factor
diff --git a/petab/__init__.py b/petab/__init__.py
@@ -21,6 +21,7 @@
 from .problem import *  # noqa: F403, F401, E402
 from .sampling import *  # noqa: F403, F401, E402
 from .sbml import *  # noqa: F403, F401, E402
+from .simulate import *  # noqa: F403, F401, E402
 from .yaml import *  # noqa: F403, F401, E402
 from .version import __version__  # noqa: F401, E402
 from .format_version import __format_version__  # noqa: F401, E402
diff --git a/petab/simulate.py b/petab/simulate.py
@@ -0,0 +1,208 @@
+import abc
+import numpy as np
+import pathlib
+import pandas as pd
+import petab
+import shutil
+import sympy as sp
+import tempfile
+from typing import Dict, Optional, Union
+
+
+class Simulator(abc.ABC):
+    """
+    Base class that specific simulators should inherit.
+    Specific simulators should minimally implement the
+    `simulate_without_noise` method.
+    Example (AMICI): https://bit.ly/33SUSG4
+
+    Attributes:
+        noise_formulas:
+            The formulae that will be used to calculate the scale of noise
+            distributions.
+        petab_problem:
+            A PEtab problem, which will be simulated.
+        rng:
+            A NumPy random generator, used to sample from noise distributions.
+        temporary_working_dir:
+            Whether `working_dir` is a temporary directory, which can be
+            deleted without significant consequence.
+        working_dir:
+            All simulator-specific output files will be saved here. This
+            directory and its contents may be modified and deleted, and
+            should be considered ephemeral.
+    """
+    def __init__(self,
+                 petab_problem: petab.Problem,
+                 working_dir: Optional[Union[pathlib.Path, str]] = None):
+        """
+        Initialize the simulator with sufficient information to perform a
+        simulation. If no working directory is specified, a temporary one is
+        created.
+
+        Arguments:
+            petab_problem:
+                A PEtab problem.
+            working_dir:
+                All simulator-specific output files will be saved here. This
+                directory and its contents may be modified and deleted, and
+                should be considered ephemeral.
+        """
+        self.petab_problem = petab_problem
+
+        self.temporary_working_dir = False
+        if working_dir is None:
+            working_dir = tempfile.mkdtemp()
+            self.temporary_working_dir = True
+        if not isinstance(working_dir, pathlib.Path):
+            working_dir = pathlib.Path(working_dir)
+        self.working_dir = working_dir
+        self.working_dir.mkdir(parents=True, exist_ok=True)
+
+        self.noise_formulas = petab.calculate.get_symbolic_noise_formulas(
+            self.petab_problem.observable_df)
+        self.rng = np.random.default_rng()
+
+    def remove_working_dir(self, force: bool = False, **kwargs) -> None:
+        """
+        Remove the simulator working directory and all files within (see the
+        `__init__` method arguments).
+
+        Arguments:
+            force:
+                If True, the working directory is removed regardless of
+                whether it is a temporary directory.
+        """
+        if force or self.temporary_working_dir:
+            shutil.rmtree(self.working_dir, **kwargs)
+            if self.working_dir.is_dir():
+                print('Failed to remove the working directory: '
+                      + str(self.working_dir))
+        else:
+            print('By default, specified working directories are not removed. '
+                  'Please call this method with `force=True`, or manually '
+                  f'delete the working directory: {self.working_dir}')
+
+    @abc.abstractmethod
+    def simulate_without_noise(self) -> pd.DataFrame:
+        """
+        Simulate a PEtab problem. This is an abstract method that should be
+        implemented in a simulation package. Links to examples of this are in
+        the class docstring.
+
+        Returns:
+            Simulated data, as a PEtab measurements table, which should be
+            equivalent to replacing all values in the `petab.C.MEASUREMENT`
+            column of the measurements table (of the PEtab problem supplied to
+            the `__init__` method), with simulated values.
+        """
+
+    def simulate(
+            self,
+            noise: bool = False,
+            noise_scaling_factor: float = 1,
+            **kwargs
+    ) -> pd.DataFrame:
+        """Simulate a PEtab problem, optionally with noise.
+
+        Arguments:
+            noise: If True, noise is added to simulated data.
+            noise_scaling_factor:
+                A multiplier of the scale of the noise distribution.
+
+        Returns:
+            Simulated data, as a PEtab measurements table.
+        """
+        simulation_df = self.simulate_without_noise(**kwargs)
+        if noise:
+            simulation_df = self.add_noise(simulation_df, noise_scaling_factor)
+        return simulation_df
+
+    def add_noise(
+            self,
+            simulation_df: pd.DataFrame,
+            noise_scaling_factor: float = 1,
+    ) -> pd.DataFrame:
+        """Add noise to simulated data.
+
+        Arguments:
+            simulation_df:
+                A PEtab measurements table that contains simulated data.
+            noise_scaling_factor:
+                A multiplier of the scale of the noise distribution.
+
+        Returns:
+            Simulated data with noise, as a PEtab measurements table.
+        """
+        simulation_df_with_noise = simulation_df.copy()
+        simulation_df_with_noise[petab.C.MEASUREMENT] = [
+            sample_noise(
+                self.petab_problem,
+                row,
+                row[petab.C.MEASUREMENT],
+                self.noise_formulas,
+                self.rng,
+                noise_scaling_factor,
+            )
+            for _, row in simulation_df_with_noise.iterrows()
+        ]
+        return simulation_df_with_noise
+
+
+def sample_noise(
+        petab_problem: petab.Problem,
+        measurement_row: pd.Series,
+        simulated_value: float,
+        noise_formulas: Optional[Dict[str, sp.Expr]] = None,
+        rng: Optional[np.random.Generator] = None,
+        noise_scaling_factor: float = 1,
+) -> float:
+    """Generate a sample from a PEtab noise distribution.
+
+    Arguments:
+        petab_problem:
+            The PEtab problem used to generate the simulated value.
+            Instance of `petab.Problem`.
+        measurement_row:
+            The row in the PEtab problem measurement table that corresponds
+            to the simulated value.
+        simulated_value:
+            A simulated value without noise.
+        noise_formulas:
+            Processed noise formulas from the PEtab observables table, in the
+            form output by the `petab.calculate.get_symbolic_noise_formulas`
+            method.
+        rng:
+            A NumPy random generator.
+        noise_scaling_factor:
+            A multiplier of the scale of the noise distribution.
+
+    Returns:
+        The sample from the PEtab noise distribution.
+    """
+    if noise_formulas is None:
+        noise_formulas = petab.calculate.get_symbolic_noise_formulas(
+            petab_problem.observable_df)
+    if rng is None:
+        rng = np.random.default_rng()
+
+    noise_value = petab.calculate.evaluate_noise_formula(
+        measurement_row,
+        noise_formulas,
+        petab_problem.parameter_df,
+        simulated_value
+    )
+
+    # default noise distribution is petab.C.NORMAL
+    noise_distribution = (
+        petab_problem
+        .observable_df
+        .loc[measurement_row[petab.C.OBSERVABLE_ID]]
+        .get(petab.C.NOISE_DISTRIBUTION, petab.C.NORMAL)
+    )
+
+    # below is e.g.: `np.random.normal(loc=simulation, scale=noise_value)`
+    return getattr(rng, noise_distribution)(
+        loc=simulated_value,
+        scale=noise_value * noise_scaling_factor
+    )
diff --git a/tests/test_simulate.py b/tests/test_simulate.py
@@ -0,0 +1,169 @@
+"""Tests for petab/simulate.py."""
+import functools
+import numpy as np
+import pandas as pd
+import pathlib
+import pytest
+import scipy.stats
+from typing import Callable
+
+import petab
+from petab.C import MEASUREMENT
+
+
+class TestSimulator(petab.simulate.Simulator):
+    """Dummy simulator."""
+
+    __test__ = False
+
+    def simulate_without_noise(self) -> pd.DataFrame:
+        """Dummy simulation."""
+        return self.petab_problem.measurement_df
+
+
+@pytest.fixture
+def petab_problem() -> petab.Problem:
+    """Create a PEtab problem for use in tests."""
+    petab_yaml_path = pathlib.Path(__file__).parent.absolute() / \
+        '../doc/example/example_Fujita/Fujita.yaml'
+    return petab.Problem.from_yaml(str(petab_yaml_path))
+
+
+def test_remove_working_dir(petab_problem):
+    """Test creation and removal of a non-empty temporary working directory."""
+    simulator = TestSimulator(petab_problem)
+    # The working directory exists
+    assert pathlib.Path(simulator.working_dir).is_dir()
+    synthetic_data_df = simulator.simulate()
+    synthetic_data_df.to_csv(f'{simulator.working_dir}/test.csv', sep='\t')
+    simulator.remove_working_dir()
+    # The (non-empty) working directory is removed
+    assert not pathlib.Path(simulator.working_dir).is_dir()
+
+    # Test creation and removal of a specified working directory
+    working_dir = 'tests/test_simulate_working_dir'
+    simulator = TestSimulator(petab_problem, working_dir=working_dir)
+    # The working directory is as specified
+    assert working_dir == str(simulator.working_dir)
+    # The working directory exists
+    assert pathlib.Path(simulator.working_dir).is_dir()
+    # A user-specified working directory should not be removed unless
+    # `force=True`.
+    simulator.remove_working_dir()
+    # The user-specified working directory is not removed without `force=True`
+    assert pathlib.Path(simulator.working_dir).is_dir()
+    simulator.remove_working_dir(force=True)
+    # The user-specified working directory is removed with `force=True`
+    assert not pathlib.Path(simulator.working_dir).is_dir()
+
+    # Test creation and removal of a specified non-empty working directory
+    simulator = TestSimulator(petab_problem, working_dir=working_dir)
+    synthetic_data_df = simulator.simulate()
+    synthetic_data_df.to_csv(f'{simulator.working_dir}/test.csv', sep='\t')
+    simulator.remove_working_dir(force=True)
+    # The non-empty, user-specified directory is removed with `force=True`
+    assert not pathlib.Path(simulator.working_dir).is_dir()
+
+
+def test_add_noise(petab_problem):
+    """Test the noise generating method."""
+
+    tested_noise_distributions = {'normal', 'laplace'}
+    assert set(petab.C.NOISE_MODELS) == tested_noise_distributions, (
+        'The noise generation methods have only been tested for '
+        f'{tested_noise_distributions}. Please edit this test '
+        'to include this distribution in its tested distributions. The '
+        'appropriate SciPy distribution will need to be added to '
+        '`petab_numpy2scipy_distribution` in `_test_add_noise`.'
+    )
+
+    for distribution in tested_noise_distributions:
+        petab_problem.observable_df[petab.C.NOISE_DISTRIBUTION] = distribution
+        _test_add_noise(petab_problem)
+
+
+def _test_add_noise(petab_problem) -> None:
+    """Test the noise generating method."""
+    n_samples = 100
+    noise_scaling_factor = 0.5
+    ks_1samp_pvalue_threshold = 0.05
+    minimum_fraction_above_threshold = 0.9
+    petab_numpy2scipy_distribution = {
+        'normal': 'norm',
+        'laplace': 'laplace',
+    }
+
+    simulator = TestSimulator(petab_problem)
+    synthetic_data_df = simulator.simulate()
+
+    # Generate samples of noisy data
+    samples = []
+    for _ in range(n_samples):
+        samples.append(
+            simulator.add_noise(
+                synthetic_data_df,
+                noise_scaling_factor=noise_scaling_factor,
+            )[MEASUREMENT]
+        )
+    samples = np.array(samples)
+
+    expected_noise_values = [
+        noise_scaling_factor *
+        petab.calculate.evaluate_noise_formula(
+            row,
+            simulator.noise_formulas,
+            petab_problem.parameter_df,
+            row[MEASUREMENT],
+        )
+        for _, row in synthetic_data_df.iterrows()
+    ]
+    expected_noise_distributions = [
+        petab_problem
+        .observable_df
+        .loc[row[petab.C.OBSERVABLE_ID]]
+        .get(petab.C.NOISE_DISTRIBUTION, petab.C.NORMAL)
+        for _, row in synthetic_data_df.iterrows()
+    ]
+
+    def row2cdf(row, index) -> Callable:
+        """
+        Get and customize the appropriate cumulative distribution function from
+        `scipy.stats`.
+
+        Arguments:
+            row:
+                A row from a PEtab measurement table.
+            index:
+                The index of the row.
+
+        Returns:
+            The appropriate SciPy cumulative distribution function, setup to
+            produce noisy simulated data.
+        """
+        return functools.partial(
+            getattr(
+                scipy.stats,
+                petab_numpy2scipy_distribution[
+                    expected_noise_distributions[index]]
+            ).cdf, loc=row[MEASUREMENT], scale=expected_noise_values[index])
+
+    # Test whether the distribution of the samples is equal to the expected
+    # distribution, for each measurement.
+    results = []
+    for index, row in synthetic_data_df.iterrows():
+        r = scipy.stats.ks_1samp(
+            samples[:, index],
+            row2cdf(row, index)
+        )
+        results.append(r)
+    observed_fraction_above_threshold = (
+        sum(r.pvalue > ks_1samp_pvalue_threshold for r in results) /
+        len(results)
+    )
+    # Sufficient distributions of measurement samples are sufficiently similar
+    # to the expected distribution
+    assert (
+        observed_fraction_above_threshold > minimum_fraction_above_threshold)
+
+    simulator.remove_working_dir()
+    assert not pathlib.Path(simulator.working_dir).is_dir()
