From 50b47d94439d0bc0d154ec027fa5a572b4f2b027 Mon Sep 17 00:00:00 2001 From: TannazHajiMohammadloo Date: Thu, 30 May 2024 15:02:40 +0000 Subject: [PATCH 1/9] source model modififcation for change in coverage --- src/pyelq/component/source_model.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/pyelq/component/source_model.py b/src/pyelq/component/source_model.py index 55a12c7..e6c9c8a 100644 --- a/src/pyelq/component/source_model.py +++ b/src/pyelq/component/source_model.py @@ -542,8 +542,12 @@ def initialise_dispersion_model(self, sensor_object: SensorGroup): def screen_coverage(self): """Screen the initial source map for coverage.""" + threshold_function = np.quantile + quantile = 0.95 in_coverage_area = self.dispersion_model.compute_coverage( - self.coupling, coverage_threshold=self.coverage_threshold + self.coupling, coverage_threshold=self.coverage_threshold, + threshold_function=threshold_function, + q=quantile ) self.coupling = self.coupling[:, in_coverage_area] all_locations = self.dispersion_model.source_map.location.to_array() @@ -622,8 +626,12 @@ 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) + threshold_function = np.quantile + quantile = 0.95 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=threshold_function, + q=quantile ) if not in_cov_area: logp_pr_g_cr = 1e10 @@ -681,8 +689,13 @@ 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) + threshold_function = np.quantile + quantile = 0.95 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=threshold_function, + q=quantile + ) if not in_cov_area: prop_state = deepcopy(current_state) From 69ec565268a86b685e25293233003e082b9372a1 Mon Sep 17 00:00:00 2001 From: "t.hajimohammadloo" Date: Fri, 21 Jun 2024 13:43:32 +0200 Subject: [PATCH 2/9] change in the source_model gaussian_plume.compute coverage and test gaussian_plume for reflect changes in the coverage function --- src/pyelq/component/source_model.py | 26 +++++++++----------- src/pyelq/dispersion_model/gaussian_plume.py | 4 +-- tests/test_gaussian_plume.py | 10 ++++---- 3 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/pyelq/component/source_model.py b/src/pyelq/component/source_model.py index e6c9c8a..b3e3096 100644 --- a/src/pyelq/component/source_model.py +++ b/src/pyelq/component/source_model.py @@ -472,6 +472,7 @@ 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.""" @@ -542,13 +543,12 @@ def initialise_dispersion_model(self, sensor_object: SensorGroup): def screen_coverage(self): """Screen the initial source map for coverage.""" - threshold_function = np.quantile - quantile = 0.95 + # threshold_function = np.quantile + # quantile = 0.95 in_coverage_area = self.dispersion_model.compute_coverage( self.coupling, coverage_threshold=self.coverage_threshold, - threshold_function=threshold_function, - q=quantile - ) + threshold_function=self.threshold_function + ) self.coupling = self.coupling[:, in_coverage_area] all_locations = self.dispersion_model.source_map.location.to_array() screened_locations = all_locations[in_coverage_area, :] @@ -626,12 +626,11 @@ 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) - threshold_function = np.quantile - quantile = 0.95 + # threshold_function = np.quantile + # quantile = 0.95 in_cov_area = self.dispersion_model.compute_coverage( prop_state["A"][:, -1], coverage_threshold=self.coverage_threshold, - threshold_function=threshold_function, - q=quantile + threshold_function=self.threshold_function ) if not in_cov_area: logp_pr_g_cr = 1e10 @@ -689,14 +688,11 @@ 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) - threshold_function = np.quantile - quantile = 0.95 + # threshold_function = np.quantile + # quantile = 0.95 in_cov_area = self.dispersion_model.compute_coverage( prop_state["A"][:, update_column], coverage_threshold=self.coverage_threshold, - threshold_function=threshold_function, - q=quantile - - ) + threshold_function=self.threshold_function) if not in_cov_area: prop_state = deepcopy(current_state) return prop_state diff --git a/src/pyelq/dispersion_model/gaussian_plume.py b/src/pyelq/dispersion_model/gaussian_plume.py index 3bb3be7..dad2c95 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. @@ -585,7 +585,7 @@ def compute_coverage( 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/test_gaussian_plume.py b/tests/test_gaussian_plume.py index a919048..a0f578e 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,14 @@ 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]))) From 89d22ceee97c279a36fbbbf93a18e31c9cf7b60c Mon Sep 17 00:00:00 2001 From: "t.hajimohammadloo" Date: Tue, 25 Jun 2024 11:31:18 +0200 Subject: [PATCH 3/9] adding the threshold function propertu to the source model to use in the coverage funtion in gaussian plume. Change the calling of coverage function based on the change and adapt gaussian plume accordingly. --- src/pyelq/component/source_model.py | 13 +++++-------- src/pyelq/dispersion_model/gaussian_plume.py | 5 ++--- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/src/pyelq/component/source_model.py b/src/pyelq/component/source_model.py index b3e3096..b3236e1 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 some single value that defines the + threshold' 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.max. Defaults to np.quantile. + """ dispersion_model: GaussianPlume = field(init=False, default=None) @@ -543,8 +547,6 @@ def initialise_dispersion_model(self, sensor_object: SensorGroup): def screen_coverage(self): """Screen the initial source map for coverage.""" - # threshold_function = np.quantile - # quantile = 0.95 in_coverage_area = self.dispersion_model.compute_coverage( self.coupling, coverage_threshold=self.coverage_threshold, threshold_function=self.threshold_function @@ -626,12 +628,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) - # threshold_function = np.quantile - # quantile = 0.95 in_cov_area = self.dispersion_model.compute_coverage( prop_state["A"][:, -1], coverage_threshold=self.coverage_threshold, - threshold_function=self.threshold_function - ) + threshold_function=self.threshold_function) if not in_cov_area: logp_pr_g_cr = 1e10 else: @@ -688,8 +687,6 @@ 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) - # threshold_function = np.quantile - # quantile = 0.95 in_cov_area = self.dispersion_model.compute_coverage( prop_state["A"][:, update_column], coverage_threshold=self.coverage_threshold, threshold_function=self.threshold_function) diff --git a/src/pyelq/dispersion_model/gaussian_plume.py b/src/pyelq/dispersion_model/gaussian_plume.py index dad2c95..46b8683 100644 --- a/src/pyelq/dispersion_model/gaussian_plume.py +++ b/src/pyelq/dispersion_model/gaussian_plume.py @@ -574,9 +574,8 @@ 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. 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. kwargs (dict, optional): Keyword arguments required for the threshold function. From 02f37f951be5f0d4452a757326695a98ca0c76e2 Mon Sep 17 00:00:00 2001 From: code_reformat <> Date: Tue, 25 Jun 2024 09:51:16 +0000 Subject: [PATCH 4/9] Automatic reformat of code Signed-off-by: code_reformat <> --- src/pyelq/component/source_model.py | 18 +++++++++++------- tests/test_gaussian_plume.py | 8 ++++++-- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/src/pyelq/component/source_model.py b/src/pyelq/component/source_model.py index b3236e1..1defaf9 100644 --- a/src/pyelq/component/source_model.py +++ b/src/pyelq/component/source_model.py @@ -477,6 +477,7 @@ class SourceModel(Component, SourceGrouping, SourceDistribution): 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.""" @@ -548,9 +549,8 @@ 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, - threshold_function=self.threshold_function - ) + 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() screened_locations = all_locations[in_coverage_area, :] @@ -629,8 +629,10 @@ 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, - threshold_function=self.threshold_function) + prop_state["A"][:, -1], + coverage_threshold=self.coverage_threshold, + threshold_function=self.threshold_function, + ) if not in_cov_area: logp_pr_g_cr = 1e10 else: @@ -688,8 +690,10 @@ 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, - threshold_function=self.threshold_function) + 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) return prop_state diff --git a/tests/test_gaussian_plume.py b/tests/test_gaussian_plume.py index a0f578e..6dbc37b 100644 --- a/tests/test_gaussian_plume.py +++ b/tests/test_gaussian_plume.py @@ -606,8 +606,12 @@ def test_compute_coverage(): 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=lambda x:np.mean(x, axis=0), 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=lambda x:np.mean(x, axis=0), 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]))) From 61283bda76e7eba8886ffb3349754ed1daad6b0e Mon Sep 17 00:00:00 2001 From: "Matthew.Jones" Date: Tue, 2 Jul 2024 09:30:55 +0200 Subject: [PATCH 5/9] Updated version number in pyproject.toml. --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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/" From 2b1fb4cab9c3abe3d4756823cffe3912d182d5ec Mon Sep 17 00:00:00 2001 From: "Matthew.Jones" Date: Tue, 2 Jul 2024 09:34:31 +0200 Subject: [PATCH 6/9] Formatted & updated docstring for threshold_function. --- src/pyelq/component/source_model.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pyelq/component/source_model.py b/src/pyelq/component/source_model.py index 1defaf9..ad19a97 100644 --- a/src/pyelq/component/source_model.py +++ b/src/pyelq/component/source_model.py @@ -448,9 +448,9 @@ 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 some single value that defines the - threshold' 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.max. Defaults to np.quantile. + 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. """ From 3ca9e6c58474662f7a1529ad53eb0f1d37c15b51 Mon Sep 17 00:00:00 2001 From: "Matthew.Jones" Date: Tue, 2 Jul 2024 10:09:43 +0200 Subject: [PATCH 7/9] Added a unit test specifically for the threshold function. --- .vscode/settings.json | 3 ++- tests/component/test_source_model.py | 8 ++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) 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/tests/component/test_source_model.py b/tests/component/test_source_model.py index 1e0507e..16ae1f7 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. From 02a8fe8082046ea3aa1d786e61dd883ceb4e20a6 Mon Sep 17 00:00:00 2001 From: code_reformat <> Date: Tue, 2 Jul 2024 08:11:17 +0000 Subject: [PATCH 8/9] Automatic reformat of code Signed-off-by: code_reformat <> --- tests/component/test_source_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/component/test_source_model.py b/tests/component/test_source_model.py index 16ae1f7..f1937e1 100644 --- a/tests/component/test_source_model.py +++ b/tests/component/test_source_model.py @@ -103,7 +103,7 @@ 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 threshold_value.shape == (1,) assert np.allclose(threshold_value, np.quantile(random_vars, 0.95)) From 0f5e24ba8c3328e0551a5890d475d16baa6cf1e6 Mon Sep 17 00:00:00 2001 From: "Matthew.Jones" Date: Tue, 2 Jul 2024 11:25:09 +0200 Subject: [PATCH 9/9] Docstring updates in gaussian_plume.py. --- src/pyelq/dispersion_model/gaussian_plume.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pyelq/dispersion_model/gaussian_plume.py b/src/pyelq/dispersion_model/gaussian_plume.py index 46b8683..e476a27 100644 --- a/src/pyelq/dispersion_model/gaussian_plume.py +++ b/src/pyelq/dispersion_model/gaussian_plume.py @@ -575,9 +575,9 @@ def compute_coverage( Args: couplings (np.ndarray): Array of coupling values. Dimensions: n_datapoints x n_sources. threshold_function (Callable): Callable function which returns some single value that defines the - maximum or 'threshold' coupling. + 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: