Improved const correctness. Simplified MCMC iteration method.

src/mcmcs.cpp

@@ -17,7 +17,7 @@ namespace fs = std::filesystem;
 
 
 // ================================================================================
-// Interface for different MCMC methods
+// MCMC Interface
 //
 // ================================================================================
 
@@ -33,10 +33,10 @@ MCMC::MCMC(
     n_burn_iters(100),
     n_sample_iters(900),
     n_total_iters(n_burn_iters + n_sample_iters),
-    acceptance_rate(0.0),
+    acceptance_rate(-1.0),
     acceptance_trace(n_total_iters),
     logposterior_trace(n_total_iters),
-    particle_trace(n_total_iters)
+    particle_trace(n_total_iters, params.K)  // TODO: Here is where space gets allocated. Double check.
 {};
 
 
@@ -87,6 +87,10 @@ Particle MCMC::get_map_particle() const
 }
 
 
+MCMC::~MCMC()
+{}
+
+
 // ================================================================================
 // Concrete MCMC methods
 //
@@ -98,64 +102,75 @@ Particle MCMC::get_map_particle() const
 
 
 MetropolisHastings::MetropolisHastings(
-        const Parameters& params, 
-        const Model& model, 
-        ProposalEngine& proposal_engine)
+    const Parameters& params, 
+    const Model& model, 
+    ProposalEngine& proposal_engine)
     : MCMC(params, model, proposal_engine)
 {};
 
 
+void MetropolisHastings::run_burn()
+{
+    // Random initialisation
+    ix = 0;
+    particle_trace[ix] = proposal_engine.create_particle();
+    logposterior_trace[ix] = model.calc_logposterior(particle_trace[ix]);
+    acceptance_rate = 1.0;
+    acceptance_trace[ix] = 1.0;
+    ++ix;
+
+    // Run burn-in iterations
+    run_iterations(n_burn_iters - 1);
+}
+
+
+void MetropolisHastings::run_sampling()
+{
+    run_iterations(n_sample_iters);
+}
+
+
 void MetropolisHastings::run_iterations(int n)
 {
     // Check that ix > 0
     if (ix == 0) {
-        throw std::invalid_argument("Must have initialised beforee running iterations.");
+        throw std::invalid_argument("Initialise the particles before iterating.");
     }
+
+    // Current state
+    Particle particle = particle_trace[ix - 1];
+    double logposterior = logposterior_trace[ix - 1];
 
     // Iterate
     int N = ix + n;
     for (; ix < N; ++ix) {
 
         // Propose
-        Particle proposed_particle = proposal_engine.propose_particle(particle_trace[ix-1]);
+        Particle proposed_particle = proposal_engine.propose_particle(particle);
         double proposed_logposterior = model.calc_logposterior(proposed_particle);
 
-        // Compute Hastings ratio
-        double A = exp(proposed_logposterior - logposterior_trace[ix-1]);
-        double u = U(rng.engine); 
+        // Compute acceptance probability
+        // TODO: Important to verify in the same base, else biased
+        double A = std::exp(proposed_logposterior - logposterior); 
+        double u = U(rng.engine);
+
+        // Accept
         if (u < A) {
-            logposterior_trace[ix] = proposed_logposterior;
-            particle_trace[ix] = proposed_particle;
-        } else {
-            logposterior_trace[ix] = logposterior_trace[ix-1];
-            particle_trace[ix] = particle_trace[ix-1];
+            particle = proposed_particle; // TODO: No longer need proposed particle, should be a move operation
+            logposterior = proposed_logposterior;
         }
 
+        // Store
+        particle_trace[ix] = particle;
+        logposterior_trace[ix] = logposterior;
+
         // Track expected acceptance rate
         acceptance_rate += (A < 1.0 ? A : 1.0);
         acceptance_trace[ix] = acceptance_rate / ix;
     }
 }
 
 
-void MetropolisHastings::run_burn()
-{
-    // Initialise
-    particle_trace[ix] = proposal_engine.create_particle();
-    logposterior_trace[ix] = model.calc_logposterior(particle_trace[ix]);
-    acceptance_trace[ix] = 0.0;
-    ++ix;
-
-    run_iterations(n_burn_iters - 1);
-}
-
-
-void MetropolisHastings::run_sampling()
-{
-    run_iterations(n_sample_iters);
-}
-
-
 void MetropolisHastings::run()
 {
     run_burn();

src/mcmcs.hpp

@@ -30,9 +30,9 @@ class MCMC
 
     // Parameters
     int ix;                                     // Iteration index
-    int n_burn_iters;                           // Number of burn-in iterations
-    int n_sample_iters;                         // Number of sampling iterations
-    int n_total_iters;
+    const int n_burn_iters;                           // Number of burn-in iterations
+    const int n_sample_iters;                         // Number of sampling iterations
+    const int n_total_iters;
     double acceptance_rate;                     // Accept rate until `ix`
 
     // Storage
@@ -54,12 +54,13 @@ class MCMC
     // Abstract
     virtual void run() = 0;
 
-    // Concrete
     void write_output(
         const string& output_dir, 
         const ParticleWriter& particle_writer) const;
 
     Particle get_map_particle() const;
+
+    virtual ~MCMC();
 };
 
 
@@ -89,3 +90,63 @@ class MetropolisHastings : public MCMC
     void run() override;
 };
 
+
+
+// --------------------------------------------------------------------------------
+// Parallel Tempering
+// --------------------------------------------------------------------------------
+
+
+// class ParallelTempering : public MCMC
+// {
+// private:
+//     /**
+//     * Define a temperature level
+//     */
+//     struct TemperatureLevel
+//     {
+//         double beta;
+//         Particle* particle_ptr;
+//         double loglikelihood;  // Needed for swaps; Thermodynamic Integration (TI)
+//         double beta_logposterior;   // Needed for within-level MH updates
+
+//         TemperatureLevel();    // TODO: Hmm.. what goes in constructor?
+//     };
+
+//     std::vector<ParallelTempering::TemperatureLevel>  static create_temp_levels(
+//         std::vector<Particle>& particles, 
+//         double lambda=0.5
+//     );
+
+//     void run_iterations(int n);
+//     void run_burn();
+//     void run_sampling();
+
+// public:
+//     // DATA MEMBERS
+//     const int n_temps;                    // Number of temperature levels
+//     const int swap_freq;                  // Number of iterations per swap attempt
+//     std::vector<Particle> particles;      // Particle for each temperature
+//     std::vector<TemperatureLevel> temps;  // Temperature information
+//     MatrixXd loglikelihoods;              // Loglikelihoods; for TI; TODO: should be array?
+
+//     //  Recording swaps
+//     double n_swap_attempts;               // No. swap attempts; same for all pairs of levels
+//     ArrayXd n_swaps;                      // No. swaps for each pair, up to current `ix`
+//     MatrixXd swap_rates;                  // Rate of swapping for each pair of temp. levels
+
+//     ParallelTempering(
+//         const Parameters& params, 
+//         const Model& model, 
+//         ProposalEngine& proposal_engine,
+//         int n_temps
+//     );
+
+//     void run() override;
+
+//     // Concrete
+//     void write_output(
+//         const std::string& output_dir, 
+//         const ParticleWriter& particle_writer) const override;
+// };
+