Steps towards dealing with parallel processing.

admin/yacc.sh

@@ -1,5 +1,15 @@
 #!/usr/bin/env bash
-# Copyright (C) 2010 - 2014 by Pedro Mendes, Virginia Tech Intellectual 
+# Copyright (C) 2019 - 2021 by Pedro Mendes, Rector and Visitors of the 
+# University of Virginia, University of Heidelberg, and University 
+# of Connecticut School of Medicine. 
+# All rights reserved. 
+
+# Copyright (C) 2017 - 2018 by Pedro Mendes, Virginia Tech Intellectual 
+# Properties, Inc., University of Heidelberg, and University of 
+# of Connecticut School of Medicine. 
+# All rights reserved. 
+
+# Copyright (C) 2010 - 2016 by Pedro Mendes, Virginia Tech Intellectual 
 # Properties, Inc., University of Heidelberg, and The University 
 # of Manchester. 
 # All rights reserved. 
@@ -42,6 +52,7 @@ sed -e 's/'$FILE_PREFIX'parse/yyparse/g' \
     -e '/#define yylex/d' \
     -e '/int yyparse (.*);/d' \
     -e 's/'$FILE_PREFIX'.tab.cpp/'$TARGET_FILE_C'/g' \
+    -e 's/'$FILE_PREFIX'.tab.hpp/'$TARGET_FILE_H'/g' \
     -e 's/int yydebug;/int yydebug = YYDEBUG;/' \
     -e '/getenv()/d' \
     ${TARGET_FILE_C} > $$.tmp && \

copasi/OpenMP/CContext.h

@@ -1,4 +1,4 @@
-// Copyright (C) 2020 by Pedro Mendes, Rector and Visitors of the
+// Copyright (C) 2020 - 2021 by Pedro Mendes, Rector and Visitors of the
 // University of Virginia, University of Heidelberg, and University
 // of Connecticut School of Medicine.
 // All rights reserved.
@@ -34,7 +34,7 @@ typedef int MPI_Win;
 template < class Data > class CContext
 {
 public:
-  CContext();
+  CContext(const bool & parallel = true);
   CContext(const CContext & src);
   ~CContext();
   CContext & operator = (const CContext & rhs);
@@ -59,18 +59,21 @@ template < class Data > class CContext
 protected:
   Data * mMasterData;
   Data * mThreadData;
+  bool mParallel;
   size_t mSize;
 };
 
-template < class Data > CContext< Data >::CContext()
+template < class Data > CContext< Data >::CContext(const bool & parallel)
   : mMasterData(NULL)
   , mThreadData(NULL)
+  , mParallel(parallel)
   , mSize(0)
 {}
 
 template < class Data > CContext< Data >::CContext(const CContext & src)
   : mMasterData(NULL)
   , mThreadData(NULL)
+  , mParallel(src.mParallel)
   , mSize(0)
 {
   init();
@@ -137,7 +140,7 @@ template < class Data > void CContext< Data >::init()
     return;
 
   mMasterData = new Data();
-  mSize = omp_get_max_threads();
+  mSize = mParallel ? omp_get_max_threads() : 1;
 
   if (mSize > 1)
     {

copasi/OpenMP/CPointerContext.h

@@ -17,7 +17,7 @@ template < class Data > class CPointerContext : public CContext< Data * >
 
   CPointerContext(const CPointerContext & src) = delete;
 
-  CPointerContext(Data * pMaster);
+  CPointerContext(const bool & parallel);
 
   ~CPointerContext();
 
@@ -28,9 +28,10 @@ template < class Data > class CPointerContext : public CContext< Data * >
   void setMaster(Data * pMaster);
 };
 
-template < class Data > CPointerContext< Data >::CPointerContext(Data * pMaster)
-  : Base()
+template < class Data > CPointerContext< Data >::CPointerContext(const bool & parallel)
+  : Base(parallel)
 {
+
   Base::init();
   Base::master() = NULL;
 
@@ -42,8 +43,6 @@ template < class Data > CPointerContext< Data >::CPointerContext(Data * pMaster)
       for (; pIt != pEnd; ++pIt)
         *pIt = NULL;
     }
-
-  setMaster(pMaster);
 }
 
 template < class Data > CPointerContext< Data >::~CPointerContext()

copasi/OpenMP/CPointerMathContext.h

@@ -17,7 +17,7 @@ template < class Data > class CPointerMathContext : public CPointerContext< Data
 
   CPointerMathContext(const CPointerMathContext & src) = delete;
 
-  CPointerMathContext(Data * pMaster);
+  CPointerMathContext(const bool & parallel);
 
   ~CPointerMathContext();
 
@@ -28,8 +28,8 @@ template < class Data > class CPointerMathContext : public CPointerContext< Data
   void setMathContext(CMathContext & Context);
 };
 
-template < class Data > CPointerMathContext< Data >::CPointerMathContext(Data * pMaster)
-  : CPointerContext< Data >(pMaster)
+template < class Data > CPointerMathContext< Data >::CPointerMathContext(const bool & parallel)
+  : CPointerContext< Data >(parallel)
 {}
 
 template < class Data > CPointerMathContext< Data >::~CPointerMathContext()

copasi/OpenMP/CRandomContext.h

@@ -0,0 +1,14 @@
+// Copyright (C) 2021 by Pedro Mendes, Rector and Visitors of the
+// University of Virginia, University of Heidelberg, and University
+// of Connecticut School of Medicine.
+// All rights reserved.
+
+#ifndef COPASI_CRandomContext
+#define COPASI_CRandomContext
+
+#include "copasi/OpenMP/CPointerContext.h"
+#include "copasi/randomGenerator/CRandom.h"
+
+typedef CPointerContext< CRandom > CRandomContext;
+
+#endif // COPASI_CRandomContext

copasi/optimization/COptMethod.cpp

@@ -43,13 +43,15 @@
 
 COptMethod::COptMethod(const CDataContainer * pParent,
                        const CTaskEnum::Method & methodType,
-                       const CTaskEnum::Task & taskType)
+                       const CTaskEnum::Task & taskType,
+                       const bool & parallel)
   : CCopasiMethod(pParent, methodType, taskType)
   , mpParentTask(NULL)
+  , mParallel(parallel)
+  , mMathContext(parallel)
+  , mProblemContext(parallel)
   , mLogVerbosity(0)
   , mMethodLog()
-  , mProblemContext(NULL)
-  , mMathContext(NULL)
 {
   assertParameter("Log Verbosity", CCopasiParameter::Type::UINT, (unsigned C_INT32) 0, eUserInterfaceFlag::editable);
 }
@@ -58,18 +60,42 @@ COptMethod::COptMethod(const COptMethod & src,
                        const CDataContainer * pParent)
   : CCopasiMethod(src, pParent)
   , mpParentTask(src.mpParentTask)
+  , mParallel(src.mParallel)
+  , mMathContext(src.mParallel)
+  , mProblemContext(src.mParallel)
   , mLogVerbosity(src.mLogVerbosity)
   , mMethodLog(src.mMethodLog)
-  , mProblemContext(src.mProblemContext.master())
-  , mMathContext(src.mMathContext.master())
 {
+  mMathContext.setMaster(src.mMathContext.master());
+  mProblemContext.setMaster(src.mProblemContext.master());
   mProblemContext.setMathContext(mMathContext);
 }
 
 //YOHE: seems "virtual" cannot be outside of class declaration
 COptMethod::~COptMethod()
 {}
 
+// static
+std::pair< C_FLOAT64, bool > COptMethod::objectiveValue(COptProblem * pProblem, const CVectorCore< C_FLOAT64 > & parameters)
+{
+  std::pair< C_FLOAT64, bool > Result;
+
+  pProblem->setParameters(parameters);
+  Result.second = pProblem->calculate();
+  Result.first = pProblem->getCalculateValue();
+
+  return Result;
+}
+
+// static
+void COptMethod::reflect(COptProblem * pProblem, const C_FLOAT64 & bestValue, C_FLOAT64 & objectiveValue)
+{
+  if (objectiveValue < bestValue
+      && (!pProblem->checkParametricConstraints()
+          || !pProblem->checkFunctionalConstraints()))
+    objectiveValue = bestValue + bestValue - objectiveValue;
+}
+
 void COptMethod::setProblem(COptProblem * pProblem)
 {
   mProblemContext.setMaster(pProblem);
@@ -87,6 +113,13 @@ bool COptMethod::initialize()
   if (mProblemContext.master() == NULL)
     return false;
 
+  COptProblem **ppProblem = mProblemContext.beginThread();
+  COptProblem **ppProblemEnd = mProblemContext.endThread();
+
+  for (; ppProblem != ppProblemEnd; ++ppProblem)
+    if (mProblemContext.isThread(ppProblem))
+      (*ppProblem)->initialize();
+
   mpParentTask = dynamic_cast<COptTask *>(getObjectParent());
 
   if (!mpParentTask) return false;

copasi/optimization/COptMethod.h

@@ -67,6 +67,15 @@ class COptMethod : public CCopasiMethod
 
   //data member
 protected:
+  /**
+   * The parent task
+   */
+  COptTask * mpParentTask;
+
+  /**
+   * Boolean indicating whether this method con use parallel execution.
+   */
+  bool mParallel;
 
   /**
    * A thread specific math container
@@ -78,8 +87,8 @@ class COptMethod : public CCopasiMethod
    */
   COptProblemContext mProblemContext;
 
-protected:
-  COptTask * mpParentTask;
+  /**
+   * A best objective value found
 
   /**
    * Define the current verbosity for the log
@@ -103,11 +112,13 @@ class COptMethod : public CCopasiMethod
    * Specific constructor
    * @param const CDataContainer * pParent
    * @param const CTaskEnum::Method & methodType
-   * @param const CTaskEnum::Task & taskType (default: optimization)
+   * @param const CTaskEnum::Task & taskType
+   * @param const bool & parallel
    */
   COptMethod(const CDataContainer * pParent,
              const CTaskEnum::Method & methodType,
-             const CTaskEnum::Task & taskType = CTaskEnum::Task::optimization);
+             const CTaskEnum::Task & taskType,
+             const bool & parallel);
 
   /**
    * Copy constructor
@@ -159,6 +170,22 @@ class COptMethod : public CCopasiMethod
   const COptLog &getMethodLog() const;
 
 protected:
+  /**
+   * Calculate the objective value for the provided parameter set
+   * @param COptProblem * pProblem
+   * @param const CVectorCore< C_FLOAT64 > & parameters
+   * @return std::pair< C_FLOAT64 objectiveValue, bool continue >
+   */
+  static std::pair< C_FLOAT64, bool > objectiveValue(COptProblem * pProblem, const CVectorCore< C_FLOAT64 > & parameters);
+
+  /**
+   * Reflect the objective value if it is outside the parametric or functional domain
+   * @param COptProblem * pProblem
+   * @param const C_FLOAT64 & bestValue
+   * @param C_FLOAT64 & objectiveValue
+   */
+  static void reflect(COptProblem * pProblem, const C_FLOAT64 & bestValue, C_FLOAT64 & objectiveValue);
+
   /**
    * Signal that the math container has changed
    */

copasi/optimization/COptMethodCoranaWalk.cpp

@@ -37,7 +37,7 @@
 COptMethodCoranaWalk::COptMethodCoranaWalk(const CDataContainer * pParent,
     const CTaskEnum::Method & methodType,
     const CTaskEnum::Task & taskType)
-  : COptMethod(pParent, methodType, taskType)
+  : COptMethod(pParent, methodType, taskType, false)
   , mTemperature(1.0)
   , mhIterations(C_INVALID_INDEX)
   , mIterations(100)
@@ -113,7 +113,7 @@ bool COptMethodCoranaWalk::optimise()
     }
 
   if (minstep > std::numeric_limits< C_FLOAT64 >::epsilon())
-    minstep /= i;
+    minstep /= mVariableSize;
   else
     minstep = 100 * std::numeric_limits< C_FLOAT64 >::epsilon();
 

copasi/optimization/COptMethodDE.cpp

@@ -34,7 +34,7 @@
 COptMethodDE::COptMethodDE(const CDataContainer * pParent,
                            const CTaskEnum::Method & methodType,
                            const CTaskEnum::Task & taskType):
-  COptPopulationMethod(pParent, methodType, taskType),
+  COptPopulationMethod(pParent, methodType, taskType, false),
   mpPermutation(NULL),
   mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
   mMutationVarians(0.1),
@@ -143,12 +143,12 @@ bool COptMethodDE::replicate()
           const COptItem & OptItem = *mProblemContext.master()->getOptItemList()[j];
           C_FLOAT64 & mut = (*mIndividuals[i])[j];
 
-          size_t r = mpRandom->getRandomU(mPopulationSize - 1);
+          size_t r = mRandomContext.master()->getRandomU(mPopulationSize - 1);
 
           if (r < 0.6 * mPopulationSize)
             {
               mut = (*mIndividuals[i - mPopulationSize])[j] *
-                    mpRandom->getRandomNormal(1, mMutationVarians);
+                    mRandomContext.master()->getRandomNormal(1, mMutationVarians);
             }
 
           else
@@ -193,10 +193,10 @@ bool COptMethodDE::replicate()
               const COptItem & OptItem = *mProblemContext.master()->getOptItemList()[j];
               C_FLOAT64 & mut = (*mIndividuals[i - 2 * mPopulationSize])[j];
 
-              size_t r = mpRandom->getRandomU(mPopulationSize - 1);
+              size_t r = mRandomContext.master()->getRandomU(mPopulationSize - 1);
 
               if (r < 0.6 * mPopulationSize)
-                mut *= mpRandom->getRandomNormal(1, mMutationVarians);
+                mut *= mRandomContext.master()->getRandomNormal(1, mMutationVarians);
 
               // force it to be within the bounds
               switch (OptItem.checkConstraint(mut))
@@ -279,15 +279,15 @@ bool COptMethodDE::creation(size_t first, size_t last)
               {
                 // determine if linear or log scale
                 if ((mn < 0.0) || (mx <= 0.0))
-                  mut = mn + mpRandom->getRandomCC() * (mx - mn);
+                  mut = mn + mRandomContext.master()->getRandomCC() * (mx - mn);
                 else
                   {
                     la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
 
                     if (la < 1.8)
-                      mut = mn + mpRandom->getRandomCC() * (mx - mn);
+                      mut = mn + mRandomContext.master()->getRandomCC() * (mx - mn);
                     else
-                      mut = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mpRandom->getRandomCC());
+                      mut = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min())) + la * mRandomContext.master()->getRandomCC());
                   }
               }
 
@@ -348,7 +348,7 @@ bool COptMethodDE::initialize()
       setValue("Population Size", mPopulationSize);
     }
 
-  mpPermutation = new CPermutation(mpRandom, mPopulationSize);
+  mpPermutation = new CPermutation(mRandomContext.master(), mPopulationSize);
 
   mIndividuals.resize(3 * mPopulationSize);