--- trunk/gadget/optinfo.h	2016/05/25 16:36:33	19
+++ trunk/gadget/optinfo.h	2017/04/07 09:20:55	20
@@ -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,10 +42,13 @@
    * \brief This is the function used to call the optimisation algorithms
    */
   virtual void OptimiseLikelihood() {};
+#ifdef _OPENMP
   /**
    * \brief This is the function used to call the optimisation algorithms parallelized with OpenMP of the reproducible version
    */
   virtual void OptimiseLikelihoodOMP() {};
+  virtual void OptimiseLikelihoodREP() {};
+#endif
   /**
    * \brief This function set the seeds used in SA
    * \param val array of unsigned int with the seeds
@@ -121,11 +125,13 @@
    * \brief This is the function that will calculate the likelihood score using the Hooke & Jeeves optimiser
    */
   virtual void OptimiseLikelihood();
-#ifdef SPECULATIVE
+#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:
   /**
@@ -211,11 +217,14 @@
    * \brief This is the function that will calculate the likelihood score using the Simulated Annealing optimiser
    */
   virtual void OptimiseLikelihood();
-#ifdef SPECULATIVE
+#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
   /**
    * \brief This function calculate a new valor for the parameter l
@@ -315,11 +324,14 @@
    * \brief This is the function that will calculate the likelihood score using the BFGS optimiser
    */
   virtual void OptimiseLikelihood();
-#ifdef SPECULATIVE
+#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:
   /**
@@ -369,4 +381,184 @@
   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 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