tristanz / ScytheMCMC / changeset / 9844e0828091

Up to file-list ScytheBugs.txt:

@@ -12,8 +12,23 @@ Bugs:
   Compiling the help using the latest doxygen seems to correct the problem,
   although macros don't get documented really.
 - Docs for scythe::min refer to max.
 - What is linesearch? Docs are not very clear.
 - Reading in from file is fragile.  outsheet in stata on OSX does not
   use the proper newline characters.  newlines across dos/unix/osx differ.
 - col_interchange doesn't sound like it is what it is in help file.  It swaps
 	rows rather than reordering.
 - qnorm1 has return type "double double sd double" in docs.
 Missing / Potential Features:
 - rtnorm is parameterized using the variance whereas rnorm is parameterized
   with sd.  Possibility of errors.
 - Fast versions of functions proportional to log density.  In most MCMC
   applications sum(log(scythe::dnorm)), for instance, takes up a significant
   amount of time.  However we usually just need the log of something
@@ -64,4 +79,3 @@ Missing / Potential Features:
   etc. Then add a -DHAS_MKL flag. See:
 	http://www.intel.com/software/products/mkl/docs/WebHelp/mkl.htm

Up to file-list src/mcmc.h:

 #include <iostream>
 #include <fstream>
 #include <scythestat/rng/mersenne.h>
 #include <scythestat/rng/lecuyer.h>
 #include <scythestat/distributions.h>
 #include <scythestat/matrix.h>
 #include <scythestat/rng.h>
 #include <scythestat/smath.h>
 #include <scythestat/stat.h>
 #include <scythestat/ide.h>
 #include <scythestat/defs.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <sys/stat.h>
 #include "Simple/SimpleOpt.h"
 #include "Simple/SimpleIni.h"
 // #include <mkl_vml_functions.h>
 // #include <mkl_vml_defines.h>
 #include "scythe_extensions.h"
 /// Scythe MCMC version.
 #define SCYTHE_MCMC_VERSION "0.1"
@@ -118,9 +125,9 @@ typedef pair<int, int> ipair;
 double dInf = std::numeric_limits<double>::infinity();
 double iInf = std::numeric_limits<int>::infinity();
 /// Global Mersenne-Twister random number generator.
 mersenne myrng;
 /// Global L'Ecuyer random number generator.
 lecuyer myrng;
 //mersenne myrng;
 /*! \brief Basic MCMC options.
+ *
@@ -137,7 +144,7 @@ struct MCMCOptions {
   /// Chains.
   int chains;
   /// Seed for random number generator.
-  int random_seed;
+  unsigned long random_seed [6];
   /// Config file (.ini style).
   string config_file;
   /// Out file to save parameter chains.
@@ -189,6 +196,18 @@ inline string to_csv (const T& t) {
   return s;
+}
 /*! \brief Convert any streaming type to string
+ *
  */
 template <class T>
 inline string to_string (const T& t) {
   std::stringstream ss;
   ss << t;
   string s = ss.str();
   return s;
+}
 /*! \brief Shows command line usage.
    Diplays:
@@ -202,11 +221,11 @@ inline string to_csv (const T& t) {
        --chains=arg                 number of chains
        --sample-size=arg            retained sample size
        --burnin=arg                 burn in period
        --thin=arg                   thinning iterval
        --random-seed=arg            random number seed
        --thin=arg                   thinning iterval (1 = no thinning)
        --random-seed=arg            random number seed (0 uses current timestamp)
      Example: ./base_model --config-file=base_model.cfg --out-file=out.txt --chains=1 --sample-size=1000
               --burnin=0 --thin=1 --random-seed=12345
      Example: ./command --config-file=command.ini --out-file=out.txt --chains=1 --sample-size=1000
               --burnin=0 --thin=1 --random-seed=0
 </pre>
 */
 void ShowUsage() {
@@ -217,11 +236,11 @@ void ShowUsage() {
        << "    --chains=arg \t\t number of chains" << endl
        << "    --sample-size=arg \t\t retained sample size" << endl
        << "    --burnin=arg \t\t burn in period" << endl
        << "    --thin=arg \t\t\t thinning iterval (0 = no thinning)" << endl
        << "    --random-seed=arg \t\t random number seed" << endl << endl;
        << "    --thin=arg \t\t\t thinning iterval (1 = no thinning)" << endl
        << "    --random-seed=arg \t\t random number seed (0 uses current timestamp, 1 attempts to load from lecuyer.seed file)" << endl << endl;
   cout << "  Example: ./base_model --config-file=base_model.cfg --out-file=out.txt --chains=1 --sample-size=1000" << endl
        << "\t   --burnin=0 --thin=1 --random-seed=12345" << endl << endl;
   cout << "  Example: ./command --config-file=command.ini --out-file=out.txt --chains=1 --sample-size=1000" << endl
        << "\t   --burnin=0 --thin=1 --random-seed=0" << endl << endl;
+}
 /*! \brief Parses command line options.
@@ -285,7 +304,36 @@ void DefaultStartUp(int argc, char* argv
         break;
       case OPT_RANDOM_SEED:
         setflags++;
-        mcmc_options.random_seed = atoi(args.OptionArg());
+        unsigned long seed = atoi(args.OptionArg());
 				cout << "TEST" << seed<<endl;
 				if(seed == 0) {
 					seed = time(NULL);
 					for(int i = 0; i < 6; ++i) {
 						mcmc_options.random_seed[i] = seed;
+					}
+				}
 				else if (seed == 1) {
 					// try to load seed from file
 					cout << "Try loading file..." << endl;
 					if (FileExists("lecuyer.seed")) {
 						Matrix<unsigned long, Col> seedmat("lecuyer.seed");
 						for(int i = 0; i < 6; ++i) {
 							mcmc_options.random_seed[i] = seedmat(i);
+						}
+					}
 					else {
 						cout << "No seed exists..." << endl;
 						seed = time(NULL);
 						for(int i = 0; i < 6; ++i) {
 							mcmc_options.random_seed[i] = seed;
+						}
+					}
+				}
 				else {
 					for(int i = 0; i < 6; ++i) {
 						mcmc_options.random_seed[i] = seed;
+					}
+				}
         break;
+      }
     } else {
@@ -309,8 +357,8 @@ void DefaultStartUp(int argc, char* argv
     cout << "Burn in cannot be negative.\n";
     exit(1);
+  }
   if (mcmc_options.thin < 0) {
     cout << "Thinning interval cannot be negative.\n";
   if (mcmc_options.thin < 1) {
     cout << "Thinning interval must be positive.\n";
     exit(1);
+  }
   if (!FileExists(mcmc_options.config_file)) {
@@ -319,8 +367,9 @@ void DefaultStartUp(int argc, char* argv
+  }
   // Seed random number generator.
   myrng.initialize(mcmc_options.random_seed);
   //myrng.initialize(mcmc_options.random_seed);
 	myrng.SetSeed(mcmc_options.random_seed);
   // Print mcmc options.
   cout << "Arguments:" << endl;
   cout << "  Config file: " << mcmc_options.config_file << endl
@@ -328,7 +377,13 @@ void DefaultStartUp(int argc, char* argv
        << "  Desired sample size: " << mcmc_options.sample_size << endl
        << "  Burn in period: " << mcmc_options.burnin << endl
        << "  Thinning interval: " << mcmc_options.thin << endl
-       << "  Random number seed: " << mcmc_options.random_seed << endl << endl;
+       << "  Random number seed: "
  						<< mcmc_options.random_seed[0] << " "
  						<< mcmc_options.random_seed[1] << " "
 						<< mcmc_options.random_seed[2] << " "
 						<< mcmc_options.random_seed[3] << " "
 						<< mcmc_options.random_seed[4] << " "
 						<< mcmc_options.random_seed[5] << endl << endl;
+}
@@ -455,7 +510,7 @@ public:
   /// Save a new value of the parameter.
   /// \param new_value New value to save.
-  virtual void Save(double new_value) {};
+  virtual void Save(ParameterStorageType new_value) {};
   /// Parameter is tracked / saved.
   /// \return True if parameter is tracked.
@@ -505,7 +560,7 @@ public:
+   *
    * Executes the main step function, such as taking a MH step or deterministic step.
    */
-  virtual void DoStep() {}
+   virtual void DoStep() {}
   /*! \brief Set starting value.
+   *
    *  Calls the parameter's Parameter::Starting function and then saves the result using Parameter::Save.
@@ -562,7 +617,7 @@ public:
    *  Draws from the parameter's Parameter::RandomPosterior function and then saves the result using Parameter::Save.
    */
   void DoStep() {
-    parameter_.Save(parameter_.RandomPosterior());
+		parameter_.Save(parameter_.RandomPosterior());
+  }
   void Start() {
@@ -614,12 +669,9 @@ public:
   /*! \brief Set function step starting value.
+   *
    * \note Because the starting value uses the normal function calculation, it is important to add
    * parameters in an order that makes the function feasible.  StartingValue functions are called in the order
    * they are added.  Any implemention of a StartingValue method for a FunctionStep will not be used.
    */
   void Start() {
-    DoStep();
+		parameter_.StartingValue();
+  }
   string ParameterLabel() {
@@ -815,7 +867,7 @@ public:
+   *
    *  Calculates the value of Parameter::Function and then saves the result using Parameter::Save.
    */
-  void DoStep() {
+   void DoStep() {
     // Draw a new parameter
     double new_value = proposal_.Draw(parameter_.Value());
     double old_value = parameter_.Value();
@@ -878,14 +930,16 @@ public:
    *  \param lower Lower bound on the support of the distribution (default = -Infinity)
    *  \param upper Upper bound on the support of the distribution (default = Infinity)
    */
   SliceStep(ParameterType parameter, double w = 1.0, double lower = -dInf, double upper = dInf) :
       parameter_(parameter), w_(w), lower_(lower), upper_(upper) {}
   SliceStep(ParameterType parameter, double w, double lower, double upper) :
       parameter_(parameter), w_(w), lower_(lower), upper_(upper) {
 		step_count_ = 0;
+	}
   /*! \brief Take a slice sampling step for the parameter.
+   *
    */
   void DoStep() {
     // Find the log density at the initial point.
    void DoStep() {
 	  // Find the log density at the initial point.
     double x0 = parameter_.Value();
     double gx0 = parameter_.LogDensity(x0);
@@ -945,7 +999,10 @@ public:
         L = x1;
+      }
+    }
 		// Adapt
 		//step_count_++;
 		//w_ = (step_count_/(step_count_+1)) * w_ + (1/(step_count_+1))*(R - L);
     // Save the point sampled
     parameter_.Save(x1);
+  }
@@ -972,6 +1029,7 @@ private:
   double w_;
   double lower_;
   double upper_;
 	int step_count_;
 };
 /*! \brief MCMC sampler.
@@ -1015,23 +1073,25 @@ public:
    * \param number_of_iterations Number of iterations to run sampler.
    * \param progress Display a progress bar.
    */
   void Iterate(int number_of_iterations, bool progress = false) {
     if (progress) {
       boost::progress_display show_progress(number_of_iterations);
       for (int iter = 0; iter < number_of_iterations; ++iter) {
         for (int i = 0; i < steps_.size(); ++i) {
           steps_[i].DoStep();
+        }
         ++show_progress;
+      }
     } else {
       for (int iter = 0; iter < number_of_iterations; ++iter) {
         for (int i = 0; i < steps_.size(); ++i) {
           steps_[i].DoStep();
+        }
+      }
+    }
+  }
 	 void Iterate(int number_of_iterations, bool progress = false) {
 	   if(progress) {
 	     boost::progress_display show_progress(number_of_iterations);
 	     for(int iter = 0; iter < number_of_iterations; ++iter) {
 	       for(int i = 0; i < steps_.size(); ++i) {
 	         steps_[i].DoStep();
+	       }
 	       ++show_progress;
+	     }
+	   }
 	   else {
 	     for(int iter = 0; iter < number_of_iterations; ++iter) {
 	       for(int i = 0; i < steps_.size(); ++i) {
 	         steps_[i].DoStep();
+	       }
+	     }
+	   }
+ 	}
   /*! \brief Run sampler.
+   *
    * Run the MCMC sampling procedure, including burning in, sampling, thinning, displaying progress,
@@ -1075,7 +1135,7 @@ public:
     cout << endl << "Sampling..." << endl;
     boost::progress_display show_progress(sample_size_);
     for (int i = 0; i < sample_size_; ++i) {
-      Iterate(thin_ + 1);
+      Iterate(thin_);
       // Save data
       for (int k = 0; k < tracks_.size(); ++k) {
         if (k == tracks_.size() - 1) {
@@ -1089,7 +1149,18 @@ public:
+    }
     outfile.close();
 		// Save lecuyer seed for future use
 		Matrix<unsigned long, Col> seedmat(6, 1, false);
 		unsigned long seed[6];
 		myrng.GetState(seed);
 		cout << "Saving seed: " << endl << seed[0] << " " << seed[1] << " " << seed[2]
 			<< " " << seed[3] << " " << seed[4] << " " << seed[5] << endl;
 		for(int i = 0; i < 6; ++i) {
 			seedmat(i) = seed[i];
+		}
 		seedmat.save("lecuyer.seed", 'o');
     cout << "Total elapsed time: " << timer.elapsed() << " seconds" << endl;
+  }
   /*! \brief Number of steps in one sampler iteration.
@@ -1142,3 +1213,4 @@ double fast_log_dnorm(const double x, co
+}

commit 13:	9844e0828091
parent 12:	fd1ece924824
branch:	default
tags:	tip

tristanz / ScytheMCMC (http://bytebucket.org/tristanz/scythemcmc/wiki/html/index.html)