From eefd9e649ace32cd1e6dbc8c71e828ade5d57baa Mon Sep 17 00:00:00 2001 From: Jasper Everink Date: Tue, 24 Jun 2025 13:37:47 +0200 Subject: [PATCH 1/3] Fix initial_point for JointDistribution --- cuqi/experimental/mcmc/_cwmh.py | 7 ++----- cuqi/experimental/mcmc/_mh.py | 3 ++- cuqi/experimental/mcmc/_sampler.py | 2 ++ 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/cuqi/experimental/mcmc/_cwmh.py b/cuqi/experimental/mcmc/_cwmh.py index 84ee1e1068..9238eaf00a 100644 --- a/cuqi/experimental/mcmc/_cwmh.py +++ b/cuqi/experimental/mcmc/_cwmh.py @@ -93,12 +93,9 @@ def scale(self, value): self._scale = value def validate_target(self): - if not isinstance(self.target, cuqi.density.Density): - raise ValueError( - "Target should be an instance of "+\ - f"{cuqi.density.Density.__class__.__name__}") # Fail when there is no log density, which is currently assumed to be the case in case NaN is returned. - if np.isnan(self.target.logd(self._get_default_initial_point(self.dim))): + initial_point = self._get_default_initial_point(self.dim) + if (isinstance(initial_point, list) and np.isnan(self.target.logd(*initial_point))) or (not isinstance(initial_point, list) and np.isnan(self.target.logd(initial_point))): raise ValueError("Target does not have valid logd") def validate_proposal(self): diff --git a/cuqi/experimental/mcmc/_mh.py b/cuqi/experimental/mcmc/_mh.py index 1322d9b7dd..923fdea9a8 100644 --- a/cuqi/experimental/mcmc/_mh.py +++ b/cuqi/experimental/mcmc/_mh.py @@ -33,7 +33,8 @@ def _initialize(self): def validate_target(self): # Fail only when there is no log density, which is currently assumed to be the case in case NaN is returned. - if np.isnan(self.target.logd(self._get_default_initial_point(self.dim))): + initial_point = self._get_default_initial_point(self.dim) + if (isinstance(initial_point, list) and np.isnan(self.target.logd(*initial_point))) or (not isinstance(initial_point, list) and np.isnan(self.target.logd(initial_point))): raise ValueError("Target does not have valid logd") def validate_proposal(self): diff --git a/cuqi/experimental/mcmc/_sampler.py b/cuqi/experimental/mcmc/_sampler.py index 90e11f2eeb..21c69e150a 100644 --- a/cuqi/experimental/mcmc/_sampler.py +++ b/cuqi/experimental/mcmc/_sampler.py @@ -388,6 +388,8 @@ def _ensure_initialized(self): def _get_default_initial_point(self, dim): """ Return the default initial point for the sampler. Defaults to an array of ones. """ + if isinstance(dim, list): + return [np.ones(d) for d in dim] return np.ones(dim) def __repr__(self): From cc6e8360e20f88a1e2b0a2a16d18ff3288bef0c8 Mon Sep 17 00:00:00 2001 From: Jasper Everink Date: Wed, 25 Jun 2025 10:33:45 +0200 Subject: [PATCH 2/3] Use StackedJointDistribution in ProposalBasedSampler --- cuqi/distribution/_joint_distribution.py | 30 ++++++++++++++++-------- cuqi/experimental/mcmc/_cwmh.py | 3 +-- cuqi/experimental/mcmc/_mh.py | 3 +-- cuqi/experimental/mcmc/_sampler.py | 20 ++++++++++++++++ 4 files changed, 42 insertions(+), 14 deletions(-) diff --git a/cuqi/distribution/_joint_distribution.py b/cuqi/distribution/_joint_distribution.py index 682e962735..e952138108 100644 --- a/cuqi/distribution/_joint_distribution.py +++ b/cuqi/distribution/_joint_distribution.py @@ -5,6 +5,7 @@ from cuqi.distribution import Distribution, Posterior from cuqi.likelihood import Likelihood from cuqi.geometry import Geometry, _DefaultGeometry1D +from cuqi.samples import Samples, JointSamples import cuqi import numpy as np # for splitting array. Can avoid. @@ -311,16 +312,7 @@ def geometry(self): def logd(self, stacked_input): """ Return the un-normalized log density function stacked joint density. """ - - # Split the stacked input into individual inputs and call superclass - split_indices = np.cumsum(super().dim) # list(accumulate(super().dim)) - inputs = np.split(stacked_input, split_indices[:-1]) - names = self.get_parameter_names() - - # Create keyword arguments - kwargs = dict(zip(names, inputs)) - - return super().logd(**kwargs) + return super().logd(**self._unstack_input(stacked_input)) def logpdf(self, stacked_input): return self.logd(stacked_input) @@ -330,6 +322,24 @@ def _sample(self, Ns=1): def __repr__(self): return "_Stacked"+super().__repr__() + + def _unstack_input(self, stacked_input): + # Split the stacked input into individual inputs and call superclass + split_indices = np.cumsum(super().dim) # list(accumulate(super().dim)) + inputs = np.split(stacked_input, split_indices[:-1]) + names = self.get_parameter_names() + + # Create keyword arguments + return dict(zip(names, inputs)) + + def _unstack_samples(self, stacked_samples): + split_indices = np.cumsum(super().dim) + split_samples = np.split(stacked_samples.samples, split_indices[:-1], axis = 0) + names = self.get_parameter_names() + geometries = [dist.geometry for dist in self._distributions] + + unstacked_samples = [Samples(samples, geometry) for samples, geometry in zip(split_samples, geometries)] + return JointSamples(zip(names, unstacked_samples)) class MultipleLikelihoodPosterior(JointDistribution, Distribution): diff --git a/cuqi/experimental/mcmc/_cwmh.py b/cuqi/experimental/mcmc/_cwmh.py index 9238eaf00a..ae48f3748b 100644 --- a/cuqi/experimental/mcmc/_cwmh.py +++ b/cuqi/experimental/mcmc/_cwmh.py @@ -94,8 +94,7 @@ def scale(self, value): def validate_target(self): # Fail when there is no log density, which is currently assumed to be the case in case NaN is returned. - initial_point = self._get_default_initial_point(self.dim) - if (isinstance(initial_point, list) and np.isnan(self.target.logd(*initial_point))) or (not isinstance(initial_point, list) and np.isnan(self.target.logd(initial_point))): + if np.isnan(self.target.logd(self._get_default_initial_point(self.dim))): raise ValueError("Target does not have valid logd") def validate_proposal(self): diff --git a/cuqi/experimental/mcmc/_mh.py b/cuqi/experimental/mcmc/_mh.py index 923fdea9a8..1322d9b7dd 100644 --- a/cuqi/experimental/mcmc/_mh.py +++ b/cuqi/experimental/mcmc/_mh.py @@ -33,8 +33,7 @@ def _initialize(self): def validate_target(self): # Fail only when there is no log density, which is currently assumed to be the case in case NaN is returned. - initial_point = self._get_default_initial_point(self.dim) - if (isinstance(initial_point, list) and np.isnan(self.target.logd(*initial_point))) or (not isinstance(initial_point, list) and np.isnan(self.target.logd(initial_point))): + if np.isnan(self.target.logd(self._get_default_initial_point(self.dim))): raise ValueError("Target does not have valid logd") def validate_proposal(self): diff --git a/cuqi/experimental/mcmc/_sampler.py b/cuqi/experimental/mcmc/_sampler.py index 21c69e150a..51ea45dc73 100644 --- a/cuqi/experimental/mcmc/_sampler.py +++ b/cuqi/experimental/mcmc/_sampler.py @@ -477,6 +477,26 @@ def initialize(self): self._is_initialized = True + @Sampler.target.setter + def target(self, value): + """ Set the target density. Runs validation of the target. """ + self._target = value + + if isinstance(self._target, cuqi.distribution.JointDistribution): + self._target = self._target._as_stacked() + + if self._target is not None: + self.validate_target() + + def get_samples(self) -> Samples: + """ Return the samples. The internal data-structure for the samples is a dynamic list so this creates a copy. """ + + samples = Samples(np.array(self._samples).T, self.target.geometry) + if isinstance(self.target, cuqi.distribution.JointDistribution): + return self.target._unstack_samples(samples) + return Samples(samples) + + @abstractmethod def validate_proposal(self): """ Validate the proposal distribution. """ From 3092bf6abdbaca6a5f8bf9d6cc412d2bb9551850 Mon Sep 17 00:00:00 2001 From: Jasper Everink Date: Wed, 25 Jun 2025 10:59:54 +0200 Subject: [PATCH 3/3] Fix error --- cuqi/experimental/mcmc/_sampler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cuqi/experimental/mcmc/_sampler.py b/cuqi/experimental/mcmc/_sampler.py index 51ea45dc73..be110235bb 100644 --- a/cuqi/experimental/mcmc/_sampler.py +++ b/cuqi/experimental/mcmc/_sampler.py @@ -494,7 +494,7 @@ def get_samples(self) -> Samples: samples = Samples(np.array(self._samples).T, self.target.geometry) if isinstance(self.target, cuqi.distribution.JointDistribution): return self.target._unstack_samples(samples) - return Samples(samples) + return samples @abstractmethod