#include <COptMethodSS.h>

Inheritance diagram for COptMethodSS:

Collaboration diagram for COptMethodSS:

Public Member Functions
	COptMethodSS (const COptMethodSS &src, const CCopasiContainer *pParent=NULL)

virtual bool	elevateChildren ()

virtual bool	optimise ()

virtual	~COptMethodSS ()

Public Member Functions inherited from COptMethod
	COptMethod (const COptMethod &src, const CCopasiContainer *pParent=NULL)

bool	isBounded (void)

virtual bool	isValidProblem (const CCopasiProblem *pProblem)

void	setProblem (COptProblem *problem)

virtual	~COptMethod ()

Public Member Functions inherited from CCopasiMethod
	CCopasiMethod (const CCopasiMethod &src, const CCopasiContainer *pParent=NULL)

const CCopasiMethod::SubType &	getSubType () const

const CCopasiTask::Type &	getType () const

virtual void	load (CReadConfig &configBuffer, CReadConfig::Mode mode=CReadConfig::SEARCH)

virtual void	print (std::ostream *ostream) const

virtual void	printResult (std::ostream *ostream) const

virtual bool	setCallBack (CProcessReport *pCallBack)

virtual	~CCopasiMethod ()

Public Member Functions inherited from CCopasiParameterGroup
bool	addGroup (const std::string &name)

bool	addParameter (const CCopasiParameter &parameter)

bool	addParameter (const std::string &name, const CCopasiParameter::Type type)

template<class CType >
bool	addParameter (const std::string &name, const CCopasiParameter::Type type, const CType &value)

void	addParameter (CCopasiParameter *pParameter)

CCopasiParameterGroup *	assertGroup (const std::string &name)

template<class CType >
CCopasiParameter *	assertParameter (const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)

index_iterator	beginIndex () const

name_iterator	beginName () const

	CCopasiParameterGroup (const CCopasiParameterGroup &src, const CCopasiContainer *pParent=NULL)

	CCopasiParameterGroup (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &objectType="ParameterGroup")

void	clear ()

index_iterator	endIndex () const

name_iterator	endName () const

CCopasiParameterGroup *	getGroup (const std::string &name)

const CCopasiParameterGroup *	getGroup (const std::string &name) const

CCopasiParameterGroup *	getGroup (const size_t &index)

const CCopasiParameterGroup *	getGroup (const size_t &index) const

size_t	getIndex (const std::string &name) const

std::string	getKey (const std::string &name) const

std::string	getKey (const size_t &index) const

virtual const std::string &	getName (const size_t &index) const

virtual const CObjectInterface *	getObject (const CCopasiObjectName &cn) const

CCopasiParameter *	getParameter (const std::string &name)

const CCopasiParameter *	getParameter (const std::string &name) const

CCopasiParameter *	getParameter (const size_t &index)

const CCopasiParameter *	getParameter (const size_t &index) const

CCopasiParameter::Type	getType (const std::string &name) const

CCopasiParameter::Type	getType (const size_t &index) const

std::string	getUniqueParameterName (const CCopasiParameter *pParameter) const

const CCopasiParameter::Value &	getValue (const std::string &name) const

const CCopasiParameter::Value &	getValue (const size_t &index) const

CCopasiParameter::Value &	getValue (const std::string &name)

CCopasiParameter::Value &	getValue (const size_t &index)

CCopasiParameterGroup &	operator= (const CCopasiParameterGroup &rhs)

bool	removeParameter (const std::string &name)

bool	removeParameter (const size_t &index)

template<class CType >
bool	setValue (const std::string &name, const CType &value)

template<class CType >
bool	setValue (const size_t &index, const CType &value)

size_t	size () const

bool	swap (const size_t &iFrom, const size_t &iTo)

bool	swap (index_iterator &from, index_iterator &to)

virtual	~CCopasiParameterGroup ()

Public Member Functions inherited from CCopasiParameter
	CCopasiParameter (const CCopasiParameter &src, const CCopasiContainer *pParent=NULL)

	CCopasiParameter (const std::string &name, const Type &type, const void pValue=NULL, const CCopasiContainer pParent=NULL, const std::string &objectType="Parameter")

virtual CCopasiObjectName	getCN () const

virtual const std::string &	getKey () const

virtual std::string	getObjectDisplayName (bool regular=true, bool richtext=false) const

const CCopasiParameter::Type &	getType () const

const Value &	getValue () const

Value &	getValue ()

virtual void *	getValuePointer () const

CCopasiObject *	getValueReference () const

bool	isValidValue (const C_FLOAT64 &value) const

bool	isValidValue (const C_INT32 &value) const

bool	isValidValue (const unsigned C_INT32 &value) const

bool	isValidValue (const bool &value) const

bool	isValidValue (const std::string &value) const

bool	isValidValue (const CCopasiObjectName &value) const

bool	isValidValue (const std::vector< CCopasiParameter * > &value) const

CCopasiParameter &	operator= (const CCopasiParameter &rhs)

template<class CType >
bool	setValue (const CType &value)

bool	setValue (const std::vector< CCopasiParameter * > &value)

virtual	~CCopasiParameter ()

Public Member Functions inherited from CCopasiContainer
virtual bool	add (CCopasiObject *pObject, const bool &adopt=true)

	CCopasiContainer (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=CCopasiObject::Container)

	CCopasiContainer (const CCopasiContainer &src, const CCopasiContainer *pParent=NULL)

virtual std::string	getChildObjectUnits (const CCopasiObject *pObject) const

virtual const objectMap &	getObjects () const

virtual std::string	getUnits () const

virtual const CCopasiObject *	getValueObject () const

virtual bool	remove (CCopasiObject *pObject)

virtual	~CCopasiContainer ()

Public Member Functions inherited from CCopasiObject
void	addDirectDependency (const CCopasiObject *pObject)

	CCopasiObject (const CCopasiObject &src, const CCopasiContainer *pParent=NULL)

void	clearDirectDependencies ()

void	clearRefresh ()

bool	dependsOn (DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const

void	getAllDependencies (DataObjectSet &dependencies, const DataObjectSet &context) const

virtual const DataObjectSet &	getDirectDependencies (const DataObjectSet &context=DataObjectSet()) const

CCopasiContainer *	getObjectAncestor (const std::string &type) const

CCopasiDataModel *	getObjectDataModel ()

const CCopasiDataModel *	getObjectDataModel () const

const std::string &	getObjectName () const

CCopasiContainer *	getObjectParent () const

const std::string &	getObjectType () const

virtual const CObjectInterface::ObjectSet &	getPrerequisites () const

virtual Refresh *	getRefresh () const

UpdateMethod *	getUpdateMethod () const

bool	hasCircularDependencies (DataObjectSet &candidates, DataObjectSet &verified, const DataObjectSet &context) const

bool	hasUpdateMethod () const

bool	isArray () const

bool	isContainer () const

bool	isDataModel () const

bool	isMatrix () const

bool	isNameVector () const

bool	isNonUniqueName () const

virtual bool	isPrerequisiteForContext (const CObjectInterface *pObject, const CMath::SimulationContextFlag &context, const CObjectInterface::ObjectSet &changedObjects) const

bool	isReference () const

bool	isRoot () const

bool	isSeparator () const

bool	isStaticString () const

bool	isValueBool () const

bool	isValueDbl () const

bool	isValueInt () const

bool	isValueInt64 () const

bool	isValueString () const

bool	isVector () const

virtual bool	mustBeDeleted (const DataObjectSet &deletedObjects) const

void	removeDirectDependency (const CCopasiObject *pObject)

void	setDirectDependencies (const DataObjectSet &directDependencies)

bool	setObjectName (const std::string &name)

virtual bool	setObjectParent (const CCopasiContainer *pParent)

void	setObjectValue (const C_FLOAT64 &value)

void	setObjectValue (const C_INT32 &value)

void	setObjectValue (const bool &value)

template<class CType >
void	setRefresh (CType pType, void(CType::method)(void))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_FLOAT64 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_INT32 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const bool &))

virtual	~CCopasiObject ()

Public Member Functions inherited from CObjectInterface
	CObjectInterface ()

virtual	~CObjectInterface ()

Private Member Functions
bool	childLocalMin (void)

virtual bool	cleanup ()

bool	closerChild (C_INT32 i, C_INT32 j, C_FLOAT64 dist)

bool	closerRefSet (C_INT32 i, C_INT32 j, C_FLOAT64 dist)

bool	combination (void)

	COptMethodSS (const CCopasiContainer *pParent=NULL)

bool	creation (void)

bool	evaluate (const CVector< C_FLOAT64 > &individual)

bool	inforefset (C_INT32 type, C_INT32 element)

virtual bool	initialize ()

void	initObjects ()

bool	localmin (CVector< C_FLOAT64 > &solution, C_FLOAT64 &fval)

bool	randomize (C_INT32 i)

bool	serializechildren (size_t first, size_t last)

bool	serializepool (size_t first, size_t last)

bool	serializeprob (void)

bool	serializerefset (C_INT32 first, C_INT32 last)

bool	serializevector (CVector< C_FLOAT64 > x, C_FLOAT64 xval)

void	sortRefSet (C_INT32 lower, C_INT32 upper)

Private Attributes
size_t	mBestIndex

C_FLOAT64	mBestValue

std::vector< CVector < C_FLOAT64 > * >	mChild

bool	mChildrenGenerated

CVector< C_FLOAT64 >	mChildVal

C_FLOAT64	mCloseValue

C_FLOAT64	mEvaluationValue

std::vector< CVector< C_INT32 > * >	mFreq

size_t	mhIterations

unsigned C_INT32	mIteration

unsigned C_INT32	mIterations

unsigned C_INT32	mLocalFreq

unsigned C_INT32	mLocalIter

unsigned C_INT32	mLocalStored

COptMethod *	mpLocalMinimizer

std::vector< CVector < C_FLOAT64 > * >	mPool

size_t	mPoolSize

CVector< C_FLOAT64 >	mPoolVal

COptProblem *	mpOptProblemLocal

C_INT32	mPopulationSize

CRandom *	mpRandom

CVector< C_FLOAT64 >	mProb

std::vector< CVector < C_FLOAT64 > * >	mRefSet

CVector< C_FLOAT64 >	mRefSetVal

CVector< C_INT32 >	mStuck

unsigned C_INT32	mVariableSize

Friends
COptMethod *	COptMethod::createMethod (CCopasiMethod::SubType subType)

Additional Inherited Members
Public Types inherited from CCopasiMethod
enum	SubType { unset = 0, RandomSearch, RandomSearchMaster, SimulatedAnnealing, CoranaWalk, DifferentialEvolution, ScatterSearch, GeneticAlgorithm, EvolutionaryProgram, SteepestDescent, HybridGASA, GeneticAlgorithmSR, HookeJeeves, LevenbergMarquardt, NelderMead, SRES, Statistics, ParticleSwarm, Praxis, TruncatedNewton, Newton, deterministic, LSODAR, directMethod, stochastic, tauLeap, adaptiveSA, hybrid, hybridLSODA, hybridODE45, DsaLsodar, tssILDM, tssILDMModified, tssCSP, mcaMethodReder, scanMethod, lyapWolf, sensMethod, EFMAlgorithm, EFMBitPatternTreeAlgorithm, EFMBitPatternAlgorithm, Householder, crossSectionMethod, linearNoiseApproximation }

Public Types inherited from CCopasiParameterGroup
typedef parameterGroup::iterator	index_iterator

typedef CCopasiContainer::objectMap::iterator	name_iterator

typedef std::vector < CCopasiParameter * >	parameterGroup

Public Types inherited from CCopasiParameter
enum	Type { DOUBLE = 0, UDOUBLE, INT, UINT, BOOL, GROUP, STRING, CN, KEY, FILE, EXPRESSION, INVALID }

Public Types inherited from CCopasiContainer
typedef std::multimap < std::string, CCopasiObject * >	objectMap

Public Types inherited from CCopasiObject
typedef std::set< const CCopasiObject * >	DataObjectSet

typedef std::vector< Refresh * >	DataUpdateSequence

Public Types inherited from CObjectInterface
typedef std::set< const CObjectInterface * >	ObjectSet

typedef std::vector < CObjectInterface * >	UpdateSequence

Static Public Member Functions inherited from COptMethod
static COptMethod *	createMethod (CCopasiMethod::SubType subType=CCopasiMethod::RandomSearch)

Static Public Member Functions inherited from CCopasiObject
static std::vector< Refresh * >	buildUpdateSequence (const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())

static void	setRenameHandler (CRenameHandler *rh)

Static Public Attributes inherited from CCopasiMethod
static const std::string	SubTypeName []

static const char *	XMLSubType []

Static Public Attributes inherited from CCopasiParameter
static const std::string	TypeName []

static const char *	XMLType []

Static Public Attributes inherited from CCopasiContainer
static const std::vector < CCopasiContainer * >	EmptyList

Protected Types inherited from CCopasiObject
enum	Flag { Container = 0x1, Vector = 0x2, Matrix = 0x4, NameVector = 0x8, Reference = 0x10, ValueBool = 0x20, ValueInt = 0x40, ValueInt64 = 0x80, ValueDbl = 0x100, NonUniqueName = 0x200, StaticString = 0x400, ValueString = 0x800, Separator = 0x1000, ModelEntity = 0x2000, Array = 0x4000, DataModel = 0x8000, Root = 0x10000, Gui = 0x20000 }

Protected Member Functions inherited from COptMethod
	COptMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)

Protected Member Functions inherited from CCopasiMethod
	CCopasiMethod (const CCopasiTask::Type &taskType, const SubType &subType, const CCopasiContainer *pParent=NULL)

Protected Member Functions inherited from CCopasiParameterGroup
	CCopasiParameterGroup ()

Protected Member Functions inherited from CCopasiContainer
template<class CType >
CCopasiObject *	addMatrixReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addObjectReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addVectorReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

void	initObjects ()

Protected Member Functions inherited from CCopasiObject
	CCopasiObject ()

	CCopasiObject (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=0)

Protected Attributes inherited from COptMethod
const bool	mBounds

const std::vector< COptItem * > *	mpOptContraints

const std::vector< COptItem * > *	mpOptItem

COptProblem *	mpOptProblem

COptTask *	mpParentTask

const std::vector < UpdateMethod * > *	mpSetCalculateVariable

Protected Attributes inherited from CCopasiMethod
CProcessReport *	mpCallBack

Protected Attributes inherited from CCopasiParameter
std::string	mKey

CCopasiObject *	mpValueReference

size_t	mSize

Value	mValue

Protected Attributes inherited from CCopasiContainer
objectMap	mObjects

Static Protected Attributes inherited from CCopasiObject
static CRenameHandler *	smpRenameHandler = NULL

Detailed Description

Definition at line 20 of file COptMethodSS.h.

Constructor & Destructor Documentation

COptMethodSS::COptMethodSS	(	const COptMethodSS &	src,
		const CCopasiContainer *	pParent = `NULL`
	)

Copy Constructor

Parameters

const	COptMethodSS & src
const	CCopasiContainer * pParent (default: NULL)

Definition at line 59 of file COptMethodSS.cpp.

References initObjects().

                                                             :
   COptMethod(src, pParent),
   mIterations(0),
   mPopulationSize(0),
   mVariableSize(0),
   mRefSet(0),
   mRefSetVal(0),
   mPool(0),
   mPoolVal(0),
   mPoolSize(0),
   mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
   mIteration(0),
   mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
   mBestIndex(C_INVALID_INDEX),
   mpRandom(NULL),
   mpOptProblemLocal(NULL),
   mpLocalMinimizer(NULL)
 {
   // remove eventual existing parameters from the release version.
   initObjects();
 }

COptMethodSS::~COptMethodSS ( )

virtual

Destructor

Definition at line 82 of file COptMethodSS.cpp.

References cleanup().

83 {cleanup();}

COptMethodSS::cleanup

virtual bool cleanup()

Definition: COptMethodSS.cpp:267

COptMethodSS::COptMethodSS ( const CCopasiContainer * pParent = NULL )

private

Default Constructor

Parameters

const CCopasiContainer * pParent (default: NULL)

Definition at line 30 of file COptMethodSS.cpp.

References CCopasiParameterGroup::addParameter(), C_INT32, initObjects(), CRandom::mt19937, and CCopasiParameter::UINT.

                                                           :
   COptMethod(CCopasiTask::optimization, CCopasiMethod::ScatterSearch, pParent),
   mIterations(0),
   mPopulationSize(0),
   mVariableSize(0),
   mRefSet(0),
   mRefSetVal(0),
   mPool(0),
   mPoolVal(0),
   mPoolSize(0),
   mEvaluationValue(std::numeric_limits< C_FLOAT64 >::max()),
   mIteration(0),
   mBestValue(std::numeric_limits< C_FLOAT64 >::max()),
   mBestIndex(C_INVALID_INDEX),
   mpRandom(NULL),
   mpOptProblemLocal(NULL),
   mpLocalMinimizer(NULL)
 {
   addParameter("Number of Iterations", CCopasiParameter::UINT, (unsigned C_INT32) 200);
 // we no longer give the user choice of rng, we use the mersenne twister!
 // but in DEBUG versions we should still have access to it
 #ifdef COPASI_DEBUG
   addParameter("Random Number Generator", CCopasiParameter::UINT, (unsigned C_INT32) CRandom::mt19937);
   addParameter("Seed", CCopasiParameter::UINT, (unsigned C_INT32) 0);
 #endif
 
   initObjects();
 }

Member Function Documentation

bool COptMethodSS::childLocalMin ( void )

private

minimize the best child if within criteria

Returns: bool continue

Definition at line 1017 of file COptMethodSS.cpp.

References C_FLOAT64, C_INT32, closerChild(), localmin(), mChild, mChildVal, mCloseValue, mLocalIter, mLocalStored, mPool, mPoolVal, mPopulationSize, and mStuck.

Referenced by optimise().

 {
   C_INT32 i, best;
   C_FLOAT64 bestVal, fvalmin;
   bool Running = true;
 
   // signal nothing found yet
   best = -1;
   bestVal = std::numeric_limits<C_FLOAT64>::infinity();
 
   // find the best child
   for (i = 0; i < mPopulationSize; i++)
     {
       if ((mStuck[i] == 0) && (mChildVal[i] < bestVal))
         {
           bestVal = mChildVal[i];
           best = i;
         }
     }
 
   // no child in this iteration? exit now
   if (best == -1) return true;
 
   // check if this child is not close to previous ones
   for (i = 0; i < mLocalStored; i++)
     {
       // is the other one like me?
       if (closerChild(best, i, mCloseValue))
         {
           // it is too close, exit now
           return true;
         }
     }
 
   // store the initial position
   *(mPool[mLocalStored]) = *(mChild[best]);
   mPoolVal[mLocalStored] = mChildVal[best]; //bestVal;
   mLocalStored++;
 
   // do local minimization on it
   Running &= localmin(*(mChild[best]), mChildVal[best]);
 
   // store the result
   *(mPool[mLocalStored]) = *(mChild[best]);
   mPoolVal[mLocalStored] = mChildVal[best];
   mLocalStored++;
 
   // clear the local optimisation counter
   mLocalIter = 1;
 
   return Running;
 }

bool COptMethodSS::cleanup ( )

privatevirtual

Cleanup arrays and pointers.

Returns: bool success

Reimplemented from COptMethod.

Definition at line 267 of file COptMethodSS.cpp.

References mChild, mFreq, mpLocalMinimizer, mPool, mpOptProblemLocal, mpRandom, mRefSet, and pdelete.

Referenced by initialize(), optimise(), and ~COptMethodSS().

 {
   size_t i;
 
   pdelete(mpRandom);
 
   pdelete(mpOptProblemLocal);
 
   pdelete(mpLocalMinimizer);
 
   for (i = 0; i < mRefSet.size(); i++)
     pdelete(mRefSet[i]);
 
   for (i = 0; i < mChild.size(); i++)
     pdelete(mChild[i]);
 
   for (i = 0; i < mPool.size(); i++)
     pdelete(mPool[i]);
 
   for (i = 0; i < mFreq.size(); i++)
     pdelete(mFreq[i]);
 
   return true;
 }

bool COptMethodSS::closerChild	(	C_INT32	i,
		C_INT32	j,
		C_FLOAT64	dist
	)

private

check if two child vectors are close (one vector in mChild, the other in mPool, where previous child are stored)

Definition at line 799 of file COptMethodSS.cpp.

References C_FLOAT64, C_INT32, mChild, mPool, and mVariableSize.

Referenced by childLocalMin().

 {
   C_FLOAT64 mx;
 
   for (C_INT32 k = 0; k < mVariableSize; k++)
     {
       mx = (fabs((*mChild[i])[k]) + fabs((*mPool[j])[k])) / 2.0;
 
       if (fabs((*mChild[i])[k] - (*mPool[j])[k]) / mx > dist) return false;
     }
 
   return true;
 }

bool COptMethodSS::closerRefSet	(	C_INT32	i,
		C_INT32	j,
		C_FLOAT64	dist
	)

private

check if two vectors in RefSet are close

Definition at line 815 of file COptMethodSS.cpp.

References C_FLOAT64, C_INT32, mRefSet, and mVariableSize.

Referenced by optimise().

 {
   C_FLOAT64 mx;
 
   for (C_INT32 k = 0; k < mVariableSize; k++)
     {
       mx = (fabs((*mRefSet[i])[k]) + fabs((*mRefSet[j])[k])) / 2.0;
 
       if (fabs((*mRefSet[i])[k] - (*mRefSet[j])[k]) / mx > dist) return false;
     }
 
   return true;
 }

bool COptMethodSS::combination ( void )

private

create new individuals with the combination method

Returns: bool continue

Definition at line 831 of file COptMethodSS.cpp.

References C_FLOAT64, C_INT32, COptItem::checkConstraint(), evaluate(), COptItem::getLowerBoundValue(), CRandom::getRandomCC(), COptItem::getUpperBoundValue(), max, mChild, mChildrenGenerated, mChildVal, mEvaluationValue, min, mPopulationSize, mpRandom, mRefSet, mRefSetVal, mStuck, mVariableSize, and CVector< CType >::resize().

Referenced by optimise().

 {
   C_INT32 i, j, k, l;   // counters
   C_FLOAT64 mn, mx;     // for bounds on parameters
   C_FLOAT64 beta;       // bias
   C_FLOAT64 la;         // for orders of magnitude
   C_FLOAT64 alpha;      // 1 or -1
   C_FLOAT64 bm2;        // b-2
   C_FLOAT64 omatb;      // (1-alpha*beta)*0.5
   C_FLOAT64 dd;         // (x(i) - x(j) ) / 2 * (1-alpha*beta)
   C_FLOAT64 c1, c2;     // coordinates of the edges of hyperectangle
   CVector< C_FLOAT64 > xnew, xpr;
   C_FLOAT64 xnewval, xprval; // to hold temp value of "parent" in go-beyond strategy
   C_FLOAT64 lambda = 1.0; // damping factor for go-beyond strategy
   C_INT32 improvement; // count iterations with improvement in go-beyond strategy
   bool Running = true; // flag for invalid values
 
   // make xnew large enough
   xnew.resize(mVariableSize);
   xpr.resize(mVariableSize);
   // calculate a constant
   bm2 = C_FLOAT64(mPopulationSize) - 2.0;
   // signal no children yet
   mChildrenGenerated = false;
 
   // generate children for each member of the population
   for (i = 0; (i < mPopulationSize) && Running; i++)
     {
       // keep the parent value in childval[i] so that we only accept better than that
       mChildVal[i] = mRefSetVal[i];
 
       for (j = 0; j < mPopulationSize; j++)
         {
           // no self-reproduction...
           if (i != j)
             {
               if (i < j) alpha = 1.0; else alpha = -1.0;
 
               beta = (fabs(C_FLOAT64(j) - C_FLOAT64(i)) - 1.0) / bm2;
               omatb = (1.0 + alpha * beta) * 0.5;
 
               // generate a child
               for (k = 0; k < mVariableSize; k++)
                 {
                   // get the bounds of this parameter
                   COptItem & OptItem = *(*mpOptItem)[k];
                   mn = *OptItem.getLowerBoundValue();
                   mx = *OptItem.getUpperBoundValue();
 
                   try
                     {
                       // calculate orders of magnitude of the interval
                       if (((*mRefSet[i])[k] > 0.0) && ((*mRefSet[j])[k] > 0.0))
                         {
                           la = log10((*mRefSet[i])[k]) - log10(std::max((*mRefSet[j])[k], std::numeric_limits< C_FLOAT64 >::min()));
                         }
                       else
                         la = 1.0;
 
                       dd = ((*mRefSet[i])[k] - (*mRefSet[j])[k]) * omatb;
                       // one of the box limits
                       c1 = (*mRefSet[i])[k] - dd;
 
                       // force it to be within the bounds
                       switch (OptItem.checkConstraint(c1))
                         {
                           case -1:
                             c1 = mn;
                             break;
 
                           case 1:
                             c1 = mx;
                             break;
                         }
 
                       // the other box limit
                       c2 = (*mRefSet[i])[k] + dd;
 
                       // force it to be within the bounds
                       switch (OptItem.checkConstraint(c2))
                         {
                           case -1:
                             c2 = mn;
                             break;
 
                           case 1:
                             c2 = mx;
                             break;
                         }
 
                       xnew[k] = c1 + (c2 - c1) * mpRandom->getRandomCC();
                     }
                   catch (...)
                     {
                       // if something failed leave the value intact
                       xnew[k] = (*mRefSet[i])[k];
                     }
 
                   // We need to set the value here so that further checks take
                   // account of the value.
                   (*(*mpSetCalculateVariable)[k])(xnew[k]);
                 }
 
               // calculate the child's fitness
               Running &= evaluate(xnew);
 
               // keep it if it is better than the previous one
               if (mEvaluationValue < mChildVal[i])
                 {
                   // keep a copy of this vector in mChild
                   (*mChild[i]) = xnew;
                   // and the fitness value
                   mChildVal[i] = mEvaluationValue;
                   // signal that child is better than parent
                   mStuck[i] = 0;
                   // signal we have generated a child (improvement)
                   mChildrenGenerated = true;
                 }
             }
         }
 
       // now we apply the go-beyond strategy for
       // cases where the child was an improvement
       if (mStuck[i] == 0)
         {
           // copy the parent
           xpr = (*mRefSet[i]);
           xprval = mRefSetVal[i];
           lambda = 1.0; // this is really 1/lambda so we can use mult rather than div
           improvement = 1;
 
           // while newval < childval
           for (; ;)
             {
               for (k = 0; k < mVariableSize; k++)
                 {
                   dd = (xpr[i] - (*mChild[i])[k]) * lambda;
                   xnew[k] = (*mChild[i])[k] + dd * mpRandom->getRandomCC();
                   // get the bounds of this parameter
                   COptItem & OptItem = *(*mpOptItem)[k];
 
                   // put it on the bounds if it had exceeded them
                   switch (OptItem.checkConstraint(xnew[k]))
                     {
                       case -1:
                         xnew[k] = *OptItem.getLowerBoundValue();
                         break;
 
                       case 1:
                         xnew[k] = *OptItem.getUpperBoundValue();
                         break;
                     }
 
                   // We need to set the value here so that further checks take
                   // account of the value.
                   (*(*mpSetCalculateVariable)[k])(xnew[k]);
                 }
 
               // calculate the child's fitness
               Running &= evaluate(xnew);
               xnewval = mEvaluationValue;
 
               // if there was no improvement we finish here => exit for(;;)
               if (mChildVal[i] <= xnewval) break;
 
               // old child becomes parent
               xpr = (*mChild[i]);
               xprval = mChildVal[i];
               // new child becomes child
               (*mChild[i]) = xnew;
               mChildVal[i] = xnewval;
               // mark improvement
               improvement++;
 
               if (improvement == 2)
                 {
                   lambda *= 0.5;
                   improvement = 0;
                 }
             }
         }
     }
 
   return Running;
 }

bool COptMethodSS::creation ( void )

private

Initialise the population

Returns: bool continue

Definition at line 416 of file COptMethodSS.cpp.

References C_FLOAT64, C_INT32, COptItem::checkConstraint(), evaluate(), COptItem::getLowerBoundValue(), CRandom::getRandomCC(), COptItem::getStartValue(), COptItem::getUpperBoundValue(), max, mEvaluationValue, mFreq, min, mPool, mPoolSize, mPoolVal, mPopulationSize, mpRandom, mProb, mRefSet, mRefSetVal, mStuck, mVariableSize, and sortRefSet().

Referenced by optimise().

 {
   C_INT32 i, j, k, l;   // counters
   C_FLOAT64 mn, mx, la; // for boundaries of rnd
   C_FLOAT64 sum;        // to calculate a summation
   C_FLOAT64 a;          // to hold a random number
   bool Running = true;  // flag for invalid values
 
   // initialize all entries of the Pool.
   // first 4 candidates as a latin hypercube
   for (i = 0; (i < 4) && Running; i++)
     {
       for (j = 0; j < mVariableSize; j++)
         {
           // get pointers to appropriate elements (easier reading of code)
           COptItem & OptItem = *(*mpOptItem)[j];
           C_FLOAT64 & Sol = (*mPool[i])[j];
           // calculate lower and upper bounds for this variable
           mn = *OptItem.getLowerBoundValue();
           mx = *OptItem.getUpperBoundValue();
 
           try
             {
               // calculate orders of magnitude of the interval
               la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
 
               // determine if linear or log scale
               if ((mn < 0.0) || (mx <= 0.0))
                 Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) i) * (mx - mn) * 0.25;
               else
                 {
                   if (la < 1.8)
                     Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) i) * (mx - mn) * 0.25;
                   else
                     Sol = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()))
                               + la * 0.25 * (mpRandom->getRandomCC() + (C_FLOAT64) i));
                 }
             }
           catch (...)
             {
               // TODO: this sounds a bit daft in this context, what else could be done, though?
               Sol = (mx + mn) * 0.5;
             }
 
           // force it to be within the bounds
           switch (OptItem.checkConstraint(Sol))
             {
               case - 1:
                 Sol = *OptItem.getLowerBoundValue();
                 break;
 
               case 1:
                 Sol = *OptItem.getUpperBoundValue();
                 break;
             }
 
           // We need to set the value here so that further checks take
           // account of the value.
           (*(*mpSetCalculateVariable)[j])(Sol);
         }
 
       // calculate its fitness
       Running &= evaluate(*mPool[i]);
       mPoolVal[i] = mEvaluationValue;
     }
 
   // next we add the initial guess from the user
   for (j = 0; j < mVariableSize; j++)
     {
       COptItem & OptItem = *(*mpOptItem)[j];
       C_FLOAT64 & Sol = (*mPool[i])[j];
 
       // get the vector of initial value
       Sol = OptItem.getStartValue();
 
       // force it to be within the bounds
       switch (OptItem.checkConstraint(Sol))
         {
           case - 1:
             Sol = *OptItem.getLowerBoundValue();
             break;
 
           case 1:
             Sol = *OptItem.getUpperBoundValue();
             break;
         }
 
       // We need to set the value here so that further checks take
       // account of the value.
       (*(*mpSetCalculateVariable)[j])(Sol);
     }
 
   // calculate its fitness
   Running &= evaluate(*mPool[i]);
   mPoolVal[i] = mEvaluationValue;
   i++;
 
   // the remaining entries depend on probabilities
   for (; (i < mPoolSize) && Running; i++)
     {
       for (j = 0; j < mVariableSize; j++)
         {
           // get pointers to appropriate elements (easier reading of code)
           COptItem & OptItem = *(*mpOptItem)[j];
           C_FLOAT64 & Sol = (*mPool[i])[j];
           // calculate lower and upper bounds for this variable
           mn = *OptItem.getLowerBoundValue();
           mx = *OptItem.getUpperBoundValue();
 
           for (k = 0; k < 4; k++)
             {
               for (l = 0, sum = 0.0 ; l < 4; l++)
                 sum += 1.0 / (*mFreq[j])[l];
 
               mProb[k] = 1.0 / (*mFreq[j])[k] / sum;
 
               // we only want the cumulative probabilities
               if (k > 0) mProb[k] += mProb[k - 1];
             }
 
           a = mpRandom->getRandomCC();
 
           for (k = 0; k < 4; k++)
             {
               // note that the original is <= but numerically < is essentially the same and faster
               if (a < mProb[k])
                 {
                   try
                     {
                       // calculate orders of magnitude of the interval
                       la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
 
                       // determine if linear or log scale
                       if ((mn < 0.0) || (mx <= 0.0))
                         Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) k) * (mx - mn) * 0.25;
                       else
                         {
                           if (la < 1.8)
                             Sol = mn + (mpRandom->getRandomCC() + (C_FLOAT64) k) * (mx - mn) * 0.25;
                           else
                             Sol = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()))
                                       + la * 0.25 * (mpRandom->getRandomCC() + (C_FLOAT64) k));
                         }
                     }
                   catch (...)
                     {
                       // TODO: this sounds a bit daft in this context, what else could be done, though?
                       Sol = (mx + mn) * 0.5;
                     }
 
                   // force it to be within the bounds
                   switch (OptItem.checkConstraint(Sol))
                     {
                       case - 1:
                         Sol = *OptItem.getLowerBoundValue();
                         break;
 
                       case 1:
                         Sol = *OptItem.getUpperBoundValue();
                         break;
                     }
 
                   // We need to set the value here so that further checks take
                   // account of the value.
                   (*(*mpSetCalculateVariable)[j])(Sol);
                   // increase the frequency
                   (*mFreq[j])[k] += 1;
                   break;
                 }
             }
         }
 
       // calculate its fitness
       Running &= evaluate(*mPool[i]);
       mPoolVal[i] = mEvaluationValue;
     }
 
   // at this point the pool is formed
   // now we partially sort the pool for the h = b/2 top elements,
   // where b is mPopulationSize, This is done in place with heap sort
   CVector< C_FLOAT64 > *tempvec;
   C_FLOAT64 tempval;
   C_INT32 parent, child;
   C_INT32 h = mPopulationSize / 2;
 
   // top h are the heap, we first have to sort it
   for (i = 1; i < h; i++)
     {
       // bubble-up
       child = i;
 
       for (;;)
         {
           if (child == 0) break;
 
           parent = floor((double)(child - 1) / 2);
 
           if (mPoolVal[child] < mPoolVal[parent])
             {
               // swap with parent
               tempval = mPoolVal[child];
               mPoolVal[child] = mPoolVal[parent];
               mPoolVal[parent] = tempval;
               tempvec = mPool[child];
               mPool[child] = mPool[parent];
               mPool[parent] = tempvec;
               // make parent the new child
               child = parent;
             }
           else break;
         }
     }
 
   for (i = h; i < mPoolSize; i++)
     {
       child = 0;
 
       // check if this element is smaller than any of the leafs
       for (size_t leaf = h / 2; leaf < h; leaf++)
         {
           if (mPoolVal[i] < mPoolVal[leaf])
             {
               // keep if this leaf is worse than previous
               if (mPoolVal[child] < mPoolVal[leaf])
                 child = leaf;
             }
         }
 
       if (child == 0) continue;
 
       if (mPoolVal[i] < mPoolVal[child])
         {
           // swap i and h+j
           tempval = mPoolVal[child];
           mPoolVal[child] = mPoolVal[i];
           mPoolVal[i] = tempval;
           tempvec = mPool[child];
           mPool[child] = mPool[i];
           mPool[i] = tempvec;
 
           // now bubble-up
           for (;;)
             {
               if (child == 0) break;
 
               parent = floor((double)(child - 1) / 2);
 
               if (mPoolVal[child] < mPoolVal[parent])
                 {
                   // swap with parent
                   tempval = mPoolVal[child];
                   mPoolVal[child] = mPoolVal[parent];
                   mPoolVal[parent] = tempval;
                   tempvec = mPool[child];
                   mPool[child] = mPool[parent];
                   mPool[parent] = tempvec;
                   // make parent the new child
                   child = parent;
                 }
               else break;
             }
         }
     }
 
   // since some leafs are not in order in the array, we now do
   // a bubble sort (note: this is best case for bubble sort)
   // we use j because we do not want to change the value of h
   j = h;
 
   do
     {
       k = 0;
 
       for (i = 0; i <= j; i++)
         {
           if (mPoolVal[i] > mPoolVal[i + 1])
             {
               tempval = mPoolVal[i];
               mPoolVal[i] = mPoolVal[i + 1];
               mPoolVal[i + 1] = tempval;
               tempvec = mPool[i];
               mPool[i] = mPool[i + 1];
               mPool[i + 1] = tempvec;
               k = i;
             }
         }
 
       j = k;
     }
   while (j > 1);
 
   // at this point the pool is formed and partially sorted
   // now we build the RefSet:
   // 1st b/2 are the best ones in the Pool (sorted already)
   // the 2nd b/2 are random (or later can be made diverse by
   // maximising the Euclidean distances)
   for (i = 0; i < mPopulationSize; i++)
     {
       (*mRefSet[i]) = (*mPool[i]); // copy the whole vector
       mRefSetVal[i] = mPoolVal[i]; // keep the value
       mStuck[i] = 1;               // initialize the mStuck counter
     }
 
   // RefSet needs to be fully sorted. Note that top half is sorted
   // and we are garanteed that it contains the best values
   // so it is only bottom half that needs sorting.
   sortRefSet(h, mPopulationSize);
   // we're done
   return Running;
 }

bool COptMethodSS::elevateChildren ( )

virtual

This methods must be called to elevate subgroups to derived objects. The default implementation does nothing.

Returns: bool success

Reimplemented from CCopasiParameterGroup.

Definition at line 86 of file COptMethodSS.cpp.

References CCopasiParameterGroup::removeParameter().

 {
 #ifndef COPASI_DEBUG
   removeParameter("Random Number Generator");
   removeParameter("Seed");
 #endif
   return true;
 }

bool COptMethodSS::evaluate ( const CVector< C_FLOAT64 > & individual )

private

Evaluate the fitness of one individual

Parameters

const CVector< C_FLOAT64 > & individual

Returns: bool continue

Definition at line 332 of file COptMethodSS.cpp.

References COptProblem::calculate(), COptProblem::checkFunctionalConstraints(), COptProblem::getCalculateValue(), mEvaluationValue, and COptMethod::mpOptProblem.

Referenced by combination(), creation(), and randomize().

 {
   bool Running = true;
 
   // We do not need to check whether the parametric constraints are fulfilled
   // since the parameters are created within the bounds.
 
   // evaluate the fitness
   Running &= mpOptProblem->calculate();
 
   // check whether the functional constraints are fulfilled
   if (!mpOptProblem->checkFunctionalConstraints())
     mEvaluationValue = std::numeric_limits<C_FLOAT64>::infinity();
   else
     mEvaluationValue = mpOptProblem->getCalculateValue();
 
   return Running;
 }

bool COptMethodSS::inforefset	(	C_INT32	type,
		C_INT32	element
	)

private

write informative messages about the refset

Parameters

size_t type

Returns: bool success

bool COptMethodSS::initialize ( )

privatevirtual

Initialize arrays and pointer.

Returns: bool success

Reimplemented from COptMethod.

Definition at line 100 of file COptMethodSS.cpp.

References CProcessReport::addItem(), C_FLOAT64, C_INT32, cleanup(), CRandom::createGenerator(), COptMethod::createMethod(), CProcessReport::finishItem(), CCopasiObject::getObjectParent(), CCopasiParameter::getValue(), CCopasiMethod::HookeJeeves, COptProblem::initialize(), COptMethod::initialize(), COptProblem::initializeSubtaskBeforeOutput(), CCopasiMethod::LevenbergMarquardt, mBestValue, mChild, mChildrenGenerated, mChildVal, mCloseValue, mFreq, mhIterations, mIteration, mIterations, mLocalFreq, CCopasiMethod::mpCallBack, mpLocalMinimizer, mPool, mPoolSize, mPoolVal, COptMethod::mpOptItem, COptMethod::mpOptProblem, mpOptProblemLocal, mPopulationSize, mpRandom, mProb, mRefSet, mRefSetVal, mStuck, mVariableSize, CCopasiParameter::Value::pUINT, CVector< CType >::resize(), COptProblem::setCalculateStatistics(), COptProblem::setCallBack(), CCopasiObject::setObjectParent(), COptMethod::setProblem(), COptProblem::setRandomizeStartValues(), and CCopasiParameterGroup::setValue().

Referenced by optimise().

 {
   cleanup();
 
   size_t i;
 
   if (!COptMethod::initialize())
     {
       if (mpCallBack)
         mpCallBack->finishItem(mhIterations);
 
       return false;
     }
 
   // get total number of iterations
   mIterations = * getValue("Number of Iterations").pUINT;
   // set current iteration to zero
   mIteration = 0;
 
   if (mpCallBack)
     mhIterations =
       mpCallBack->addItem("Current Iteration",
                           mIteration,
                           & mIterations);
 
   mIteration++;
 
 #ifdef COPASI_DEBUG
   mpRandom =
     CRandom::createGenerator(* (CRandom::Type *) getValue("Random Number Generator").pUINT,
                              * getValue("Seed").pUINT);
 #else
   // Use the default random number generator which is the Mersenne Twister
   mpRandom = CRandom::createGenerator();
 #endif
 
   mCloseValue = 0.001;
 
   // frequency of local optimiser; this is hardcoded as 20 per Jose Egea,
   mLocalFreq = 20;
 
   // get number of variables in the problem
   mVariableSize = mpOptItem->size();
 
   // calculate the number of individuals in population
   mPopulationSize = (C_INT32) ceil(1.0 + sqrt(1.0 + 40.0 * mVariableSize) / 2.0);
 
   if (mPopulationSize % 2 != 0) mPopulationSize++;
 
   CFitProblem * pFitProblem = dynamic_cast< CFitProblem * >(mpOptProblem);
 
   if (pFitProblem != NULL)
     {
       // this is a least squares problem (param estimation)
       // let's use our favorite lsq method
       mpLocalMinimizer = COptMethod::createMethod(CCopasiMethod::LevenbergMarquardt);
       // the intermediate local minimizations use a rather relaxed tolerance
       mpLocalMinimizer->setValue("Tolerance", (C_FLOAT64) 1.e-003);
       // TODO: not sure if we should let this one go that long...
       mpLocalMinimizer->setValue("Iteration Limit", (C_INT32) 2000);
     }
   else
     {
       // this is a generic optimisation problem
       // let's use Hooke and Jeeves
       mpLocalMinimizer = COptMethod::createMethod(CCopasiMethod::HookeJeeves);
       // with a rather relaxed tolerance (1e-3) for intermediate minimizations
       mpLocalMinimizer->setValue("Tolerance", (C_FLOAT64) 1.e-003);
       mpLocalMinimizer->setValue("Iteration Limit", (C_INT32) 50);
       mpLocalMinimizer->setValue("Rho", (C_FLOAT64) 0.2);
     }
 
   // local minimization problem (starts as a copy of the current problem)
   if (pFitProblem != NULL)
     {
       // this is a least squares problem (param estimation)
       mpOptProblemLocal = new CFitProblem(*pFitProblem, getObjectParent());
     }
   else
     {
       // this is a generic optimisation problem
       mpOptProblemLocal = new COptProblem(*mpOptProblem, getObjectParent());
     }
 
   // the local optimization method should not have a callback
   mpOptProblemLocal->setCallBack(NULL);
 
   // set object parent (this is needed or else initialize() will fail)
   mpLocalMinimizer->setObjectParent(getObjectParent());
   // we also have to initialize the subtask
   mpOptProblemLocal->initializeSubtaskBeforeOutput();
   // initialize it
   mpOptProblemLocal->initialize();
   // no statistics to be calculated in the local problems
   mpOptProblemLocal->setCalculateStatistics(false);
   // do not randomize the initial values
   mpOptProblemLocal->setRandomizeStartValues(false);
   mpLocalMinimizer->setProblem(mpOptProblemLocal);
 
   // create matrix for the RefSet (population)
   mRefSet.resize(mPopulationSize);
 
   for (i = 0; i < mPopulationSize; i++)
     mRefSet[i] = new CVector< C_FLOAT64 >(mVariableSize);
 
   // create vector for function values (of RefSet members)
   mRefSetVal.resize(mPopulationSize);
   mRefSetVal = std::numeric_limits<C_FLOAT64>::infinity();
 
   // create vector for counting stuck iterations
   mStuck.resize(mPopulationSize);
   mStuck = 0;
 
   // create matrix for the RefSet children
   mChild.resize(mPopulationSize);
 
   for (i = 0; i < mPopulationSize; i++)
     mChild[i] = new CVector< C_FLOAT64 >(mVariableSize);
 
   // create vector for function values (of child members)
   mChildVal.resize(mPopulationSize);
   mChildVal = std::numeric_limits<C_FLOAT64>::infinity();
 
   // we have not generated any children yet
   mChildrenGenerated = false;
 
   // create matrix for the pool of diverse solutions
   // this will also be used to store the initial and
   // final positions of local optimizations
   if (10 * mVariableSize > 2 * mIterations / mLocalFreq)
     {
       mPoolSize = 10 * mVariableSize;
     }
   else
     {
       mPoolSize = 2 * mIterations / mLocalFreq;
     }
 
   mPool.resize(mPoolSize);
 
   for (i = 0; i < mPoolSize; i++)
     mPool[i] = new CVector< C_FLOAT64 >(mVariableSize);
 
   mPoolVal.resize(mPoolSize);
   mPoolVal = std::numeric_limits<C_FLOAT64>::infinity();
 
   // best is infinity, so anything will improve it
   mBestValue = std::numeric_limits<C_FLOAT64>::infinity();
 
   // array for frequencies
   // initialized at 1
   mFreq.resize(mVariableSize);
 
   for (i = 0; i < mVariableSize; i++)
     {
       mFreq[i] = new CVector< C_INT32 >(4);
       *mFreq[i] = 1;
     }
 
   // vector for probabilities
   // initialized at 0
   mProb.resize(4);
   mProb = 0.0;
 
   return true;
 }

void COptMethodSS::initObjects ( )

private

Initialize contained objects.

Definition at line 95 of file COptMethodSS.cpp.

References CCopasiContainer::addObjectReference(), mIteration, and CCopasiObject::ValueInt.

Referenced by COptMethodSS().

 {
   addObjectReference("Current Iteration", mIteration, CCopasiObject::ValueInt);
 }

bool COptMethodSS::localmin	(	CVector< C_FLOAT64 > &	solution,
		C_FLOAT64 &	fval
	)

private

Find a local minimum

Parameters

CVector<	C_FLOAT64 > & solution
C_FLOAT64	& fval

Returns: bool continue

Definition at line 295 of file COptMethodSS.cpp.

References C_INT32, COptProblem::getFunctionEvaluations(), COptProblem::getOptItemList(), COptProblem::getSolutionValue(), COptProblem::getSolutionVariables(), COptProblem::incrementEvaluations(), mpLocalMinimizer, COptMethod::mpOptProblem, mpOptProblemLocal, mVariableSize, COptMethod::optimise(), COptProblem::reset(), and COptProblem::resetEvaluations().

Referenced by childLocalMin(), and optimise().

 {
   bool Running = true;
   unsigned C_INT32 i;
 
   mpOptProblemLocal->reset();
 
   // first we set up the problem
   // (optmethod and optproblem already setup in initialization)
   // let's get the list of parameters
   std::vector<COptItem *> optitem = mpOptProblemLocal->getOptItemList();
 
   // and set them to the values passed in solution
   for (i = 0; i < mVariableSize; i++)
     {
       optitem[i]->setStartValue(solution[i]);
     }
 
   // reset the function counter of the local minimizer
   mpOptProblemLocal->resetEvaluations();
 
   // run it
   Running &= mpLocalMinimizer->optimise();
   // add the function evaluations taken in local to the global problem
   mpOptProblem->incrementEvaluations(mpOptProblemLocal->getFunctionEvaluations());
   // pass the results on to the calling parameters
   fval = mpOptProblemLocal->getSolutionValue();
 
   for (i = 0; i < mVariableSize; i++)
     {
       solution[i] = mpOptProblemLocal->getSolutionVariables()[i];
     }
 
   return Running;
 }

bool COptMethodSS::optimise ( void )

virtual

Execute the optimization algorithm calling simulation routine when needed. It is noted that this procedure can give feedback of its progress by the callback function set with SetCallback. @ return success;

Reimplemented from COptMethod.

Definition at line 1070 of file COptMethodSS.cpp.

References C_FLOAT64, childLocalMin(), cleanup(), closerRefSet(), combination(), creation(), COutputInterface::DURING, CProcessReport::finishItem(), COptProblem::getSolutionVariables(), initialize(), localmin(), mBestValue, mChild, mChildrenGenerated, mChildVal, mCloseValue, mhIterations, mIteration, mIterations, mLocalFreq, mLocalIter, mLocalStored, CCopasiMethod::mpCallBack, mpLocalMinimizer, mPoolSize, COptMethod::mpOptProblem, mPopulationSize, COptMethod::mpParentTask, mRefSet, mRefSetVal, mStuck, CCopasiTask::output(), CProcessReport::progressItem(), randomize(), COptProblem::setSolution(), CCopasiParameterGroup::setValue(), and sortRefSet().

 {
   bool Running = true;
   bool needsort;
   size_t i, j;
   C_FLOAT64 mx, mn, la;
 
   if (!initialize())
     {
       // initialisation failed, we exit
       if (mpCallBack)
         mpCallBack->finishItem(mhIterations);
 
       return false;
     }
 
   mIteration = 0;
 
   // create the Pool of diverse candidate solutions
   Running &= creation();
 
   // best value is (always) at position zero
   // store that value
   mBestValue = mRefSetVal[0];
   // set it upstream
   Running &= mpOptProblem->setSolution(mBestValue, *mRefSet[0]);
   // We found a new best value let's report it.
   mpParentTask->output(COutputInterface::DURING);
 
   // test if the user wants to stop, and do so if needed
   if (!Running)
     {
       if (mpCallBack)
         mpCallBack->finishItem(mhIterations);
 
       cleanup();
       return true;
     }
 
   // mPool is now going to be used to keep track of initial and final
   // points of local minimizations (to avoid running them more than once)
   mPoolSize = 2 * mIterations / mLocalFreq;
   // reset the number of stored minimizations
   mLocalStored = 0;
   // reset the counter for local minimisation
   mLocalIter = 1;
 
   // run the mIterations (and count the creation as being the first)
   for (mIteration = 1; mIteration < mIterations && Running; mIteration++)
     {
       // check for stagnation or similarity
       needsort = false;
 
       for (i = 0; i < mPopulationSize; i++)
         {
           // are we stuck? (20 iterations)
           if (mStuck[i] == 19)
             {
               // substitute this one by a random guess
               Running &= randomize(i);
               needsort = true;
               mStuck[i] = 1;
             }
           else
             {
               // check if another RefSet member is similar to us (relative dist 0.1%)
               for (j = i + 1; j < mPopulationSize; j++)
                 {
                   // is the other one like me?
                   if (closerRefSet(i, j, mCloseValue))
                     //if (distRefSet(i, j) < 0.01)
                     {
                       // randomize the other one because it has worse value
                       Running &= randomize(j);
                       needsort = true;
                       mStuck[j] = 1;
                     }
                 }
             }
         }
 
       // sort the RefSet if needed
       if (needsort) sortRefSet(0, mPopulationSize);
 
       // create children by combination
       Running &= combination();
 
       // check if we have to run a local search
       if (mLocalIter >= mLocalFreq && mChildrenGenerated)
         {
           // carry out a local search
           Running &= childLocalMin();
         }
       else
         {
           // count this
           mLocalIter++;
         }
 
       // substitute the parents for children or increment stuck counter
       needsort = false;
 
       for (i = 0; i < mPopulationSize; i++)
         {
           // check if child was better than parent
           if (mStuck[i] == 0)
             {
               // copy the child into the population
               (*mRefSet[i]) = (*mChild[i]);
               // keep its obj funct value
               mRefSetVal[i] = mChildVal[i];
               // and reset the stuck counter
               mStuck[i] = 1;
               needsort = true;
             }
           else
             {
               // nothing to do, increment parent's stuck counters
               mStuck[i]++;
             }
         }
 
       // sort the RefSet if needed
       if (needsort) sortRefSet(0, mPopulationSize);
 
       // have we made any progress?
       if (mRefSetVal[0] < mBestValue)
         {
           // and store that value
           mBestValue = mRefSetVal[0];
           Running &= mpOptProblem->setSolution(mBestValue, *mRefSet[0]);
           // We found a new best value lets report it.
           mpParentTask->output(COutputInterface::DURING);
         }
 
       if (mpCallBack)
         Running &= mpCallBack->progressItem(mhIterations);
     }
 
   // end of loop for iterations
 
   // the best ever might not be what is on position 0, so bring it back
   *mRefSet[0] = mpOptProblem->getSolutionVariables();
 
   // now let's do a final local minimisation with a tighter tolerance
 
   if (Running)
     {
       mpLocalMinimizer->setValue("Tolerance", (C_FLOAT64) 1.e-006);
       Running &= localmin(*(mRefSet[0]), mRefSetVal[0]);
     }
 
   // has it improved?
   if (mRefSetVal[0] < mBestValue)
     {
       // and store that value
       mBestValue = mRefSetVal[0];
       Running &= mpOptProblem->setSolution(mBestValue, *mRefSet[0]);
       // We found a new best value lets report it.
       mpParentTask->output(COutputInterface::DURING);
     }
 
   if (mpCallBack)
     mpCallBack->finishItem(mhIterations);
 
   cleanup();
   return true;
 }

bool COptMethodSS::randomize ( C_INT32 i )

private

randomize an element of the RefSet

Definition at line 352 of file COptMethodSS.cpp.

References C_FLOAT64, C_INT32, COptItem::checkConstraint(), evaluate(), COptItem::getLowerBoundValue(), CRandom::getRandomCC(), COptItem::getUpperBoundValue(), max, mEvaluationValue, min, mpRandom, mRefSet, mRefSetVal, mStuck, and mVariableSize.

Referenced by optimise().

 {
   C_FLOAT64 mn, mx, la; // for boundaries of rnd
   bool Running = true;  // flag for invalid values
 
   for (C_INT32 j = 0; j < mVariableSize; j++)
     {
       // get pointers to appropriate elements (easier reading of code)
       COptItem & OptItem = *(*mpOptItem)[j];
       C_FLOAT64 & Sol = (*mRefSet[i])[j];
       // calculate lower and upper bounds for this variable
       mn = *OptItem.getLowerBoundValue();
       mx = *OptItem.getUpperBoundValue();
 
       try
         {
           // calculate orders of magnitude of the interval
           la = log10(mx) - log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()));
 
           // determine if linear or log scale
           if ((mn < 0.0) || (mx <= 0.0))
             Sol = mn + mpRandom->getRandomCC() * (mx - mn);
           else
             {
               if (la < 1.8)
                 Sol = mn + mpRandom->getRandomCC() * (mx - mn);
               else
                 Sol = pow(10.0, log10(std::max(mn, std::numeric_limits< C_FLOAT64 >::min()))
                           + la * mpRandom->getRandomCC());
             }
         }
       catch (...)
         {
           // if there were errors, let's just stay with the midpoint
           Sol = (mx + mn) * 0.5;
         }
 
       // force it to be within the bounds
       switch (OptItem.checkConstraint(Sol))
         {
           case - 1:
             Sol = *OptItem.getLowerBoundValue();
             break;
 
           case 1:
             Sol = *OptItem.getUpperBoundValue();
             break;
         }
 
       // We need to set the value here so that further checks take
       // account of the value.
       (*(*mpSetCalculateVariable)[j])(Sol);
     }
 
   // calculate its fitness
   Running = evaluate(*mRefSet[i]);
   mRefSetVal[i] = mEvaluationValue;
   // reset the stuck flag
   mStuck[i] = 1;
   return Running;
 }

bool COptMethodSS::serializechildren	(	size_t	first,
		size_t	last
	)

private

Serialize the children for debug

Parameters

size_t	first
size_t	last (default: population size)

Returns: bool success

bool COptMethodSS::serializepool	(	size_t	first,
		size_t	last
	)

private

Serialize the pool for debug

Parameters

size_t	first
size_t	last (default: pool size)

Returns: bool success

bool COptMethodSS::serializeprob ( void )

private

Serialize the probability table for debug

Returns: bool success

bool COptMethodSS::serializerefset	(	C_INT32	first,
		C_INT32	last
	)

private

Serialize the refset for debug

Parameters

size_t	first
size_t	last (default: population size)

Returns: bool success

bool COptMethodSS::serializevector	(	CVector< C_FLOAT64 >	x,
		C_FLOAT64	xval
	)

private

Serialize a vector to the refset file ,for debug

Parameters

x	the vector
xval	function value for that vector

Returns: bool success

void COptMethodSS::sortRefSet	(	C_INT32	lower,
		C_INT32	upper
	)

private

sort the RefSet array between position lower and upper

Definition at line 729 of file COptMethodSS.cpp.

References C_FLOAT64, C_INT32, mRefSet, mRefSetVal, and mStuck.

Referenced by creation(), and optimise().

 {
   C_INT32 i, j, k;
   C_INT32 parent, child;
   CVector< C_FLOAT64 > *tempvec;
   C_FLOAT64 tempval;
 
   // Use heap sort
   for (i = lower + 1; i < upper; i++)
     {
       // bubble-up
       child = i;
 
       for (;;)
         {
           if (child == 0) break;
 
           parent = floor((double)(child - 1) / 2);
 
           if (mRefSetVal[child] < mRefSetVal[parent])
             {
               // swap with parent
               tempval = mRefSetVal[child];
               mRefSetVal[child] = mRefSetVal[parent];
               mRefSetVal[parent] = tempval;
               tempval = mStuck[child];
               mStuck[child] = mStuck[parent];
               mStuck[parent] = tempval;
               tempvec = mRefSet[child];
               mRefSet[child] = mRefSet[parent];
               mRefSet[parent] = tempvec;
               // make parent the new child
               child = parent;
             }
           else break;
         }
     }
 
   // since some leafs are not in order in the array, we now do
   // a bubble sort (note: this is best case for bubble sort)
   j = upper - 1;
 
   do
     {
       k = lower;
 
       for (i = lower; i < j; i++)
         {
           if (mRefSetVal[i] > mRefSetVal[i + 1])
             {
               tempval = mRefSetVal[i];
               mRefSetVal[i] = mRefSetVal[i + 1];
               mRefSetVal[i + 1] = tempval;
               tempval = mStuck[i];
               mStuck[i] = mStuck[i + 1];
               mStuck[i + 1] = tempval;
               tempvec = mRefSet[i];
               mRefSet[i] = mRefSet[i + 1];
               mRefSet[i + 1] = tempvec;
               k = i;
             }
         }
 
       j = k;
     }
   while (j > lower);
 }