diff --git a/.vscode/settings.json b/.vscode/settings.json index ad7af29..b061595 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,6 +1,7 @@ { "python.testing.pytestArgs": [ - "tests" + "tests", + "--no-cov" ], "python.testing.unittestEnabled": false, "python.testing.pytestEnabled": true diff --git a/pyproject.toml b/pyproject.toml index 466be3c..9be3da8 100644 --- a/pyproject.toml +++ b/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/" diff --git a/src/pyelq/component/source_model.py b/src/pyelq/component/source_model.py index 55a12c7..ad19a97 100644 --- a/src/pyelq/component/source_model.py +++ b/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) diff --git a/src/pyelq/dispersion_model/gaussian_plume.py b/src/pyelq/dispersion_model/gaussian_plume.py index 3bb3be7..e476a27 100644 --- a/src/pyelq/dispersion_model/gaussian_plume.py +++ b/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 diff --git a/tests/component/test_source_model.py b/tests/component/test_source_model.py index 1e0507e..f1937e1 100644 --- a/tests/component/test_source_model.py +++ b/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. diff --git a/tests/test_gaussian_plume.py b/tests/test_gaussian_plume.py index a919048..6dbc37b 100644 --- a/tests/test_gaussian_plume.py +++ b/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])))