Header-only patch

.gitignore

@@ -1,6 +1,7 @@
 .Rproj.user
 .Rhistory
 .RData
+.DS_store
 src/*.o
 src/*.so
 src/*.dll

R/RcppExports.R

@@ -5,16 +5,16 @@ blockpfGaussianOpt_impl <- function(data, part, lag) {
     .Call(`_RcppSMC_blockpfGaussianOpt_impl`, data, part, lag)
 }
 
-LinReg_impl <- function(Data, lNumber) {
-    .Call(`_RcppSMC_LinReg_impl`, Data, lNumber)
+LinRegLA_adapt_impl <- function(Data, lNumber, resampTol, tempTol) {
+    .Call(`_RcppSMC_LinRegLA_adapt_impl`, Data, lNumber, resampTol, tempTol)
 }
 
 LinRegLA_impl <- function(Data, intemps, lNumber) {
     .Call(`_RcppSMC_LinRegLA_impl`, Data, intemps, lNumber)
 }
 
-LinRegLA_adapt_impl <- function(Data, lNumber, resampTol, tempTol) {
-    .Call(`_RcppSMC_LinRegLA_adapt_impl`, Data, lNumber, resampTol, tempTol)
+LinReg_impl <- function(Data, lNumber) {
+    .Call(`_RcppSMC_LinReg_impl`, Data, lNumber)
 }
 
 nonLinPMMH_impl <- function(data, lNumber, lMCMCits) {

inst/include/history.h

@@ -40,8 +40,20 @@ namespace smc {
         /// true if the particle system was resampled during the described iteration.
         unsigned int Resampled : 1;
     public:
-        ///Create a new set of history flags corresponding to the specified properties
-        historyflags(int wasResampled);
+      // ///Create a new set of history flags corresponding to the specified properties
+      // historyflags(int wasResampled);
+
+      /// This constructor produces an initialised historyflags instance.
+      ///
+      /// \param wasResampled An indicator which should be nonzero if the particle 
+      /// system was resampled during the iteration being described
+      historyflags(int wasResampled)
+      {
+	if(wasResampled)
+	  Resampled = 1;
+	else
+	  Resampled = 0;
+      }
 
         ///This function returns true if the flag set indicates that the ensemble was resampled during the described iteration.
         int WasResampled(void) {return Resampled;}

inst/include/population.h

@@ -39,7 +39,15 @@
 namespace smc {
 
     /// A stable calculation of the log sum of the weights, used in ESS calculations
-    double stableLogSumWeights(const arma::vec & logw);
+    /// This function performs a stable calculation of the log sum of the weights, which is useful for
+    /// normalising weights, calculating the effective sample size and estimating the normalising constant.
+    ///
+    /// \param logw The log weights of interest.
+    inline double stableLogSumWeights(const arma::vec & logw){
+        double dMaxWeight = arma::max(logw);
+        double sum = arma::sum(exp(logw - dMaxWeight));
+        return (dMaxWeight + log(sum));
+    }
 
     /// A template class for the particles of an SMC algorithm
     template <class Space> class population

inst/include/sampler.h

@@ -319,8 +319,7 @@ namespace smc {
             historyelement<Space> histel;
             histel.Set(N, pPopulation, nAccepted, nRepeats, historyflags(nResampled));
             History.push_back(histel);          
-        }  
-
+        }
         return;
     }
 
@@ -460,7 +459,7 @@ namespace smc {
     void sampler<Space,Params>::IterateBack(void)
     {
         if(htHistoryMode == HistoryType::NONE)
-        throw SMC_EXCEPTION(SMCX_MISSING_HISTORY, "An attempt to undo an iteration was made; unforunately, the system history has not been stored.");
+        throw SMC_EXCEPTION(SMCX_MISSING_HISTORY, "An attempt to undo an iteration was made; unfortunately, the system history has not been stored.");
 
         History.pop_back();
         historyelement<Space> recent = History.back();

inst/include/smc-exception.h

@@ -51,17 +51,42 @@ namespace smc {
     long lCode;   //!< A numerical code indicating the nature of the exception generated.
     char const * szMessage; //!< A human-readable explanation of the cause of the exception.   
 
-    exception(char const *, long, long, char const *);
+    //! Generate an SMCTC Exception class with the specified initialisation.
+
+    //! This constructor fills the four elements of the class with their specified values.
+    //! It is used to allow a single-line command to create and throw an exception.
+    //!
+    //! \param szN The name of the source file generating the exception.
+    //! \param lL The line in that file responsible for the exception.
+    //! \param lC The numerical code identifying the exception.
+    //! \param szM An textual explanation of the problem.
+    exception(char const * szN, long lL, long lC, char const * szM)
+    {
+      szFile = szN; 
+      lLine = lL; 
+      lCode = lC; 
+      szMessage = szM;
+    }
   };
 }
 
+
+
 namespace std {
   /// Produce a human-readable display of the state of an smc::exception class using the stream operator.
 
-  /// Produce a human-readable display of the state of an smc::exception class using the stream operator.
-  /// \param os The output stream to which the display should be made.
-  /// \param e  The exception which is to be displayed.
-  std::ostream & operator<< (std::ostream &, smc::exception &);
+  ///Display a human-readable version of an SMC exception.
+
+  /// \param os The stream to write to.
+  /// \param e The exception class to display.
+  /// \return os
+  inline std::ostream & operator<< (std::ostream & os, smc::exception & e)
+  {
+    os << "SMC Exception: #" << e.lCode << endl;
+    os << "Error occured in file " << e.szFile << " at line " << e.lLine << "." << endl;
+    os << "Details: " << endl << '\t' << e.szMessage << endl;
+    return os;
+  }
 }
 
 #endif

inst/include/staticModelAdapt.h

@@ -76,5 +76,98 @@ namespace smc {
         /// Returns the empirical covariance matrix based on the current weighted particle set.
         arma::mat GetEmpCov(void) const {return empCov;}
     };
+
+      /// Computes the difference between the conditional ESS given the specified temperature difference and the desired conditional ESS.
+    inline double staticModelAdapt::CESSdiff(const arma::vec & logweight, const arma::vec & loglike, double tempDiff, double desiredCESS){
+        double logsum1 = stableLogSumWeights(logweight + tempDiff*loglike);
+        double logsum2 = stableLogSumWeights(logweight + 2*tempDiff*loglike);
+
+        return exp(log(static_cast<double>(logweight.n_rows)) + 2*logsum1 - logsum2) - desiredCESS;
+    }
+
+    ///Performs the bisection method to find the temperature within (temp_curr,1) which gives the desired conditional ESS.
+    inline double staticModelAdapt::bisection(double curr, const arma::vec & logweight, const arma::vec & loglike, double desiredCESS, double epsilon){
+        double a = curr;
+        double b = 1.0;
+        double f_a = CESSdiff(logweight,loglike,a-curr,desiredCESS);
+        double f_b = CESSdiff(logweight,loglike,b-curr,desiredCESS);
+        if (f_a*f_b>0){
+            Rcpp::stop("Bisection method to choose the next temperature failed");
+        } else{
+            double m, f_m, err;
+            m = (a+b)/2.0;
+            f_m = CESSdiff(logweight,loglike,m-curr,desiredCESS);
+            err = 10.0;
+            while (err > epsilon){
+                if (f_m<0.0){
+                    b = m;
+                    f_b = f_m;
+                } else{
+                    a = m;
+                    f_a = f_m;
+                }
+                m = (a+b)/2.0;
+                f_m = CESSdiff(logweight,loglike,m-curr,desiredCESS);
+                err = std::abs(f_m);
+            }
+            return m;
+        }
+    }
+
+    /// Chooses the next temperature such that a desired conditional ESS is maintained.
+    ///
+    /// \param logweight An armadillo vector containing the logarithm of the current particle weights.
+    /// \param loglike An armadillo vector containing the log likelihood of the current particle values.
+    /// \param desiredCESS The target conditional ESS for the next temperature (generally fixed).
+    /// \param epsilon The tolerance for the bisection method (maximum difference between desired and actual conditional ESS).
+    inline void staticModelAdapt::ChooseTemp(const arma::vec & logweight, const arma::vec & loglike, double desiredCESS, double epsilon) {
+        double temp_curr = temp.back();
+        if (CESSdiff(logweight,loglike,1.0-temp_curr,desiredCESS)>=-epsilon){
+            temp.push_back(1.0);
+        } else {
+            temp.push_back(bisection(temp_curr, logweight, loglike, desiredCESS, epsilon));
+        }
+    }
+
+    /// Calculates the empirical covariance matrix based on the current weighted particle set.
+    ///
+    /// \param theta An [Nxd] armadillo matrix of doubles for the current particle values, where N is
+    /// the number of particles and d is the dimension of the parameter.
+    /// \param logweight An armadillo vector containing the logarithm of the current particle weights.
+    inline void staticModelAdapt::calcEmpCov(const arma::mat & theta, const arma::vec & logweight){
+        int N = logweight.n_rows;
+        arma::vec normWeights = exp(logweight - stableLogSumWeights(logweight));
+
+        arma::mat diff = theta - arma::ones(N,1)*arma::mean(theta,0);
+        empCov = diff.t()*diagmat(normWeights)*diff;
+    }
+
+    /// Calculates the Cholesky decomposition of the empirical covariance matrix based on the current weighted particle set.
+    ///
+    /// \param theta An [Nxd] armadillo matrix of doubles for the current particle values, where N is
+    /// the number of particles and d is the dimension of the parameter
+    /// \param logweight An armadillo vector of the logarithm of the current particle weights
+    inline void staticModelAdapt::calcCholCov(const arma::mat & theta, const arma::vec logweight){
+        calcEmpCov(theta,logweight);
+        cholCov = arma::chol(empCov);
+    }
+
+    /// Calculates the number of MCMC repeats based on the results of the most recent set of MCMC moves.
+    ///
+    /// \param acceptProb The proportion of accepted MCMC steps in the most recent iteration.
+    /// \param desiredAcceptProb The desired probability of a successful move for each particle.
+    /// \param initialN The initial number of MCMC repeats.
+    /// \param maxReps The maximum number of MCMC repeats.
+    inline int staticModelAdapt::calcMcmcRepeats(double acceptProb, double desiredAcceptProb, int initialN, int maxReps){
+        if (acceptProb + 1.0 <= 1e-9){
+            return initialN;
+        } else if (acceptProb - 1.0 >= -1e-9){
+            return 1;
+        } else if (acceptProb <= 1e-9){
+            return maxReps;
+        } else {
+            return std::min(maxReps,static_cast<int>(std::ceil(log(1.0-desiredAcceptProb)/log(1.0-acceptProb))));
+        }
+    }
 }
 #endif

src/RcppExports.cpp

@@ -19,15 +19,17 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// LinReg_impl
-Rcpp::List LinReg_impl(arma::mat Data, unsigned long lNumber);
-RcppExport SEXP _RcppSMC_LinReg_impl(SEXP DataSEXP, SEXP lNumberSEXP) {
+// LinRegLA_adapt_impl
+Rcpp::List LinRegLA_adapt_impl(arma::mat Data, unsigned long lNumber, double resampTol, double tempTol);
+RcppExport SEXP _RcppSMC_LinRegLA_adapt_impl(SEXP DataSEXP, SEXP lNumberSEXP, SEXP resampTolSEXP, SEXP tempTolSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::mat >::type Data(DataSEXP);
     Rcpp::traits::input_parameter< unsigned long >::type lNumber(lNumberSEXP);
-    rcpp_result_gen = Rcpp::wrap(LinReg_impl(Data, lNumber));
+    Rcpp::traits::input_parameter< double >::type resampTol(resampTolSEXP);
+    Rcpp::traits::input_parameter< double >::type tempTol(tempTolSEXP);
+    rcpp_result_gen = Rcpp::wrap(LinRegLA_adapt_impl(Data, lNumber, resampTol, tempTol));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -44,17 +46,15 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// LinRegLA_adapt_impl
-Rcpp::List LinRegLA_adapt_impl(arma::mat Data, unsigned long lNumber, double resampTol, double tempTol);
-RcppExport SEXP _RcppSMC_LinRegLA_adapt_impl(SEXP DataSEXP, SEXP lNumberSEXP, SEXP resampTolSEXP, SEXP tempTolSEXP) {
+// LinReg_impl
+Rcpp::List LinReg_impl(arma::mat Data, unsigned long lNumber);
+RcppExport SEXP _RcppSMC_LinReg_impl(SEXP DataSEXP, SEXP lNumberSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< arma::mat >::type Data(DataSEXP);
     Rcpp::traits::input_parameter< unsigned long >::type lNumber(lNumberSEXP);
-    Rcpp::traits::input_parameter< double >::type resampTol(resampTolSEXP);
-    Rcpp::traits::input_parameter< double >::type tempTol(tempTolSEXP);
-    rcpp_result_gen = Rcpp::wrap(LinRegLA_adapt_impl(Data, lNumber, resampTol, tempTol));
+    rcpp_result_gen = Rcpp::wrap(LinReg_impl(Data, lNumber));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -100,9 +100,9 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_RcppSMC_blockpfGaussianOpt_impl", (DL_FUNC) &_RcppSMC_blockpfGaussianOpt_impl, 3},
-    {"_RcppSMC_LinReg_impl", (DL_FUNC) &_RcppSMC_LinReg_impl, 2},
-    {"_RcppSMC_LinRegLA_impl", (DL_FUNC) &_RcppSMC_LinRegLA_impl, 3},
     {"_RcppSMC_LinRegLA_adapt_impl", (DL_FUNC) &_RcppSMC_LinRegLA_adapt_impl, 4},
+    {"_RcppSMC_LinRegLA_impl", (DL_FUNC) &_RcppSMC_LinRegLA_impl, 3},
+    {"_RcppSMC_LinReg_impl", (DL_FUNC) &_RcppSMC_LinReg_impl, 2},
     {"_RcppSMC_nonLinPMMH_impl", (DL_FUNC) &_RcppSMC_nonLinPMMH_impl, 3},
     {"_RcppSMC_pfLineartBS_impl", (DL_FUNC) &_RcppSMC_pfLineartBS_impl, 4},
     {"_RcppSMC_pfNonlinBS_impl", (DL_FUNC) &_RcppSMC_pfNonlinBS_impl, 2},