Unit tests and code coverage (#47) (#48)

* Added unit tests

* Added GitHub actions workflow for tests and coverage report

.github/workflows/deploy.yml

@@ -1,4 +1,3 @@
-
 name: Deploy
 
 on:

.github/workflows/tests.yml

@@ -0,0 +1,45 @@
+name: tests
+
+on:
+  push:
+    branches:
+      - main
+  pull_request:
+    branches:
+      - main
+      - develop
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+
+    strategy:
+      matrix:
+        python-version: ['3.8', '3.9', '3.10']
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v3
+
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v4
+      with:
+        python-version: ${{ matrix.python-version }}
+
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install .[dev]
+
+    - name: Run tests with coverage
+      run: |
+        pytest --cov=pyphenopop --cov-report xml:coverage.xml
+
+    - name: Upload coverage to Codecov
+      uses: codecov/codecov-action@v4
+      with:
+          token: ${{ secrets.CODECOV_TOKEN }}
+          file: ./coverage.xml
+          flags: unittests
+          name: codecov-umbrella
+          fail_ci_if_error: true

setup.cfg

@@ -31,11 +31,15 @@ install_requires =
     pandas >= 1.5.1
     tqdm >= 4.64.1
 
-
 python_requires = >=3.8.0
 
-
 packages = find:
 
+[options.extras_require]
+dev =
+    pytest
+    pytest-cov
+
+
 [options.packages.find]
 include = pyphenopop*

tests/test_get_optimization_bounds.py

@@ -0,0 +1,60 @@
+import unittest
+from pyphenopop.mixpopid import get_optimization_bounds
+
+
+class TestGetOptimizationBounds(unittest.TestCase):
+
+    def test_single_subpopulation(self):
+        num_subpop = 1
+        bounds_model = {
+            'alpha': (0.1, 1.0),
+            'b': (0.1, 1.0),
+            'E': (0.1, 1.0),
+            'n': (0.1, 1.0)
+        }
+        bounds_sigma_low = (1e-05, 5000.0)
+        bounds_sigma_high = (1e-05, 10000.0)
+
+        bnds, lb, ub = get_optimization_bounds(num_subpop, bounds_model, bounds_sigma_low, bounds_sigma_high)
+
+        expected_bnds = (
+            (0.1, 1.0), (0.1, 1.0), (0.1, 1.0), (0.1, 1.0),
+            (1e-05, 10000.0), (1e-05, 5000.0)
+        )
+        expected_lb = [0.1, 0.1, 0.1, 0.1, 1e-05, 1e-05]
+        expected_ub = [1.0, 1.0, 1.0, 1.0, 10000.0, 5000.0]
+
+        self.assertEqual(bnds, expected_bnds)
+        self.assertEqual(lb, expected_lb)
+        self.assertEqual(ub, expected_ub)
+
+    def test_multiple_subpopulations(self):
+        num_subpop = 3
+        bounds_model = {
+            'alpha': (0.1, 1.0),
+            'b': (0.1, 1.0),
+            'E': (0.1, 1.0),
+            'n': (0.1, 1.0)
+        }
+        bounds_sigma_low = (1e-05, 5000.0)
+        bounds_sigma_high = (1e-05, 10000.0)
+
+        bnds, lb, ub = get_optimization_bounds(num_subpop, bounds_model, bounds_sigma_low, bounds_sigma_high)
+
+        expected_bnds = (
+            (0.0, 0.5), (0.0, 0.5),
+            (0.1, 1.0), (0.1, 1.0), (0.1, 1.0), (0.1, 1.0),
+            (0.1, 1.0), (0.1, 1.0), (0.1, 1.0), (0.1, 1.0),
+            (0.1, 1.0), (0.1, 1.0), (0.1, 1.0), (0.1, 1.0),
+            (1e-05, 10000.0), (1e-05, 5000.0)
+        )
+        expected_lb = [0.0, 0.0, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 1e-05, 1e-05]
+        expected_ub = [0.5, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 10000.0, 5000.0]
+
+        self.assertEqual(bnds, expected_bnds)
+        self.assertEqual(lb, expected_lb)
+        self.assertEqual(ub, expected_ub)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_mixture_id.py

@@ -0,0 +1,54 @@
+import unittest
+import numpy as np
+import pandas as pd
+from pyphenopop.mixpopid import mixture_id
+import os
+
+
+class TestMixtureId(unittest.TestCase):
+
+    def setUp(self):
+        # Create a temporary CSV file with mock data
+        self.data_file = 'test_data.csv'
+        data = np.random.rand(10*2, 3)  # 10 timepoints, 3 concentrations
+        pd.DataFrame(data).to_csv(self.data_file, header=False, index=False)
+
+        self.max_subpop = 3
+        self.timepoints = np.linspace(0, 48, 10)  # 10 timepoints from 0 to 48 hours
+        self.concentrations = np.linspace(0.01, 10, 3)  # 3 concentrations from 0.01 to 10
+        self.num_replicates = 2
+        self.model = 'expo'
+        self.bounds_model = {'alpha': (0.0, 0.1), 'b': (0.0, 1.0), 'E': (1e-06, 15), 'n': (0.01, 10)}
+        self.bounds_sigma_high = (1e-05, 10000.0)
+        self.bounds_sigma_low = (1e-05, 5000.0)
+        self.optimizer_options = {'method': 'L-BFGS-B', 'options': {'disp': False, 'ftol': 1e-12}}
+        self.num_optim = 5
+        self.selection_method = 'BIC'
+
+    def test_mixture_id(self):
+        results = mixture_id(
+            self.max_subpop,
+            self.data_file,
+            self.timepoints,
+            self.concentrations,
+            self.num_replicates,
+            self.model,
+            self.bounds_model,
+            self.bounds_sigma_high,
+            self.bounds_sigma_low,
+            self.optimizer_options,
+            self.num_optim,
+            self.selection_method
+        )
+
+        self.assertIn('summary', results)
+        self.assertIn('estimated_num_populations', results['summary'])
+        self.assertIn('final_neg_log_likelihood', results['summary'])
+        self.assertIn('final_parameters', results['summary'])
+
+    def tearDown(self):
+        os.remove(self.data_file)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_neg_log_likelihood.py

@@ -0,0 +1,48 @@
+import unittest
+import numpy as np
+from pyphenopop.mixpopid import neg_log_likelihood
+
+
+class TestNegLogLikelihood(unittest.TestCase):
+
+    def setUp(self):
+        # Setting up common variables for tests
+        self.max_subpop = 2
+        self.parameters = np.array([0.5, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
+        measurements = np.array([[[10, 20, 30], [15, 25, 35]], [[12, 22, 32], [18, 28, 38]]])
+        measurements = measurements.reshape((3, 2, 2))
+        self.measurements = measurements
+        self.concvec = np.array([0.1, 0.2])
+        self.timevec = np.array([0, 24, 48])
+        self.num_replicates = 2
+        self.model = 'expo'
+        self.num_timepoints_high = np.array(1)
+        self.num_conc_high_noise = np.array(1)
+        self.num_noise_high = 1
+        self.num_noise_low = 1
+
+    def test_neg_log_likelihood(self):
+        # Test the neg_log_likelihood function with valid inputs
+        result = neg_log_likelihood(self.max_subpop, self.parameters, self.measurements, self.concvec, self.timevec,
+                                    self.num_replicates, self.model, self.num_timepoints_high, self.num_conc_high_noise,
+                                    self.num_noise_high, self.num_noise_low)
+        self.assertIsInstance(result, float)
+
+    def test_neg_log_likelihood_invalid_model(self):
+        # Test the neg_log_likelihood function with an invalid model
+        with self.assertRaises(NotImplementedError):
+            neg_log_likelihood(self.max_subpop, self.parameters, self.measurements, self.concvec, self.timevec,
+                               self.num_replicates, 'invalid_model', self.num_timepoints_high, self.num_conc_high_noise,
+                               self.num_noise_high, self.num_noise_low)
+
+    def test_neg_log_likelihood_invalid_parameters(self):
+        # Test the neg_log_likelihood function with invalid parameters length
+        invalid_parameters = np.array([0.5, 0.1, 0.2, 0.3, 0.4])
+        with self.assertRaises(KeyError):
+            neg_log_likelihood(self.max_subpop, invalid_parameters, self.measurements, self.concvec, self.timevec,
+                               self.num_replicates, self.model, self.num_timepoints_high, self.num_conc_high_noise,
+                               self.num_noise_high, self.num_noise_low)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_pop_expo.py

@@ -0,0 +1,65 @@
+import unittest
+import numpy as np
+from pyphenopop.mixpopid import pop_expo, rate_expo
+
+
+class TestPopExpo(unittest.TestCase):
+
+    def test_pop_expo_basic(self):
+        parameters = [0.1, 0.5, 1.0, 2.0]
+        concentrations = np.array([0.1, 1.0, 10.0])
+        timepoints = np.array([0, 24, 48])
+        expected_shape = (len(concentrations), len(timepoints))
+
+        result = pop_expo(parameters, concentrations, timepoints)
+
+        self.assertEqual(result.shape, expected_shape)
+        self.assertTrue(np.all(result >= 0), "All population counts should be non-negative")
+
+    def test_pop_expo_zero_timepoints(self):
+        parameters = [0.1, 0.5, 1.0, 2.0]
+        concentrations = np.array([0.1, 1.0, 10.0])
+        timepoints = np.array([0])
+        expected_shape = (len(concentrations), len(timepoints))
+
+        result = pop_expo(parameters, concentrations, timepoints)
+
+        self.assertEqual(result.shape, expected_shape)
+        self.assertTrue(np.all(result == 1), "Population counts at time zero should be 1")
+
+    def test_pop_expo_high_concentration(self):
+        parameters = [0.1, 0.5, 1.0, 2.0]
+        concentrations = np.array([1000.0])
+        timepoints = np.array([0, 24, 48])
+        expected_shape = (len(concentrations), len(timepoints))
+
+        result = pop_expo(parameters, concentrations, timepoints)
+
+        self.assertEqual(result.shape, expected_shape)
+        self.assertTrue(np.all(result >= 0), "All population counts should be non-negative")
+
+    def test_pop_expo_edge_case(self):
+        parameters = [0.0, 0.0, 0.0, 0.0]
+        concentrations = np.array([0.0])
+        timepoints = np.array([0])
+        expected_shape = (len(concentrations), len(timepoints))
+
+        result = pop_expo(parameters, concentrations, timepoints)
+
+        self.assertEqual(result.shape, expected_shape)
+        self.assertTrue(np.all(result == 1), "Population counts at time zero should be 1")
+
+    def test_pop_expo_large_timepoints(self):
+        parameters = [0.1, 0.5, 1.0, 2.0]
+        concentrations = np.array([0.1, 1.0, 10.0])
+        timepoints = np.array([0, 1000, 2000])
+        expected_shape = (len(concentrations), len(timepoints))
+
+        result = pop_expo(parameters, concentrations, timepoints)
+
+        self.assertEqual(result.shape, expected_shape)
+        self.assertTrue(np.all(result >= 0), "All population counts should be non-negative")
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_rate_expo.py

@@ -0,0 +1,24 @@
+import unittest
+import numpy as np
+from pyphenopop.mixpopid import rate_expo
+
+
+class TestRateExpo(unittest.TestCase):
+
+    def test_rate_expo_basic_functionality(self):
+        parameters = [0.1, 0.5, 2.0, 1.0]
+        concentrations = np.array([0.1, 1.0, 10.0])
+        expected_output = np.array([0.1 + np.log(0.5 + (1 - 0.5) / (1 + (0.1 / 2.0) ** 1.0)),
+                                    0.1 + np.log(0.5 + (1 - 0.5) / (1 + (1.0 / 2.0) ** 1.0)),
+                                    0.1 + np.log(0.5 + (1 - 0.5) / (1 + (10.0 / 2.0) ** 1.0))])
+        np.testing.assert_array_almost_equal(rate_expo(parameters, concentrations), expected_output)
+
+    def test_rate_expo_zero_concentration(self):
+        parameters = [0.1, 0.5, 2.0, 1.0]
+        concentrations = np.array([0.0])
+        expected_output = np.array([0.1 + np.log(0.5 + (1 - 0.5) / (1 + (0.0 / 2.0) ** 1.0))])
+        np.testing.assert_array_almost_equal(rate_expo(parameters, concentrations), expected_output)
+
+
+if __name__ == '__main__':
+    unittest.main()