From 9ef1a8c4a8dcea54595e996d9a8ad515463c3aba Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 12:24:51 +0100 Subject: [PATCH 01/10] feat: add output classes Add classes to wrap output of algorithms --- discretesampling/base/output/__init__.py | 3 +++ discretesampling/base/output/base.py | 8 +++++++ discretesampling/base/output/mcmc_output.py | 23 +++++++++++++++++++++ discretesampling/base/output/smc_output.py | 17 +++++++++++++++ 4 files changed, 51 insertions(+) create mode 100644 discretesampling/base/output/__init__.py create mode 100644 discretesampling/base/output/base.py create mode 100644 discretesampling/base/output/mcmc_output.py create mode 100644 discretesampling/base/output/smc_output.py diff --git a/discretesampling/base/output/__init__.py b/discretesampling/base/output/__init__.py new file mode 100644 index 0000000..1a1a943 --- /dev/null +++ b/discretesampling/base/output/__init__.py @@ -0,0 +1,3 @@ +from .base import BaseOutput # no qa +from .smc_output import SMCOutput # no qa +from .mcmc_output import MCMCOutput # no qa diff --git a/discretesampling/base/output/base.py b/discretesampling/base/output/base.py new file mode 100644 index 0000000..cd9ab98 --- /dev/null +++ b/discretesampling/base/output/base.py @@ -0,0 +1,8 @@ +from abc import abc +import numpy as np + + +class BaseOutput(abc): + def __init__(self, samples): + self.samples = samples + self.log_weights = np.repeat([-len(samples)], len(samples)) diff --git a/discretesampling/base/output/mcmc_output.py b/discretesampling/base/output/mcmc_output.py new file mode 100644 index 0000000..14b8eb7 --- /dev/null +++ b/discretesampling/base/output/mcmc_output.py @@ -0,0 +1,23 @@ +import numpy as np +from discretesampling.base.types import DiscreteVariable +from discretesampling.base.output import BaseOutput +from discretesampling.base.executor import Executor + + +class MCMCOutput(BaseOutput): + def __init__( + self, + samples: list[DiscreteVariable], + acceptance_probabilities, + include_warmup, + N=None, N_warmup=None, + exec: Executor = Executor() + ): + super().__init__(self, samples) + self.acceptance_probabilities = acceptance_probabilities + self.include_warmup = include_warmup + if N is None: + N = len(samples) + self.N = N + self.N_warmup = N_warmup + self.acceptance_rate = np.mean(acceptance_probabilities) diff --git a/discretesampling/base/output/smc_output.py b/discretesampling/base/output/smc_output.py new file mode 100644 index 0000000..69883fc --- /dev/null +++ b/discretesampling/base/output/smc_output.py @@ -0,0 +1,17 @@ +import numpy as np +from discretesampling.base.types import DiscreteVariable +from discretesampling.base.output import BaseOutput +from discretesampling.base.algorithms.smc_components import ess +from discretesampling.base.executor import Executor + + +class SMCOutput(BaseOutput): + def __init__( + self, + samples: list[DiscreteVariable], + log_weights: np.ndarray, + exec: Executor = Executor() + ): + super().__init__(self, samples) + self.log_weights = log_weights + self.ess = ess(log_weights, exec=exec) From af5f945949fe7eb9b90bcd7dfe8a0f5880746257 Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 12:44:58 +0100 Subject: [PATCH 02/10] feat: change output of MCMC and SMC Change MCMC and SMC to output MCMCOutput and SMCOutput respectively --- discretesampling/base/algorithms/mcmc.py | 20 ++++++++++++++++++-- discretesampling/base/algorithms/smc.py | 4 +++- 2 files changed, 21 insertions(+), 3 deletions(-) diff --git a/discretesampling/base/algorithms/mcmc.py b/discretesampling/base/algorithms/mcmc.py index 87d8678..43382ac 100644 --- a/discretesampling/base/algorithms/mcmc.py +++ b/discretesampling/base/algorithms/mcmc.py @@ -1,7 +1,12 @@ import math import copy +<<<<<<< HEAD from tqdm.auto import tqdm +======= +import numpy as np +>>>>>>> d8b1383 (feat: change output of MCMC and SMC) from discretesampling.base.random import RNG +from discretesampling.base.output import MCMCOutput class DiscreteVariableMCMC(): @@ -12,8 +17,12 @@ def __init__(self, variableType, target, initialProposal): self.initialProposal = initialProposal self.target = target - def sample(self, N, seed=0, verbose=True): + def sample(self, N, N_warmup=None, seed=0, include_warmup=True, verbose = True): rng = RNG(seed) + + if N_warmup is None: + N_warmup = int(N/2) + initialSample = self.initialProposal.sample(rng) current = initialSample @@ -21,6 +30,7 @@ def sample(self, N, seed=0, verbose=True): display_progress_bar = verbose progress_bar = tqdm(total=N, desc="MCMC sampling", disable=not display_progress_bar) + acceptance_probabilities = [] for i in range(N): forward_proposal = self.proposalType(current, rng) @@ -49,7 +59,13 @@ def sample(self, N, seed=0, verbose=True): pass samples.append(copy.copy(current)) + acceptance_probabilities.append(acceptance_probability) progress_bar.update(1) + if not include_warmup: + samples = samples[(N_warmup-1):N] + acceptance_probabilities = acceptance_probabilities[(N_warmup-1):N] + + results = MCMCOutput(samples, acceptance_probabilities, include_warmup, N, N_warmup) progress_bar.close() - return samples + return results diff --git a/discretesampling/base/algorithms/smc.py b/discretesampling/base/algorithms/smc.py index 9d61493..ef2b3c8 100644 --- a/discretesampling/base/algorithms/smc.py +++ b/discretesampling/base/algorithms/smc.py @@ -7,6 +7,7 @@ from discretesampling.base.algorithms.smc_components.normalisation import normalise from discretesampling.base.algorithms.smc_components.effective_sample_size import ess from discretesampling.base.algorithms.smc_components.resampling import systematic_resampling +from discretesampling.base.output import SMCOutput class DiscreteVariableSMC(): @@ -79,5 +80,6 @@ def sample(self, Tsmc, N, seed=0, verbose=True): current_particles = new_particles progress_bar.update(1) + results = SMCOutput(current_particles, logWeights, exec=exec) progress_bar.close() - return current_particles + return results From 59f1de03d53b24d5d6b3cac933645172b5d60d47 Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 17:50:14 +0100 Subject: [PATCH 03/10] feat: results output for MPI Gather particles (and weights) prior to contructing SMCOutput block --- discretesampling/base/algorithms/mcmc.py | 9 +++------ discretesampling/base/algorithms/smc.py | 4 +++- .../algorithms/smc_components/effective_sample_size.py | 3 ++- discretesampling/base/executor/executor.py | 3 +++ discretesampling/base/executor/executor_MPI.py | 5 +++++ discretesampling/base/output/base.py | 4 ++-- discretesampling/base/output/mcmc_output.py | 5 ++--- discretesampling/base/output/smc_output.py | 8 +++----- examples/additive_structure_example.py | 6 +++--- 9 files changed, 26 insertions(+), 21 deletions(-) diff --git a/discretesampling/base/algorithms/mcmc.py b/discretesampling/base/algorithms/mcmc.py index 43382ac..622e41d 100644 --- a/discretesampling/base/algorithms/mcmc.py +++ b/discretesampling/base/algorithms/mcmc.py @@ -1,10 +1,7 @@ import math import copy -<<<<<<< HEAD from tqdm.auto import tqdm -======= import numpy as np ->>>>>>> d8b1383 (feat: change output of MCMC and SMC) from discretesampling.base.random import RNG from discretesampling.base.output import MCMCOutput @@ -17,7 +14,7 @@ def __init__(self, variableType, target, initialProposal): self.initialProposal = initialProposal self.target = target - def sample(self, N, N_warmup=None, seed=0, include_warmup=True, verbose = True): + def sample(self, N, N_warmup=None, seed=0, include_warmup=True, verbose=True): rng = RNG(seed) if N_warmup is None: @@ -65,7 +62,7 @@ def sample(self, N, N_warmup=None, seed=0, include_warmup=True, verbose = True): if not include_warmup: samples = samples[(N_warmup-1):N] acceptance_probabilities = acceptance_probabilities[(N_warmup-1):N] - + results = MCMCOutput(samples, acceptance_probabilities, include_warmup, N, N_warmup) progress_bar.close() - return results + return diff --git a/discretesampling/base/algorithms/smc.py b/discretesampling/base/algorithms/smc.py index ef2b3c8..d045be0 100644 --- a/discretesampling/base/algorithms/smc.py +++ b/discretesampling/base/algorithms/smc.py @@ -80,6 +80,8 @@ def sample(self, Tsmc, N, seed=0, verbose=True): current_particles = new_particles progress_bar.update(1) - results = SMCOutput(current_particles, logWeights, exec=exec) + current_particles = self.exec.gather_all(current_particles) + logWeights = self.exec.gather(logWeights, [N]) + results = SMCOutput(current_particles, logWeights) progress_bar.close() return results diff --git a/discretesampling/base/algorithms/smc_components/effective_sample_size.py b/discretesampling/base/algorithms/smc_components/effective_sample_size.py index 01c64f4..3645e34 100755 --- a/discretesampling/base/algorithms/smc_components/effective_sample_size.py +++ b/discretesampling/base/algorithms/smc_components/effective_sample_size.py @@ -1,7 +1,8 @@ import numpy as np +from discretesampling.base.executor import Executor -def ess(logw, exec): +def ess(logw, exec: Executor = Executor()): """ Description ----------- diff --git a/discretesampling/base/executor/executor.py b/discretesampling/base/executor/executor.py index 085359b..ca790e5 100644 --- a/discretesampling/base/executor/executor.py +++ b/discretesampling/base/executor/executor.py @@ -31,3 +31,6 @@ def redistribute(self, particles, ncopies): [[particles[i]]*ncopies[i] for i in range(len(particles))] )) return particles + + def gather_all(self, particles): + return particles diff --git a/discretesampling/base/executor/executor_MPI.py b/discretesampling/base/executor/executor_MPI.py index 852d024..bacf390 100644 --- a/discretesampling/base/executor/executor_MPI.py +++ b/discretesampling/base/executor/executor_MPI.py @@ -2,6 +2,7 @@ import numpy as np from scipy.special import logsumexp from discretesampling.base.executor import Executor +from discretesampling.base.util import gather_all from discretesampling.base.executor.MPI.distributed_fixed_size_redistribution.prefix_sum import ( inclusive_prefix_sum ) @@ -61,3 +62,7 @@ def cumsum(self, x): def redistribute(self, particles, ncopies): return variable_size_redistribution(particles, ncopies, self) + + def gather_all(self, particles): + particles = gather_all(particles, exec=self) + return particles diff --git a/discretesampling/base/output/base.py b/discretesampling/base/output/base.py index cd9ab98..7fd8085 100644 --- a/discretesampling/base/output/base.py +++ b/discretesampling/base/output/base.py @@ -1,8 +1,8 @@ -from abc import abc +from abc import ABC import numpy as np -class BaseOutput(abc): +class BaseOutput(ABC): def __init__(self, samples): self.samples = samples self.log_weights = np.repeat([-len(samples)], len(samples)) diff --git a/discretesampling/base/output/mcmc_output.py b/discretesampling/base/output/mcmc_output.py index 14b8eb7..1d42539 100644 --- a/discretesampling/base/output/mcmc_output.py +++ b/discretesampling/base/output/mcmc_output.py @@ -10,10 +10,9 @@ def __init__( samples: list[DiscreteVariable], acceptance_probabilities, include_warmup, - N=None, N_warmup=None, - exec: Executor = Executor() + N=None, N_warmup=None ): - super().__init__(self, samples) + super().__init__(samples) self.acceptance_probabilities = acceptance_probabilities self.include_warmup = include_warmup if N is None: diff --git a/discretesampling/base/output/smc_output.py b/discretesampling/base/output/smc_output.py index 69883fc..9cbadd0 100644 --- a/discretesampling/base/output/smc_output.py +++ b/discretesampling/base/output/smc_output.py @@ -2,16 +2,14 @@ from discretesampling.base.types import DiscreteVariable from discretesampling.base.output import BaseOutput from discretesampling.base.algorithms.smc_components import ess -from discretesampling.base.executor import Executor class SMCOutput(BaseOutput): def __init__( self, samples: list[DiscreteVariable], - log_weights: np.ndarray, - exec: Executor = Executor() + log_weights: np.ndarray ): - super().__init__(self, samples) + super().__init__(samples) self.log_weights = log_weights - self.ess = ess(log_weights, exec=exec) + self.ess = ess(log_weights) diff --git a/examples/additive_structure_example.py b/examples/additive_structure_example.py index b8fa1de..8316e47 100644 --- a/examples/additive_structure_example.py +++ b/examples/additive_structure_example.py @@ -37,9 +37,9 @@ def func(x): # Comparison between the three methods based on sampling the true structure -mcmc = [x.discrete_set for x in asSamplesMCMC] -lkern = [x.discrete_set for x in asSamplesLk] -smc = [x.discrete_set for x in asSamplesSMC] +mcmc = [x.discrete_set for x in asSamplesMCMC.samples] +lkern = [x.discrete_set for x in asSamplesLk.samples] +smc = [x.discrete_set for x in asSamplesSMC.samples] mcmc.count(true_str) lkern.count(true_str) From de18a893cd60a273266387e80d3362610ef1e23a Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 17:52:00 +0100 Subject: [PATCH 04/10] chore: update examples to match output format --- examples/decision_tree_example.py | 4 ++-- examples/mcmc_example.py | 2 +- examples/mpi/decision_tree_example.py | 5 +---- examples/mpi/regression_tree_example.py | 4 +--- examples/regression_decision_tree_example.py | 4 ++-- examples/smc_example.py | 2 +- 6 files changed, 8 insertions(+), 13 deletions(-) diff --git a/examples/decision_tree_example.py b/examples/decision_tree_example.py index 97760fc..5f638d9 100644 --- a/examples/decision_tree_example.py +++ b/examples/decision_tree_example.py @@ -29,7 +29,7 @@ try: treeSamples = dtMCMC.sample(500) - mcmcLabels = dt.stats(treeSamples, X_test).predict(X_test, use_majority=True) + mcmcLabels = dt.stats(treeSamples.samples, X_test).predict(X_test, use_majority=True) mcmcAccuracy = [dt.accuracy(y_test, mcmcLabels)] print("MCMC mean accuracy: ", (mcmcAccuracy)) except ZeroDivisionError: @@ -40,7 +40,7 @@ try: treeSMCSamples = dtSMC.sample(10, 1000) - smcLabels = dt.stats(treeSMCSamples, X_test).predict(X_test, use_majority=True) + smcLabels = dt.stats(treeSMCSamples.samples, X_test).predict(X_test, use_majority=True) smcAccuracy = [dt.accuracy(y_test, smcLabels)] print("SMC mean accuracy: ", np.mean(smcAccuracy)) diff --git a/examples/mcmc_example.py b/examples/mcmc_example.py index 48515c7..b79ce3a 100644 --- a/examples/mcmc_example.py +++ b/examples/mcmc_example.py @@ -24,7 +24,7 @@ try: treeSamples = dtMCMC.sample(N=1000) - mcmcLabels = dt.stats(treeSamples[500:999], X_test).predict(X_test) + mcmcLabels = dt.stats(treeSamples.samples[500:999], X_test).predict(X_test) mcmc_acc = dt.accuracy(y_test, mcmcLabels) print(numpy.mean(mcmc_acc)) except ZeroDivisionError: diff --git a/examples/mpi/decision_tree_example.py b/examples/mpi/decision_tree_example.py index 267dee5..a72a673 100644 --- a/examples/mpi/decision_tree_example.py +++ b/examples/mpi/decision_tree_example.py @@ -35,10 +35,7 @@ MPI.COMM_WORLD.Barrier() end = MPI.Wtime() - if MPI.COMM_WORLD.Get_size() > 1: - treeSMCSamples = gather_all(treeSMCSamples, exec) - - smcLabels = dt.stats(treeSMCSamples, X_test).predict(X_test) + smcLabels = dt.stats(treeSMCSamples.samples, X_test).predict(X_test) smcAccuracy = dt.accuracy(y_test, smcLabels) if MPI.COMM_WORLD.Get_rank() == 0: diff --git a/examples/mpi/regression_tree_example.py b/examples/mpi/regression_tree_example.py index 408a975..157b60b 100644 --- a/examples/mpi/regression_tree_example.py +++ b/examples/mpi/regression_tree_example.py @@ -32,10 +32,8 @@ treeSMCSamples = dtSMC.sample(T, N, seed) MPI.COMM_WORLD.Barrier() end = MPI.Wtime() - if MPI.COMM_WORLD.Get_size() > 1: - treeSMCSamples = gather_all(treeSMCSamples, exec) - smcLabels = dt.RegressionStats(treeSMCSamples, X_test).predict(X_test) + smcLabels = dt.RegressionStats(treeSMCSamples.samples, X_test).predict(X_test) smcAccuracy = dt.accuracy_mse(y_test, smcLabels) if MPI.COMM_WORLD.Get_rank() == 0: diff --git a/examples/regression_decision_tree_example.py b/examples/regression_decision_tree_example.py index 0091db7..6f50a60 100644 --- a/examples/regression_decision_tree_example.py +++ b/examples/regression_decision_tree_example.py @@ -19,7 +19,7 @@ try: treeSamples = dtMCMC.sample(100) - mcmcLabels = dt.RegressionStats(treeSamples, X_test).predict(X_test, use_majority=True) + mcmcLabels = dt.RegressionStats(treeSamples.samples, X_test).predict(X_test, use_majority=True) mcmcAccuracy = [dt.accuracy_mse(y_test, mcmcLabels)] print("MCMC mean MSE: ", (mcmcAccuracy)) @@ -31,7 +31,7 @@ try: treeSMCSamples = dtSMC.sample(10, 10) - smcLabels = dt.RegressionStats(treeSMCSamples, X_test).predict(X_test, use_majority=True) + smcLabels = dt.RegressionStats(treeSMCSamples.samples, X_test).predict(X_test, use_majority=True) smcAccuracy = [dt.accuracy_mse(y_test, smcLabels)] print("SMC mean MSE: ", (smcAccuracy)) diff --git a/examples/smc_example.py b/examples/smc_example.py index 6b2e20c..b392d7b 100644 --- a/examples/smc_example.py +++ b/examples/smc_example.py @@ -25,7 +25,7 @@ try: treeSamples = dtSMC.sample(10, 10) - smcLabels = dt.stats(treeSamples, X_test).predict(X_test, use_majority=True) + smcLabels = dt.stats(treeSamples.samples, X_test).predict(X_test, use_majority=True) smc_acc = dt.accuracy(y_test, smcLabels) print(numpy.mean(smc_acc)) except ZeroDivisionError: From 6eb62b6a213bf837e0715f2b7d72024cd9c47a02 Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 17:52:38 +0100 Subject: [PATCH 05/10] chore: remove unused imports --- discretesampling/base/output/mcmc_output.py | 1 - examples/mpi/decision_tree_example.py | 1 - examples/mpi/regression_tree_example.py | 1 - 3 files changed, 3 deletions(-) diff --git a/discretesampling/base/output/mcmc_output.py b/discretesampling/base/output/mcmc_output.py index 1d42539..187062f 100644 --- a/discretesampling/base/output/mcmc_output.py +++ b/discretesampling/base/output/mcmc_output.py @@ -1,7 +1,6 @@ import numpy as np from discretesampling.base.types import DiscreteVariable from discretesampling.base.output import BaseOutput -from discretesampling.base.executor import Executor class MCMCOutput(BaseOutput): diff --git a/examples/mpi/decision_tree_example.py b/examples/mpi/decision_tree_example.py index a72a673..61b53f9 100644 --- a/examples/mpi/decision_tree_example.py +++ b/examples/mpi/decision_tree_example.py @@ -5,7 +5,6 @@ from discretesampling.base.algorithms import DiscreteVariableSMC from discretesampling.domain import decision_tree as dt from discretesampling.base.executor.executor_MPI import Executor_MPI -from discretesampling.base.util import gather_all data = datasets.load_wine() diff --git a/examples/mpi/regression_tree_example.py b/examples/mpi/regression_tree_example.py index 157b60b..aceb0b7 100644 --- a/examples/mpi/regression_tree_example.py +++ b/examples/mpi/regression_tree_example.py @@ -5,7 +5,6 @@ from discretesampling.base.algorithms import DiscreteVariableSMC from discretesampling.domain import decision_tree as dt from discretesampling.base.executor.executor_MPI import Executor_MPI -from discretesampling.base.util import gather_all data = datasets.load_diabetes() From b11c17f183bbb978f06b408957074f56c262deea Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 18:43:10 +0100 Subject: [PATCH 06/10] fix: add __eq__ for BaseOutput Compare samples and weights for equality Needed for unit/regression testing --- discretesampling/base/output/base.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/discretesampling/base/output/base.py b/discretesampling/base/output/base.py index 7fd8085..3884f38 100644 --- a/discretesampling/base/output/base.py +++ b/discretesampling/base/output/base.py @@ -6,3 +6,8 @@ class BaseOutput(ABC): def __init__(self, samples): self.samples = samples self.log_weights = np.repeat([-len(samples)], len(samples)) + + def __eq__(self, other): + samples_eq = np.array_equal(self.samples, other.samples) + weights_eq = np.array_equal(self.log_weights, other.log_weights) + return samples_eq and weights_eq From 28c013a5ebb2e98410a3a7e528ad62084300d699 Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 18:43:54 +0100 Subject: [PATCH 07/10] chore: improved default for Outputs If values not given, assume equal weights --- discretesampling/base/output/mcmc_output.py | 9 ++++++--- discretesampling/base/output/smc_output.py | 8 +++++--- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/discretesampling/base/output/mcmc_output.py b/discretesampling/base/output/mcmc_output.py index 187062f..9f93530 100644 --- a/discretesampling/base/output/mcmc_output.py +++ b/discretesampling/base/output/mcmc_output.py @@ -7,8 +7,8 @@ class MCMCOutput(BaseOutput): def __init__( self, samples: list[DiscreteVariable], - acceptance_probabilities, - include_warmup, + acceptance_probabilities=None, + include_warmup=None, N=None, N_warmup=None ): super().__init__(samples) @@ -18,4 +18,7 @@ def __init__( N = len(samples) self.N = N self.N_warmup = N_warmup - self.acceptance_rate = np.mean(acceptance_probabilities) + acceptance_rate = None + if acceptance_probabilities is not None: + acceptance_rate = np.mean(acceptance_probabilities) + self.acceptance_rate = acceptance_rate diff --git a/discretesampling/base/output/smc_output.py b/discretesampling/base/output/smc_output.py index 9cbadd0..8b85f16 100644 --- a/discretesampling/base/output/smc_output.py +++ b/discretesampling/base/output/smc_output.py @@ -1,4 +1,5 @@ import numpy as np +from typing import Union from discretesampling.base.types import DiscreteVariable from discretesampling.base.output import BaseOutput from discretesampling.base.algorithms.smc_components import ess @@ -8,8 +9,9 @@ class SMCOutput(BaseOutput): def __init__( self, samples: list[DiscreteVariable], - log_weights: np.ndarray + log_weights: Union[np.ndarray, None] = None ): super().__init__(samples) - self.log_weights = log_weights - self.ess = ess(log_weights) + if log_weights is not None: + self.log_weights = log_weights + self.ess = ess(self.log_weights) From 7ecf77bb6dcf0c33348fef4f2095bd34c709cc25 Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 18:44:16 +0100 Subject: [PATCH 08/10] test: update SMC and MCMC tests to match --- tests/test_mcmc.py | 10 ++++++++-- tests/test_smc.py | 41 +++++++++++++++++++++++++++++++++++------ 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/tests/test_mcmc.py b/tests/test_mcmc.py index 3fc21de..1f16792 100644 --- a/tests/test_mcmc.py +++ b/tests/test_mcmc.py @@ -1,12 +1,18 @@ import pytest from discretesampling.base.algorithms import DiscreteVariableMCMC from discretesampling.domain import spectrum +from discretesampling.base.output import MCMCOutput @pytest.mark.parametrize( "seed,T,expected", - [(0, 10, [spectrum.SpectrumDimension(i) for i in [31, 30, 30, 30, 30, 29, 28, 27, 28, 27]]), - (1, 10, [spectrum.SpectrumDimension(i) for i in [27, 28, 27, 26, 25, 26, 27, 26, 25, 24]])] + [ + (0, 10, MCMCOutput( + [spectrum.SpectrumDimension(i) for i in [31, 30, 30, 30, 30, 29, 28, 27, 28, 27]] + )), + (1, 10, MCMCOutput( + [spectrum.SpectrumDimension(i) for i in [27, 28, 27, 26, 25, 26, 27, 26, 25, 24]])) + ] ) def test_mcmc(seed, T, expected): diff --git a/tests/test_smc.py b/tests/test_smc.py index e576bc0..c98298d 100644 --- a/tests/test_smc.py +++ b/tests/test_smc.py @@ -1,16 +1,32 @@ import pytest +import numpy as np from discretesampling.base.algorithms import DiscreteVariableSMC from discretesampling.base.util import gather_all from discretesampling.base.executor.executor_MPI import Executor_MPI from discretesampling.domain import spectrum +from discretesampling.base.output import SMCOutput # Test reproducibility @pytest.mark.parametrize( "seed,T,N,expected", - [(0, 3, 16, [spectrum.SpectrumDimension(i) for i in [15, 13, 15, 6, 6, 14, 8, 8, 8, 6, 6, 6, 14, 14, 12, 10]]), - (1, 3, 16, [spectrum.SpectrumDimension(i) for i in [13, 15, 15, 6, 2, 6, 18, 8, 6, 6, 6, 8, 14, 12, 10, 12]])] + [ + (0, 3, 16, SMCOutput( + [spectrum.SpectrumDimension(i) for i in [15, 13, 15, 6, 6, 14, 8, 8, 8, 6, 6, 6, 14, 14, 12, 10]], + np.array([-3.3423949417970995, -2.6927167512983843, -3.3423949417970995, -2.5485352600349263, -2.5485352600349263, + -1.7791312350561306, -2.807685770214025, -2.807685770214025, -2.807685770214025, -3.2906867661536747, + -3.2906867661536747, -3.2906867661536747, -3.286288420776657, -3.286288420776657, -2.753770144750778, + -2.4895429738738315]) + )), + (1, 3, 16, SMCOutput( + [spectrum.SpectrumDimension(i) for i in [13, 15, 15, 6, 2, 6, 18, 8, 6, 6, 6, 8, 14, 12, 10, 12]], + np.array([-2.577088421206954, -3.2267666117056693, -3.2267666117056693, -2.432906929943496, -5.389379226012296, + -2.432906929943496, -3.37139647864896, -2.692057440122595, -3.1750584360622445, -3.1750584360622445, + -3.1750584360622445, -2.692057440122595, -3.1706600906852267, -2.6381418146593476, -2.3739146437824012, + -2.6381418146593476]) + ) + )] ) def test_smc(seed, T, N, expected): target = spectrum.SpectrumDimensionTarget(10, 3.4) # NB with mean 10 and variance 3.4^2 @@ -25,8 +41,22 @@ def test_smc(seed, T, N, expected): @pytest.mark.mpi @pytest.mark.parametrize( "seed,T,N,expected", - [(0, 3, 16, [spectrum.SpectrumDimension(i) for i in [15, 13, 15, 6, 6, 14, 8, 8, 8, 6, 6, 6, 14, 14, 12, 10]]), - (1, 3, 16, [spectrum.SpectrumDimension(i) for i in [13, 15, 15, 6, 2, 6, 18, 8, 6, 6, 6, 8, 14, 12, 10, 12]])] + [ + (0, 3, 16, SMCOutput( + [spectrum.SpectrumDimension(i) for i in [15, 13, 15, 6, 6, 14, 8, 8, 8, 6, 6, 6, 14, 14, 12, 10]], + np.array([-3.3423949417970995, -2.6927167512983843, -3.3423949417970995, -2.5485352600349263, -2.5485352600349263, + -1.7791312350561306, -2.807685770214025, -2.807685770214025, -2.807685770214025, -3.2906867661536747, + -3.2906867661536747, -3.2906867661536747, -3.286288420776657, -3.286288420776657, -2.753770144750778, + -2.4895429738738315]) + )), + (1, 3, 16, SMCOutput( + [spectrum.SpectrumDimension(i) for i in [13, 15, 15, 6, 2, 6, 18, 8, 6, 6, 6, 8, 14, 12, 10, 12]], + np.array([-2.577088421206954, -3.2267666117056693, -3.2267666117056693, -2.432906929943496, -5.389379226012296, + -2.432906929943496, -3.37139647864896, -2.692057440122595, -3.1750584360622445, -3.1750584360622445, + -3.1750584360622445, -2.692057440122595, -3.1706600906852267, -2.6381418146593476, -2.3739146437824012, + -2.6381418146593476]) + )) + ] ) def test_smc_MPI(seed, T, N, expected): target = spectrum.SpectrumDimensionTarget(10, 3.4) # NB with mean 10 and variance 3.4^2 @@ -34,6 +64,5 @@ def test_smc_MPI(seed, T, N, expected): exec = Executor_MPI() specSMC = DiscreteVariableSMC(spectrum.SpectrumDimension, target, initialProposal, exec=exec) - local_samples = specSMC.sample(T, N, seed=seed) - samples = gather_all(local_samples, exec) + samples = specSMC.sample(T, N, seed=seed) assert samples == expected From 69955a81f6d6134d21a026652b00c5ec37970b14 Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Wed, 27 Sep 2023 18:47:12 +0100 Subject: [PATCH 09/10] style: fix linting errors --- discretesampling/base/algorithms/mcmc.py | 3 +-- discretesampling/base/output/__init__.py | 6 +++--- tests/test_smc.py | 1 - 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/discretesampling/base/algorithms/mcmc.py b/discretesampling/base/algorithms/mcmc.py index 622e41d..09cf752 100644 --- a/discretesampling/base/algorithms/mcmc.py +++ b/discretesampling/base/algorithms/mcmc.py @@ -1,7 +1,6 @@ import math import copy from tqdm.auto import tqdm -import numpy as np from discretesampling.base.random import RNG from discretesampling.base.output import MCMCOutput @@ -65,4 +64,4 @@ def sample(self, N, N_warmup=None, seed=0, include_warmup=True, verbose=True): results = MCMCOutput(samples, acceptance_probabilities, include_warmup, N, N_warmup) progress_bar.close() - return + return results diff --git a/discretesampling/base/output/__init__.py b/discretesampling/base/output/__init__.py index 1a1a943..bef45fd 100644 --- a/discretesampling/base/output/__init__.py +++ b/discretesampling/base/output/__init__.py @@ -1,3 +1,3 @@ -from .base import BaseOutput # no qa -from .smc_output import SMCOutput # no qa -from .mcmc_output import MCMCOutput # no qa +from .base import BaseOutput # noqa +from .smc_output import SMCOutput # noqa +from .mcmc_output import MCMCOutput # noqa diff --git a/tests/test_smc.py b/tests/test_smc.py index c98298d..141eb3f 100644 --- a/tests/test_smc.py +++ b/tests/test_smc.py @@ -1,7 +1,6 @@ import pytest import numpy as np from discretesampling.base.algorithms import DiscreteVariableSMC -from discretesampling.base.util import gather_all from discretesampling.base.executor.executor_MPI import Executor_MPI from discretesampling.domain import spectrum from discretesampling.base.output import SMCOutput From ecfd2ffc2a9580c9ce94ef4255f2025eea75e29f Mon Sep 17 00:00:00 2001 From: Alex Phillips Date: Thu, 28 Sep 2023 15:54:54 +0100 Subject: [PATCH 10/10] feat: add flag to disable gathering of particles in smc --- discretesampling/base/algorithms/smc.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/discretesampling/base/algorithms/smc.py b/discretesampling/base/algorithms/smc.py index d045be0..52cf094 100644 --- a/discretesampling/base/algorithms/smc.py +++ b/discretesampling/base/algorithms/smc.py @@ -30,7 +30,7 @@ def __init__(self, variableType, target, initialProposal, self.initialProposal = initialProposal self.target = target - def sample(self, Tsmc, N, seed=0, verbose=True): + def sample(self, Tsmc, N, seed=0, gather_results=True, verbose=True): loc_n = int(N/self.exec.P) rank = self.exec.rank mvrs_rng = RNG(seed) @@ -80,8 +80,10 @@ def sample(self, Tsmc, N, seed=0, verbose=True): current_particles = new_particles progress_bar.update(1) - current_particles = self.exec.gather_all(current_particles) - logWeights = self.exec.gather(logWeights, [N]) + if gather_results: + current_particles = self.exec.gather_all(current_particles) + logWeights = self.exec.gather(logWeights, [N]) + results = SMCOutput(current_particles, logWeights) progress_bar.close() return results