Merge pull request #11 from sede-open/change_coverage

Modify threshold function for calculating coverage

.vscode/settings.json

@@ -1,6 +1,7 @@
 {
     "python.testing.pytestArgs": [
-        "tests"
+        "tests",
+        "--no-cov"
     ],
     "python.testing.unittestEnabled": false,
     "python.testing.pytestEnabled": true

pyproject.toml

@@ -8,7 +8,7 @@ build-backend = "poetry.core.masonry.api"
 
 [tool.poetry]
 name = "pyelq-sdk"
-version = "1.0.7"
+version = "1.0.8"
 description = "Package for detection, localization and quantification code."
 authors = ["Bas van de Kerkhof", "Matthew Jones", "David Randell"]
 homepage = "https://sede-open.github.io/pyELQ/"

src/pyelq/component/source_model.py

@@ -448,6 +448,10 @@ class SourceModel(Component, SourceGrouping, SourceDistribution):
         coverage_detection (float): sensor detection threshold (in ppm) to be used for coverage calculations.
         coverage_test_source (float): test source (in kg/hr) which we wish to be able to see in coverage calculation.
 
+        threshold_function (Callable): Callable function which returns a single value that defines the threshold
+            for the coupling in a lambda function form. Examples: lambda x: np.quantile(x, 0.95, axis=0),
+            lambda x: np.max(x, axis=0), lambda x: np.mean(x, axis=0). Defaults to np.quantile.
+
     """
 
     dispersion_model: GaussianPlume = field(init=False, default=None)
@@ -472,6 +476,8 @@ class SourceModel(Component, SourceGrouping, SourceDistribution):
     coverage_detection: float = 0.1
     coverage_test_source: float = 6.0
 
+    threshold_function: callable = lambda x: np.quantile(x, 0.95, axis=0)
+
     @property
     def nof_sources(self):
         """Get number of sources in the source map."""
@@ -543,7 +549,7 @@ def initialise_dispersion_model(self, sensor_object: SensorGroup):
     def screen_coverage(self):
         """Screen the initial source map for coverage."""
         in_coverage_area = self.dispersion_model.compute_coverage(
-            self.coupling, coverage_threshold=self.coverage_threshold
+            self.coupling, coverage_threshold=self.coverage_threshold, threshold_function=self.threshold_function
         )
         self.coupling = self.coupling[:, in_coverage_area]
         all_locations = self.dispersion_model.source_map.location.to_array()
@@ -623,7 +629,9 @@ def birth_function(self, current_state: dict, prop_state: dict) -> Tuple[dict, f
         prop_state = self.update_coupling_column(prop_state, int(prop_state["n_src"]) - 1)
         prop_state["alloc_s"] = np.concatenate((prop_state["alloc_s"], np.array([0], ndmin=2)), axis=0)
         in_cov_area = self.dispersion_model.compute_coverage(
-            prop_state["A"][:, -1], coverage_threshold=self.coverage_threshold
+            prop_state["A"][:, -1],
+            coverage_threshold=self.coverage_threshold,
+            threshold_function=self.threshold_function,
         )
         if not in_cov_area:
             logp_pr_g_cr = 1e10
@@ -682,7 +690,9 @@ def move_function(self, current_state: dict, update_column: int) -> dict:
         prop_state = deepcopy(current_state)
         prop_state = self.update_coupling_column(prop_state, update_column)
         in_cov_area = self.dispersion_model.compute_coverage(
-            prop_state["A"][:, update_column], coverage_threshold=self.coverage_threshold
+            prop_state["A"][:, update_column],
+            coverage_threshold=self.coverage_threshold,
+            threshold_function=self.threshold_function,
         )
         if not in_cov_area:
             prop_state = deepcopy(current_state)

src/pyelq/dispersion_model/gaussian_plume.py

@@ -561,7 +561,7 @@ def compute_coupling_ground(
 
     @staticmethod
     def compute_coverage(
-        couplings: np.ndarray, threshold_function: Callable = np.max, coverage_threshold: float = 6, **kwargs
+        couplings: np.ndarray, threshold_function: Callable, coverage_threshold: float = 6, **kwargs
     ) -> Union[np.ndarray, dict]:
         """Returns a logical vector that indicates which sources in the couplings are, or are not, within the coverage.
 
@@ -574,18 +574,17 @@ def compute_coverage(
 
         Args:
             couplings (np.ndarray): Array of coupling values. Dimensions: n_datapoints x n_sources.
-            threshold_function (Callable, optional): Callable function which returns some single value that defines the
-                                                     maximum or 'threshold' coupling. Examples: np.quantile(q=0.9),
-                                                     np.max, np.mean. Defaults to np.max.
+            threshold_function (Callable): Callable function which returns some single value that defines the
+                maximum or 'threshold' coupling. For example: np.quantile(., q=0.95)
             coverage_threshold (float, optional): The threshold value of the estimated emission rate which is
-                                                  considered to be within the coverage. Defaults to 6 kg/hr.
+                considered to be within the coverage. Defaults to 6 kg/hr.
             kwargs (dict, optional): Keyword arguments required for the threshold function.
 
         Returns:
             coverage (Union[np.ndarray, dict]): A logical array specifying which sources are within the coverage.
 
         """
-        coupling_threshold = threshold_function(couplings, axis=0, **kwargs)
+        coupling_threshold = threshold_function(couplings, **kwargs)
         no_warning_threshold = np.where(coupling_threshold <= 1e-100, 1, coupling_threshold)
         no_warning_estimated_emission_rates = np.where(coupling_threshold <= 1e-100, np.inf, 1 / no_warning_threshold)
         coverage = no_warning_estimated_emission_rates < coverage_threshold

tests/component/test_source_model.py

@@ -99,6 +99,14 @@ def test_make_sampler(source_model):
     assert isinstance(sampler_object[0], NormalNormal)
 
 
+def test_coverage_function(source_model):
+    """Test that the coverage function has defaulted correctly."""
+    random_vars = np.random.normal(0, 1, size=(10000, 1))
+    threshold_value = source_model.threshold_function(random_vars)
+    assert threshold_value.shape == (1,)
+    assert np.allclose(threshold_value, np.quantile(random_vars, 0.95))
+
+
 def test_birth_function(source_model):
     """Test the birth_function implementation, and some aspects of the reversible jump sampler.
 

tests/test_gaussian_plume.py

@@ -588,7 +588,7 @@ def test_compute_coverage():
     source_object.location.east = np.array([-10, 10])
     source_object.location.north = np.array([25, 25])
     source_object.location.up = np.array([0, 0])
-
+    threshold_function = lambda x: np.quantile(x, 0.95, axis=0)
     plume_object = GaussianPlume(source_map=source_object)
 
     couplings = np.array(
@@ -600,14 +600,18 @@ def test_compute_coverage():
         ]
     )
 
-    coverage = plume_object.compute_coverage(couplings)
+    coverage = plume_object.compute_coverage(couplings, threshold_function=threshold_function)
     assert np.all(np.equal(coverage, np.array([True, False])))
 
-    coverage = plume_object.compute_coverage(couplings, coverage_threshold=0.3)
+    coverage = plume_object.compute_coverage(couplings, threshold_function=threshold_function, coverage_threshold=0.3)
     assert np.all(np.equal(coverage, np.array([False, False])))
 
-    coverage = plume_object.compute_coverage(couplings, threshold_function=np.mean, coverage_threshold=0.3)
+    coverage = plume_object.compute_coverage(
+        couplings, threshold_function=lambda x: np.mean(x, axis=0), coverage_threshold=0.3
+    )
     assert np.all(np.equal(coverage, np.array([False, False])))
 
-    coverage = plume_object.compute_coverage(couplings, threshold_function=np.mean, coverage_threshold=6)
+    coverage = plume_object.compute_coverage(
+        couplings, threshold_function=lambda x: np.mean(x, axis=0), coverage_threshold=6
+    )
     assert np.all(np.equal(coverage, np.array([True, False])))