hops
ReversibleJumpProposal.hpp
Go to the documentation of this file.
1 #ifndef HOPS_REVERSIBLEJUMPPROPOSAL_HPP
2 #define HOPS_REVERSIBLEJUMPPROPOSAL_HPP
3 
4 #include <Eigen/Core>
5 #include <hops/MarkovChain/Proposal/ChordStepDistributions.hpp>
6 #include <hops/MarkovChain/Draw/IsCalculateLogAcceptanceProbabilityAvailable.hpp>
7 #include <hops/RandomNumberGenerator/RandomNumberGenerator.hpp>
8 #include <random>
9 #include <string>
10 #include <utility>
11 #include <vector>
12 
13 namespace {
14  std::pair<double, double>
15  distanceInCoordinateDirection(Eigen::MatrixXd A, Eigen::VectorXd b, Eigen::VectorXd x, long coordinate) {
16  Eigen::VectorXd slacks = b - A * x;
17  Eigen::VectorXd inverseDistances = A.col(coordinate).cwiseQuotient(slacks);
18  double forwardDistance = 1. / inverseDistances.maxCoeff();
19  double backwardDistance = 1. / inverseDistances.minCoeff();
20  assert(backwardDistance <= 0 && forwardDistance >= 0);
21  return std::make_pair(backwardDistance, forwardDistance);
22  }
23 }
24 
25 namespace hops {
26  template<typename MarkovChainImpl, typename Model>
27  class ReversibleJumpProposal : public MarkovChainImpl, public Model {
28  public:
29  ReversibleJumpProposal(const MarkovChainImpl &markovChainImpl,
30  const Model &model,
31  Eigen::VectorXi m_jumpIndices,
32  const typename MarkovChainImpl::StateType parameterDefaultValues) :
33  MarkovChainImpl(markovChainImpl),
34  Model(model),
35  jumpIndices_(std::move(m_jumpIndices)),
36  defaultParameterValues(parameterDefaultValues) {
37  for (long i = 0; i < parameterDefaultValues.rows(); i++) {
38  parameterActivationStates_.emplace_back(1);
39  }
40  // Starts with all optional parameters deactivated
41  typename MarkovChainImpl::StateType parameterState = MarkovChainImpl::getState();
42  for (long i = 0; i < jumpIndices_.rows(); ++i) {
43  parameterActivationStates_[jumpIndices_(i)] = 0;
44  parameterState[jumpIndices_(i)] = defaultParameterValues[jumpIndices_(i)];
45  }
46  MarkovChainImpl::setState(parameterState);
47  }
48 
49  void draw(RandomNumberGenerator &randomNumberGenerator) {
50  if (uniformRealDistribution(randomNumberGenerator) < m_modelJumpProbability) {
51  drawInModelSpace(randomNumberGenerator, parameterActivationStates_, defaultParameterValues);
52  } else {
53  drawInParameterSpace(randomNumberGenerator);
54  }
55  }
56 
57  void drawInModelSpace(RandomNumberGenerator &randomNumberGenerator,
58  std::vector<int> &parameterActivationStates,
59  typename MarkovChainImpl::StateType m_defaultValues) {
60  typename MarkovChainImpl::StateType m_proposal = MarkovChainImpl::getState();
61  std::vector<int> proposalActivationStates(parameterActivationStates_);
62 
63  double logModelJumpProbabilityDifferential = jumpModel(randomNumberGenerator,
64  proposalActivationStates,
65  m_proposal,
66  m_defaultValues);
67 
68  proposalNegativeLogLikelihood = Model::computeNegativeLogLikelihood(m_proposal);
69  double modelJumpAcceptanceProbability = stateNegativeLogLikelihood - proposalNegativeLogLikelihood +
70  logModelJumpProbabilityDifferential;
71  if (std::log(uniformRealDistribution(randomNumberGenerator)) <= modelJumpAcceptanceProbability) {
72  stateNegativeLogLikelihood = proposalNegativeLogLikelihood;
73  proposalActivationStates.swap(parameterActivationStates_);
74  MarkovChainImpl::setState(m_proposal);
75  } else {
76  // TODO update rejection stats
77  }
78  }
79 
80 
81  void drawInParameterSpace(RandomNumberGenerator &randomNumberGenerator) {
82  numberOfProposals++;
83  MarkovChainImpl::propose(randomNumberGenerator);
84  for(long i=0; i<MarkovChainImpl::m_proposal.rows(); ++i) {
85  if(!parameterActivationStates_[i]) {
86  // If parameter is not active, reset m_proposal to state
87  MarkovChainImpl::m_proposal(i) = MarkovChainImpl::state(i);
88  }
89  }
90  double acceptanceProbability = computeParameterDrawAcceptanceProbability();
91  double acceptanceChance = std::log(uniformRealDistribution(randomNumberGenerator));
92  if (acceptanceChance < acceptanceProbability) {
93  MarkovChainImpl::acceptProposal();
94  stateNegativeLogLikelihood = proposalNegativeLogLikelihood;
95  numberOfAcceptedProposals++;
96  }
97  }
98 
99  double computeParameterDrawAcceptanceProbability() {
100  double acceptanceProbability = 0;
101  if constexpr(IsCalculateLogAcceptanceProbabilityAvailable<MarkovChainImpl>::value) {
102  acceptanceProbability += MarkovChainImpl::computeLogAcceptanceProbability();
103  }
104  if (std::isfinite(acceptanceProbability)) {
105  proposalNegativeLogLikelihood = Model::computeNegativeLogLikelihood(
106  MarkovChainImpl::getProposal());
107  acceptanceProbability += stateNegativeLogLikelihood - proposalNegativeLogLikelihood;
108  }
109  return acceptanceProbability;
110  }
111 
112  double getAcceptanceRate() {
113  return static_cast<double>(numberOfAcceptedProposals) / numberOfProposals;
114  }
115 
116  typename MarkovChainImpl::StateType getState() {
117  typename MarkovChainImpl::StateType parameterState = MarkovChainImpl::getState();
118  typename MarkovChainImpl::StateType state(parameterState.rows() + 1);
119  long modelIndex = std::accumulate(parameterActivationStates_.begin(),
120  parameterActivationStates_.end(),
121  0,
122  [](unsigned int x, unsigned int y) { return (x << 1) + y; });
123  state << modelIndex, parameterState;
124  return state;
125  }
126 
127  [[nodiscard]] std::vector<std::string> getParameterNames() const {
128  // Vector is constructed on demand, because it typically is not used repeatedly.
129  std::vector<std::string> parameterNames = Model::getParameterNames();
130  std::vector<std::string> names = {"model index"};
131  names.insert(names.end(), parameterNames.begin(), parameterNames.end());
132  return names;
133  }
134 
135  private:
136  double stateNegativeLogLikelihood = 0;
137  double proposalNegativeLogLikelihood = 0;
138 
139  // fixed value from https://doi.org/10.1093/bioinformatics/btz500
140  typename MarkovChainImpl::StateType::Scalar m_modelJumpProbability = 0.5;
141  typename MarkovChainImpl::StateType::Scalar parameterActivationProbability = 0.1;
142  typename MarkovChainImpl::StateType::Scalar parameterDeactivationProbability = 0.1;
143  typename MarkovChainImpl::StateType::Scalar stepSize = 0.1;
144  std::uniform_real_distribution<double> uniformRealDistribution;
145  hops::GaussianStepDistribution<double> gaussianStepDistribution;
146 
147  typename MarkovChainImpl::StateType defaultParameterValues;
148  Eigen::VectorXi jumpIndices_;
149 
150  std::vector<int> parameterActivationStates_;
151 
152  long numberOfAcceptedProposals = 0;
153  long numberOfProposals = 0;
154 
155 
156  double jumpModel(hops::RandomNumberGenerator &randomNumberGenerator,
157  std::vector<int> &proposalActivationState,
158  typename MarkovChainImpl::StateType &m_proposal,
159  typename MarkovChainImpl::StateType &m_defaultValues) {
160  Eigen::VectorXi activationTracker = Eigen::VectorXi::Zero(jumpIndices_.size());
161  Eigen::VectorXi deactivationTracker = Eigen::VectorXi::Zero(jumpIndices_.size());
162 
163  Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> permutationMatrix(jumpIndices_.size());
164  permutationMatrix.setIdentity();
165  std::shuffle(permutationMatrix.indices().data(),
166  permutationMatrix.indices().data() + permutationMatrix.indices().size(),
167  randomNumberGenerator);
168  Eigen::VectorXi shuffledJumpIndices = permutationMatrix * jumpIndices_;
169  auto A = Model::getA();
170  auto b = Model::getB();
171  double j_fwd = 1;
172  double u_fwd = 1;
173  std::vector<double> tester;
174  for (long index = 0; index < shuffledJumpIndices.size(); ++index) {
175  // If parameter is active, sample deactivation
176  double defaultValue = m_defaultValues[shuffledJumpIndices(index)];
177  if (proposalActivationState[shuffledJumpIndices(index)]) {
178  // NOTE, WE MEASURE DISTANCE FROM DEFAULT BOTH IN ACTIVATION AND DEACTIVATION
179  // FOR NON-SQUARE REGIONS, SWITCH OUT DEFAULT FOR EVERYTHING EXCEPT THE CURRENT INDEX
180  auto[backwardsDistance, forwardsDistance] = distanceInCoordinateDirection(
181  A,
182  b,
183  m_defaultValues,
184  shuffledJumpIndices(index));
185  double width = (forwardsDistance - backwardsDistance);
186  // FOR NOW WE USE GAUSSIANS FOR BOTH ACTIVATION AND DEACTIVATION
187  double temp = gaussianStepDistribution.computeProbabilityDensity(
188  m_proposal[shuffledJumpIndices(index)] - defaultValue, stepSize * width, backwardsDistance,
189  forwardsDistance);
190 
191  double procProb = std::min(temp * parameterDeactivationProbability, 1.);
192  if (uniformRealDistribution(randomNumberGenerator) < procProb) {
193  proposalActivationState[shuffledJumpIndices(index)] = 0;
194  m_proposal[shuffledJumpIndices(index)] = defaultValue;
195  deactivationTracker(shuffledJumpIndices(index)) = 1;
196  j_fwd *= procProb;
197  } else {
198  j_fwd *= 1 - procProb;
199  }
200  } else { // If parameter is inactive, sample activation
201  if (uniformRealDistribution(randomNumberGenerator) < parameterActivationProbability) {
202  proposalActivationState[shuffledJumpIndices(index)] = 1;
203  // SAME AS ABOVE
204  auto[backwardsDistance, forwardsDistance] = distanceInCoordinateDirection(
205  A,
206  b,
207  m_defaultValues,
208  shuffledJumpIndices(index));
209 
210  double width = (forwardsDistance - backwardsDistance);
211  //THIS IS A TEMPORARY CHECK, SPECIFIC FOR THE GAMMA EXAMPLE
212  assert((std::abs(width - 0.9) < 1e-12) || (std::abs(width - 9.9) < 1e-12));
213  // RESETTING THE BOUNDS OF THE GAUSSIAN STEP DISTRIBUTION IS A WASTE (DEPENDING ON GENERALITY), WE COULD INSTEAD SAVE
214  // ONE DISTRIBUTION PER JUMPSTATE
215  m_proposal[shuffledJumpIndices(index)] += gaussianStepDistribution.draw(randomNumberGenerator,
216  stepSize * width,
217  backwardsDistance,
218  forwardsDistance);
219  assert(defaultValue + backwardsDistance <= m_proposal[shuffledJumpIndices(index)] && forwardsDistance + defaultValue >= m_proposal[shuffledJumpIndices(index)]);
220  activationTracker(shuffledJumpIndices(index)) = 1;
221  // GET THE DENSITY OF THE TAKEN STEP
222  double temp = gaussianStepDistribution.computeProbabilityDensity(
223  m_proposal[shuffledJumpIndices(index)] - defaultValue, stepSize * width, backwardsDistance,
224  forwardsDistance);
225  tester.push_back(m_proposal[shuffledJumpIndices(index)] - defaultValue);
226  tester.push_back(stepSize * width);
227  tester.push_back(temp);
228  assert(temp > 0);
229  // HERE WE HAVE TO MULTIPLY WITH THE WIDTH TO COMPENSATE FOR THE PRIOR VOLUMES
230  u_fwd *= temp * width;
231  j_fwd *= parameterActivationProbability;
232  } else {
233  j_fwd *= 1 - parameterActivationProbability;
234  }
235  }
236  }
237 
238  double j_bck = 1;
239  double u_bck = 1;
240  std::shuffle(permutationMatrix.indices().data(),
241  permutationMatrix.indices().data() + permutationMatrix.indices().size(),
242  randomNumberGenerator);
243  shuffledJumpIndices = permutationMatrix * jumpIndices_;
244 
245  // computes reverse probability (swap deactivation with activation)
246  for (long index = 0; index < shuffledJumpIndices.size(); ++index) {
247  // Note the if now checks for deactivated states
248  if (!proposalActivationState[shuffledJumpIndices(index)]) {
249  if (deactivationTracker(shuffledJumpIndices(index))) {
250  // AGAIN WASTEFUL COMPUTATION OF CONSTANT
251  auto[backwardsDistance, forwardsDistance] = distanceInCoordinateDirection(
252  A,
253  b,
254  m_defaultValues,
255  shuffledJumpIndices(index));
256  j_bck *= parameterActivationProbability;
257  double defaultValue = m_defaultValues[shuffledJumpIndices(index)];
258  // NOW TRACK PROBABILITY FROM DEFAULT TO PREVIOUS STATE, GET BETTER WAY OF RETRIEVING PAST STATE WITHOUT + 1!!!!
259  auto temp_old_state = this->getState()[shuffledJumpIndices(index) + 1];
260  assert(defaultValue + backwardsDistance <= temp_old_state && forwardsDistance + defaultValue >= temp_old_state);
261  // also here we have to take the width into account
262  double width = (forwardsDistance - backwardsDistance);
263  //THIS IS A TEMPORARY CHECK, SPECIFIC FOR THE GAMMA EXAMPLE
264  assert((std::abs(width - 0.9) < 1e-12) || (std::abs(width - 9.9) < 1e-12));
265  double temp = gaussianStepDistribution.computeProbabilityDensity(
266  this->getState()[shuffledJumpIndices(index) + 1] - defaultValue, stepSize * width, backwardsDistance,
267  forwardsDistance);
268  //tester.push_back(temp);
269  if (temp <= 0){
270  temp += 0;
271  }
272  //assert(temp > 0);
273  u_bck *= temp * width;
274  } else {
275  j_bck *= 1 - parameterActivationProbability;
276  }
277  } else {
278  double defaultValue = defaultParameterValues[shuffledJumpIndices(index)];
279  //Eigen::VectorXd shiftedProposal = m_proposal;
280  //shiftedProposal(shuffledJumpIndices(index)) -= defaultValue;
281  //NOW CHECK PROBABILITY OF DEACTIVATING THE PROPOSAL
282  auto[backwardsDistance, forwardsDistance] = distanceInCoordinateDirection(
283  A,
284  b,
285  m_defaultValues,
286  shuffledJumpIndices(index));
287  double width = (forwardsDistance - backwardsDistance);
288  //THIS IS A TEMPORARY CHECK, SPECIFIC FOR THE GAMMA EXAMPLE
289  assert((std::abs(width - 0.9) < 1e-12) || (std::abs(width - 9.9) < 1e-12));
290  //double defaultDistance = std::abs(m_proposal[shuffledJumpIndices(index)] - defaultValue);
291  //defaultDistance = defaultDistance / width;
292 
293  //double prob = std::pow(defaultDistance, parameterDeactivationProbability);
294  double temp = gaussianStepDistribution.computeProbabilityDensity(
295  m_proposal[shuffledJumpIndices(index)] - defaultValue, stepSize * width, backwardsDistance,
296  forwardsDistance);
297  tester.push_back(m_proposal[shuffledJumpIndices(index)] - defaultValue);
298  tester.push_back(stepSize * width);
299  tester.push_back(temp);
300  double prob = std::min(temp * parameterDeactivationProbability, 1.);
301  if (activationTracker(shuffledJumpIndices(index))) {
302  j_bck *= prob;
303  } else {
304  j_bck *= 1 - prob;
305  }
306  }
307  }
308  return std::log(j_bck * u_bck / (j_fwd * u_fwd));
309  }
310  };
311 }
312 
313 #endif //HOPS_REVERSIBLEJUMPPROPOSAL_HPP
hops::ReversibleJumpProposal::computeParameterDrawAcceptanceProbability
double computeParameterDrawAcceptanceProbability()
Definition: ReversibleJumpProposal.hpp:99
RandomNumberGenerator.hpp
hops::ReversibleJumpProposal::drawInParameterSpace
void drawInParameterSpace(RandomNumberGenerator &randomNumberGenerator)
Definition: ReversibleJumpProposal.hpp:81
hops::Model::computeNegativeLogLikelihood
virtual MatrixType::Scalar computeNegativeLogLikelihood(const VectorType &x) const =0
Evaluates the negative log likelihood for input x.
pcg_extras::shuffle
void shuffle(Iter from, Iter to, RandType &&rng)
Definition: pcg_extras.hpp:553
pcg_detail::engine
Definition: pcg_random.hpp:364
hops::Model::getParameterNames
virtual std::optional< std::vector< std::string > > getParameterNames() const
Definition: Model.hpp:31
hops::IsCalculateLogAcceptanceProbabilityAvailable
Definition: IsCalculateLogAcceptanceProbabilityAvailable.hpp:8
hops
Definition: CsvReader.hpp:8
ChordStepDistributions.hpp
hops::ReversibleJumpProposal::ReversibleJumpProposal
ReversibleJumpProposal(const MarkovChainImpl &markovChainImpl, const Model &model, Eigen::VectorXi m_jumpIndices, const typename MarkovChainImpl::StateType parameterDefaultValues)
Definition: ReversibleJumpProposal.hpp:29
hops::ReversibleJumpProposal::getParameterNames
std::vector< std::string > getParameterNames() const
Definition: ReversibleJumpProposal.hpp:127
hops::ReversibleJumpProposal::getState
MarkovChainImpl::StateType getState()
Definition: ReversibleJumpProposal.hpp:116
hops::GaussianStepDistribution< double >
hops::GaussianStepDistribution::draw
RealType draw(RandomNumberGenerator &randomNumberGenerator, RealType lowerLimit, RealType upperLimit)
Definition: ChordStepDistributions.hpp:29
hops::Model
Definition: Model.hpp:12
IsCalculateLogAcceptanceProbabilityAvailable.hpp
hops::ReversibleJumpProposal::draw
void draw(RandomNumberGenerator &randomNumberGenerator)
Definition: ReversibleJumpProposal.hpp:49
hops::ReversibleJumpProposal::getAcceptanceRate
double getAcceptanceRate()
Definition: ReversibleJumpProposal.hpp:112
hops::ReversibleJumpProposal::drawInModelSpace
void drawInModelSpace(RandomNumberGenerator &randomNumberGenerator, std::vector< int > &parameterActivationStates, typename MarkovChainImpl::StateType m_defaultValues)
Definition: ReversibleJumpProposal.hpp:57
hops::ReversibleJumpProposal
Definition: ReversibleJumpProposal.hpp:27
hops::GaussianStepDistribution::computeProbabilityDensity
RealType computeProbabilityDensity(RealType x, RealType m_sigma, RealType m_lowerBound, RealType m_upperBound)
Definition: ChordStepDistributions.hpp:53
Generated by   doxygen 1.8.17