--- trunk/gadget/optinfo.h	2015/07/23 19:00:38	11
+++ trunk/gadget/optinfo.h	2017/07/05 12:08:06	28
@@ -5,9 +5,10 @@
 #include "doublematrix.h"
 #include "doublevector.h"
 #include "intvector.h"
+#include "intmatrix.h"
 #include "seq_optimize_template.h"
 
-enum OptType { OPTHOOKE = 1, OPTSIMANN, OPTBFGS };
+enum OptType { OPTHOOKE = 1, OPTSIMANN, OPTBFGS, OPTPSO };
 
 /**
  * \class OptInfo
@@ -41,9 +42,17 @@
    * \brief This is the function used to call the optimisation algorithms
    */
   virtual void OptimiseLikelihood() {};
-  //FIXME doc
+#ifdef _OPENMP
+  /**
+   * \brief This is the function used to call the optimisation algorithms parallelized with OpenMP of the reproducible version
+   */
   virtual void OptimiseLikelihoodOMP() {};
-  //FIXME doc
+  virtual void OptimiseLikelihoodREP() {};
+#endif
+  /**
+   * \brief This function set the seeds used in SA
+   * \param val array of unsigned int with the seeds
+   */
   void setSeed(unsigned* val) {seed = val[0]; seedM =  val[1]; seedP = val[2];};
   /**
    * \brief This will return the type of optimisation class
@@ -67,9 +76,17 @@
    * \brief This denotes what type of optimisation class has been created
    */
   OptType type;
-  //FIXME doc
+  /**
+   * \brief This is the seed used for the calculation of the new value of the parameters
+   */
   unsigned seed;
+  /**
+   * \brief This is the seed used for the acceptance of the metropolis
+   */
   unsigned seedM;
+  /**
+   * \brief This is the seed used to change the order of the parameters
+   */
   unsigned seedP;
 };
 
@@ -108,9 +125,13 @@
    * \brief This is the function that will calculate the likelihood score using the Hooke & Jeeves optimiser
    */
   virtual void OptimiseLikelihood();
-#ifdef GADGET_OPENMP
-  //TODO doc
+#ifdef _OPENMP
+  /**
+   * \brief This is the function that will calculate the likelihood score using the Hooke & Jeeves optimiser parallelized with the reproducible version implemented OpenMP
+   */
+
   virtual void OptimiseLikelihoodOMP();
+  virtual void OptimiseLikelihoodREP();
 #endif
 private:
   /**
@@ -122,9 +143,24 @@
    * \return the best function value found from the search
    */
   double bestNearby(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param);
-  //TODO doc
-  double bestNearbyOMP(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param);
-  double bestNearbyOMP2(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param);
+  /**
+   * \brief This function implemented the reproducible version with OpenMP will calculate the best point that can be found close to the current point
+   * \param delta is the DoubleVector of the steps to take when looking for the best point
+   * \param point is the DoubleVector that will contain the parameters corresponding to the best function value found from the search
+   * \param prevbest is the current best point value
+   * \param param is the IntVector containing the order that the parameters should be searched in
+   * \return the best function value found from the search
+   */
+  double bestNearbyRepro(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param);
+  /**
+     * \brief This function implemented the speculative version with OpenMP will calculate the best point that can be found close to the current point
+     * \param delta is the DoubleVector of the steps to take when looking for the best point
+     * \param point is the DoubleVector that will contain the parameters corresponding to the best function value found from the search
+     * \param prevbest is the current best point value
+     * \param param is the IntVector containing the order that the parameters should be searched in
+     * \return the best function value found from the search
+     */
+  double bestNearbySpec(DoubleVector& delta, DoubleVector& point, double prevbest, IntVector& param);
   /**
    * \brief This is the maximum number of iterations for the Hooke & Jeeves optimisation
    */
@@ -181,17 +217,28 @@
    * \brief This is the function that will calculate the likelihood score using the Simulated Annealing optimiser
    */
   virtual void OptimiseLikelihood();
-  //FIXME doc
-#ifdef GADGET_OPENMP
+#ifdef _OPENMP
+//#ifdef SPECULATIVE
+  /**
+   * \brief This is the function that will calculate the likelihood score using the Simulated Annealing optimiser parallelized with the reproducible version implemented OpenMP
+   */
   virtual void OptimiseLikelihoodOMP();
+  virtual void OptimiseLikelihoodREP();
+//#endif
 #endif
-//FIXME doc
+  /**
+   * \brief This function calculate a new valor for the parameter l
+   * \param nvars the number of variables to be optimised
+   * \param l the parameter to change
+   * \param param IntVector with the order of the parameters
+   * \param trialx DoubleVector that storage the values of the parameters to evaluate during this iteration
+   * \param x DoubleVector that storage the best values of the parameters
+   * \param lowerb DoubleVector with the lower bounds of the variables to be optimised
+   * \param upperb DoubleVector with the upper bounds of the variables to be optimised
+   * \param vm DoubleVector with the value for the maximum step length for each parameter
+   */
   virtual void newValue(int nvars, int l, IntVector& param, DoubleVector& trialx,
   		DoubleVector& x, DoubleVector& lowerb, DoubleVector& upperb, DoubleVector& vm);
-//  //FIXME doc
-//  virtual void buildNewParams_f(Siman& seed, DoubleVector& params);
-//  virtual void adjustVm(Siman& seed);
-//  virtual void temperature(Siman& seed);
 private:
   /**
    * \brief This is the temperature reduction factor
@@ -277,9 +324,14 @@
    * \brief This is the function that will calculate the likelihood score using the BFGS optimiser
    */
   virtual void OptimiseLikelihood();
-  //FIXME doc
-#ifdef GADGET_OPENMP
+#ifdef _OPENMP
+//#ifdef SPECULATIVE
+  /**
+  * \brief This function call the sequential function. BFGS isn't implemented with OpenMP
+  */
   virtual void OptimiseLikelihoodOMP();
+  virtual void OptimiseLikelihoodREP();
+//#endif
 #endif
 private:
   /**
@@ -329,4 +381,189 @@
   double gradeps;
 };
 
+/**
+ * \class OptInfoPso
+ * \brief This is the class used for the PSO optimisation
+ *
+ * PSO or Particle Swarm Optimization 
+ *
+ * The PSO algorithm used in Gadget is derived from that presented by Kyriakos Kentzoglanakis.
+ */
+class OptInfoPso : public OptInfo  {
+public:
+  /**
+   * \brief This is the default OptInfoBFGS constructor
+   */
+  OptInfoPso();
+  /**
+   * \brief This is the default OptInfoBFGS destructor
+   */
+  ~OptInfoPso() {};
+  /**
+   * \brief This is the function used to read in the BFGS parameters
+   * \param infile is the CommentStream to read the optimisation parameters from
+   * \param text is a text string used to compare parameter names
+   */
+  virtual void read(CommentStream& infile, char* text);
+  /**
+   * \brief This function will print information from the optimisation algorithm
+   * \param outfile is the ofstream that the optimisation information gets sent to
+   * \param prec is the precision to use in the output file
+   */
+  virtual void Print(ofstream& outfile, int prec);
+  /**
+   * \brief This is the function that will calculate the likelihood score using the PSO optimiser
+   */
+  virtual void OptimiseLikelihood();
+#ifdef _OPENMP
+//#ifdef SPECULATIVE
+  /**
+  * \brief This function call the sequential function. PSO isn't implemented with OpenMP
+  */
+  virtual void OptimiseLikelihoodOMP();
+  virtual void OptimiseLikelihoodREP();
+//#endif
+#endif
+private:
+
+ /**
+  * \brief CONSTANTS: max swarm size
+  */
+ #define PSO_MAX_SIZE 100
+ 
+ /**
+  * \brief CONSTANTS: default value of w (see clerc02) 
+  */
+ #define PSO_INERTIA 0.7298 
+ 
+ 
+ 
+ /**
+  * \brief NEIGHBORHOOD SCHEMES: global best topology 
+  */
+#define PSO_NHOOD_GLOBAL 0
+ 
+ /**
+  * \brief NEIGHBORHOOD SCHEMES: ring topology 
+  */
+#define PSO_NHOOD_RING 1
+ 
+ /**
+  * \brief NEIGHBORHOOD SCHEMES: Random neighborhood topology. see http://clerc.maurice.free.fr/pso/random_topology.pdf 
+  */
+#define PSO_NHOOD_RANDOM 2
+ 
+ 
+ 
+ /**
+  * \brief INERTIA WEIGHT UPDATE FUNCTIONS 
+  */
+#define PSO_W_CONST 0
+#define PSO_W_LIN_DEC 1
+
+
+
+ /**
+  * \brief PSO SOLUTION -- Initialized by the user
+  */
+ typedef struct {
+ 
+     double error;
+     double *gbest; // should contain DIM elements!!
+ 
+ } pso_result_t;
+ 
+ /**
+  * \brief optimization goal (error threshold) 
+  */
+ double goal; 
+ 
+ /**
+  * \brief swarm size (number of particles)
+  */
+ int size; 
+ 
+ /**
+  * \brief maximum number of iterations
+  */
+ int psoiter; 
+ 
+ /**
+  * \brief cognitive coefficient
+  */
+ double c1; 
+ 
+ /**
+  * \brief social coefficient
+  */
+ double c2; 
+ 
+ /**
+  * \brief max inertia weight value
+  */
+ double w_max;
+ 
+ /**
+  * \brief min inertia weight value
+  */
+ double w_min; 
+ 
+ /**
+  * \brief whether to keep particle position within defined bounds (TRUE) or apply periodic boundary conditions (FALSE)
+  */
+ int clamp_pos; 
+ 
+ /**
+  * \brief neighborhood strategy (see PSO_NHOOD_*)
+  */
+ int nhood_strategy;
+ 
+ /**
+  * \brief neighborhood size
+  */
+ int nhood_size;  
+ 
+ /**
+  * \brief inertia weight strategy (see PSO_W_*)
+  */
+ int w_strategy;
+
+ /**
+    * \brief This is the flag to denote whether the parameters should be scaled or not (default 0, not scale)
+    */
+ int scale;
+  
+ 
+// /**
+//  * \brief seed for the generator
+//  */
+// long seed;
+ 
+ 
+ 
+
+ /**
+  * \brief return the swarm size based on dimensionality
+  */
+ int pso_calc_swarm_size(int dim);
+ 
+ double calc_inertia_const(int step);
+ double calc_inertia_lin_dec(int step);
+ void inform_global(IntMatrix& comm, DoubleMatrix& pos_nb,
+                   DoubleMatrix& pos_b, DoubleVector& fit_b,
+                   DoubleVector& gbest, int improved);
+ void inform_ring(IntMatrix& comm, DoubleMatrix& pos_nb,
+                 DoubleMatrix& pos_b, DoubleVector& fit_b,
+                 DoubleVector& gbest, int improved);
+ void inform_random(IntMatrix& comm, DoubleMatrix& pos_nb,
+                   DoubleMatrix& pos_b, DoubleVector& fit_b,
+                   DoubleVector& gbest,  int improved);
+ void inform(IntMatrix& comm, DoubleMatrix& pos_nb, DoubleMatrix& pos_b, DoubleVector& fit_b, int improved);
+ void init_comm_ring(IntMatrix& comm);
+ void init_comm_random(IntMatrix& comm) ;
+
+ typedef void (OptInfoPso::*Inform_fun)(IntMatrix&, DoubleMatrix&, DoubleMatrix&, DoubleVector&, DoubleVector&,  int); // neighborhood update function
+ typedef double (OptInfoPso::*Calc_inertia_fun)(int); // inertia weight update function
+
+};
 #endif