It works well enough to solve polynomial maxima.

14 files changed, 494 insertions, 2 deletions

diff --git a/src/config.h b/src/config.h
new file mode 100644
index 0000000..2563eb9
--- /dev/null
+++ b/src/config.h
@@ -0,0 +1,12 @@
+#ifndef GENETIC_CONFIG_H
+#define GENETIC_CONFIG_H
+
+#include <stdint.h>
+
+namespace Genetic
+{
+        typedef uint32_t Time;
+        typedef uint32_t PhenotypeId;
+};
+
+#endif
diff --git a/src/explicitsimulation.cpp b/src/explicitsimulation.cpp
new file mode 100644
index 0000000..281e7e5
--- /dev/null
+++ b/src/explicitsimulation.cpp
@@ -0,0 +1,139 @@
+#include "genetic/explicitsimulation.h"
+#include "genetic/operator.h"
+#include "genetic/fitnessfunction.h"
+
+#include <bu/random.h>
+#include <bu/sio.h>
+
+using namespace Bu;
+
+Genetic::ExplicitSimulation::ExplicitSimulation( Genetic::Operator *pOper,
+                Genetic::FitnessFunction *pFunc, int iPopSize, float fKeep,
+                float fRandom, bool bKeepBest ) :
+        pOper( pOper ),
+        pFunc( pFunc ),
+        iPopSize( iPopSize ),
+        fKeep( fKeep ),
+        fRandom( fRandom ),
+        bKeepBest( bKeepBest )
+{
+        for( int j = 0; j < iPopSize; j++ )
+        {
+                xPop.addPhenotype( pOper->random() );
+        }
+
+        updateFitness();
+}
+
+Genetic::ExplicitSimulation::~ExplicitSimulation()
+{
+        delete pOper;
+        delete pFunc;
+}
+
+void Genetic::ExplicitSimulation::timestep()
+{
+        PhenotypeList lNew;
+        
+        int iChildren = iPopSize*(1.0-fKeep-fRandom);
+
+        // Create children
+        for( int j = 0; j < iChildren; j++ )
+        {
+                PhenotypeList lParents;
+                for( int k = 0; k < pOper->parentCount(); k++ )
+                        lParents.append( xPop.getPhenotype( selectWeighted() ) );
+                lNew.append( pOper->mate( lParents ) );
+        }
+
+        // Select phenotypes for keeping
+        int iKeep = iPopSize*fKeep;
+        FitnessHash hTempFitness;
+        for( int j = 0; j < iKeep; j++ )
+        {
+                Genetic::PhenotypeId id = selectWeighted();
+                lNew.append( xPop.takePhenotype( id ) );
+                hTempFitness.insert( id, hFitness.get( id ) );
+                hFitness.erase( id );
+                dTotalFitness -= hTempFitness.get( id );
+        }
+
+        if( bKeepBest && hFitness.has( uMaxFitness ) )
+        {
+                lNew.append( xPop.takePhenotype( uMaxFitness ) );
+                hTempFitness.insert( uMaxFitness, hFitness.get( uMaxFitness ) );
+                hFitness.erase( uMaxFitness );
+                dTotalFitness -= hTempFitness.get( uMaxFitness );
+        }
+
+        // Fill in the remainder with random phenotypes
+        while( lNew.getSize() < iPopSize )
+        {
+                lNew.append( pOper->random() );
+        }
+
+        // Refill the population
+        hFitness = hTempFitness;
+        xPop.clear();
+        xPop.timestep();
+        for( PhenotypeList::iterator i = lNew.begin(); i; i++ )
+        {
+                xPop.addPhenotype( *i );
+        }
+
+        updateFitness();
+}
+
+Genetic::PhenotypeId Genetic::ExplicitSimulation::selectWeighted()
+{
+        double dSel = Bu::Random::randNorm()*dTotalFitness;
+        double dRun = 0.0;
+
+        for( FitnessHash::iterator i = hFitness.begin(); i; i++ )
+        {
+                dRun += *i;
+                if( dSel < dRun )
+                        return i.getKey();
+        }
+
+        sio << "Genetic::ExplicitSimulation::selectWeighted() - failed, picked max"
+                << sio.nl;
+
+        return uMaxFitness;
+}
+
+void Genetic::ExplicitSimulation::updateFitness()
+{
+        dMinFitness = -1.0;
+        dTotalFitness = 0.0;
+        for( Population::iterator i = xPop.begin(); i; i++ )
+        {
+                double dFitness;
+                if( hFitness.has( i.getKey() ) )
+                {
+                        dFitness = hFitness.get( i.getKey() );
+                }
+                else
+                {
+                        dFitness = (*pFunc)( *i );
+                        if( dFitness < 0.0 )
+                                dFitness = 0.0;
+                        hFitness.insert( i.getKey(), dFitness );
+                }
+                dTotalFitness += dFitness;
+                if( dMinFitness < 0.0 )
+                {
+                        dMinFitness = dMaxFitness = dFitness;
+                        uMaxFitness = i.getKey();
+                }
+                else if( dMinFitness > dFitness )
+                {
+                        dMinFitness = dFitness;
+                }
+                else if( dMaxFitness < dFitness )
+                {
+                        dMaxFitness = dFitness;
+                        uMaxFitness = i.getKey();
+                }
+        }
+}
diff --git a/src/explicitsimulation.h b/src/explicitsimulation.h
new file mode 100644
index 0000000..ebddd62
--- /dev/null
+++ b/src/explicitsimulation.h
@@ -0,0 +1,48 @@
+#ifndef GENETIC_EXPLICIT_SIMULATION_H
+#define GENETIC_EXPLICIT_SIMULATION_H
+
+#include "genetic/population.h"
+#include "genetic/config.h"
+
+namespace Genetic
+{
+        class Operator;
+        class FitnessFunction;
+
+        class ExplicitSimulation
+        {
+        public:
+                ExplicitSimulation( Operator *pOper, FitnessFunction *pFunc,
+                                int iPopSize, float fKeep, float fRandom, bool bKeepBest=true );
+                virtual ~ExplicitSimulation();
+
+                void timestep();
+
+                Genetic::PhenotypeId selectWeighted();
+
+                double getMinFitness() const { return dMinFitness; }
+                double getMaxFitness() const { return dMaxFitness; }
+
+        protected:
+                void updateFitness();
+
+        protected:
+                Population xPop;
+                Operator *pOper;
+                FitnessFunction *pFunc;
+
+        private:
+                int iPopSize;
+                float fKeep;
+                float fRandom;
+                bool bKeepBest;
+                typedef Bu::Hash<Genetic::PhenotypeId, double> FitnessHash;
+                FitnessHash hFitness;
+                double dMinFitness;
+                double dMaxFitness;
+                double dTotalFitness;
+                PhenotypeId uMaxFitness;
+        };
+};
+
+#endif
diff --git a/src/fitnessfunction.cpp b/src/fitnessfunction.cpp
new file mode 100644
index 0000000..dd2ee08
--- /dev/null
+++ b/src/fitnessfunction.cpp
@@ -0,0 +1,10 @@
+#include "genetic/fitnessfunction.h"
+
+Genetic::FitnessFunction::FitnessFunction()
+{
+}
+
+Genetic::FitnessFunction::~FitnessFunction()
+{
+}
+
diff --git a/src/fitnessfunction.h b/src/fitnessfunction.h
new file mode 100644
index 0000000..c41f733
--- /dev/null
+++ b/src/fitnessfunction.h
@@ -0,0 +1,18 @@
+#ifndef GENETIC_FITNESS_FUNCTION_H
+#define GENETIC_FITNESS_FUNCTION_H
+
+namespace Genetic
+{
+        class Phenotype;
+
+        class FitnessFunction
+        {
+        public:
+                FitnessFunction();
+                virtual ~FitnessFunction();
+
+                virtual double operator()( Phenotype *pTest )=0;
+        };
+};
+
+#endif
diff --git a/src/phenotype.cpp b/src/phenotype.cpp
index 1eca170..b965ffd 100644
--- a/src/phenotype.cpp
+++ b/src/phenotype.cpp
@@ -1,6 +1,9 @@
 #include "genetic/phenotype.h"
 
-Genetic::Phenotype::Phenotype()
+Genetic::Phenotype::Phenotype() :
+        tCreated( 0 ),
+        uId( 0 ),
+        bActive( false )
 {
 }
 
@@ -8,3 +11,19 @@ Genetic::Phenotype::~Phenotype()
 {
 }
 
+bool Genetic::Phenotype::hasProperty( const Bu::String &sKey ) const
+{
+        return hProp.has( sKey );
+}
+
+Bu::Variant Genetic::Phenotype::getProperty( const Bu::String &sKey ) const
+{
+        return hProp.get( sKey );
+}
+
+void Genetic::Phenotype::setProperty( const Bu::String &sKey,
+                Bu::Variant vValue )
+{
+        hProp.insert( sKey, vValue );
+}
+
diff --git a/src/phenotype.h b/src/phenotype.h
index 7c564d5..0eafb78 100644
--- a/src/phenotype.h
+++ b/src/phenotype.h
@@ -2,6 +2,10 @@
 #define GENETIC_PHENOTYPE_H
 
 #include <bu/string.h>
+#include <bu/hash.h>
+#include <bu/variant.h>
+
+#include "genetic/config.h"
 
 namespace Genetic
 {
@@ -29,7 +33,22 @@ namespace Genetic
 
                 virtual Bu::String toString()=0;
 
+                bool hasProperty( const Bu::String &sKey ) const;
+                Bu::Variant getProperty( const Bu::String &sKey ) const;
+                void setProperty( const Bu::String &sKey, Bu::Variant vValue );
+                void setCreated( Genetic::Time t, Genetic::PhenotypeId id )
+                        { tCreated = t; uId = id; bActive = true; }
+                Genetic::Time getCreated() const { return tCreated; }
+                Genetic::PhenotypeId getId() const { return uId; }
+
+                bool isActive() const { return bActive; }
+
         private:
+                Genetic::Time tCreated;
+                Genetic::PhenotypeId uId;
+                bool bActive;
+                typedef Bu::Hash<Bu::String, Bu::Variant> PropHash;
+                PropHash hProp;
         };
 };
 
diff --git a/src/phenotypebinary.cpp b/src/phenotypebinary.cpp
index d16588b..68f9ebf 100644
--- a/src/phenotypebinary.cpp
+++ b/src/phenotypebinary.cpp
@@ -99,3 +99,37 @@ Bu::String Genetic::PhenotypeBinary::toString()
         return sRet;
 }
 
+void Genetic::PhenotypeBinary::extractBits( uint32_t &rTarget, int iStart,
+                int iBits )
+{
+        rTarget = 0;
+        if( iBits > sizeof(rTarget)*8 )
+                iBits = sizeof(rTarget)*8;
+
+        // This is pretty much the same problem as copyFrom, so I'm doing it the
+        // same way for now.
+        for( int j = 0; j < iBits; j++ )
+        {
+                int ws = wordWithBit(j+iStart);
+                if( (aGenes[ws]&(1<<((j+iStart)%BPW))) != 0 )
+                        rTarget |= (1<<(j%BPW));
+        }
+}
+
+void Genetic::PhenotypeBinary::extractBits( uint64_t &rTarget, int iStart,
+                int iBits )
+{
+        rTarget = 0;
+        if( iBits > sizeof(rTarget)*8 )
+                iBits = sizeof(rTarget)*8;
+
+        // This is pretty much the same problem as copyFrom, so I'm doing it the
+        // same way for now.
+        for( int j = 0; j < iBits; j++ )
+        {
+                int ws = wordWithBit(j+iStart);
+                if( (aGenes[ws]&(1<<((j+iStart)%BPW))) != 0 )
+                        rTarget |= (1ll<<(j));
+        }
+}
+
diff --git a/src/phenotypebinary.h b/src/phenotypebinary.h
index 72435ed..cbae392 100644
--- a/src/phenotypebinary.h
+++ b/src/phenotypebinary.h
@@ -19,6 +19,9 @@ namespace Genetic
                                 int iCount, int iDest=-1 );
                 virtual Bu::String toString();
 
+                void extractBits( uint32_t &rTarget, int iStart, int iBits );
+                void extractBits( uint64_t &rTarget, int iStart, int iBits );
+
         private:
                 int iSize;
                 int iWords;
diff --git a/src/population.cpp b/src/population.cpp
index 85b859e..fbd6d6c 100644
--- a/src/population.cpp
+++ b/src/population.cpp
@@ -1,10 +1,85 @@
 #include "genetic/population.h"
+#include "genetic/phenotype.h"
 
-Genetic::Population::Population()
+#include <bu/mutexlocker.h>
+
+Genetic::Population::Population() :
+        uNextId( 0 ),
+        uTime( 0 )
 {
 }
 
 Genetic::Population::~Population()
 {
+        for( PhenotypeHash::iterator i = hPhenotype.begin(); i; i++ )
+                delete *i;
+}
+
+Genetic::PhenotypeId Genetic::Population::addPhenotype(
+                Genetic::Phenotype *pNew )
+{
+        Bu::MutexLocker ml( mGeneral );
+        if( !pNew->isActive() )
+                pNew->setCreated( uTime, uNextId++ );
+
+        hPhenotype.insert( pNew->getId(), pNew );
+        return pNew->getId();
+}
+
+bool Genetic::Population::hasProperty( Genetic::PhenotypeId id,
+                const Bu::String &sKey ) const
+{
+        Bu::MutexLocker ml( mGeneral );
+        return hPhenotype.get( id )->hasProperty( sKey );
+}
+
+Bu::Variant Genetic::Population::getProperty(
+                Genetic::PhenotypeId id, const Bu::String &sKey ) const
+{
+        Bu::MutexLocker ml( mGeneral );
+        return hPhenotype.get( id )->getProperty( sKey );
+}
+
+void Genetic::Population::setProperty( Genetic::PhenotypeId id,
+                const Bu::String &sKey, Bu::Variant vValue )
+{
+        Bu::MutexLocker ml( mGeneral );
+
+        hPhenotype.get( id )->setProperty( sKey, vValue );
+}
+
+void Genetic::Population::delPhenotype( PhenotypeId id )
+{
+        Bu::MutexLocker ml( mGeneral );
+        delete hPhenotype.get( id );
+        hPhenotype.erase( id );
+}
+
+Genetic::Phenotype *Genetic::Population::takePhenotype( PhenotypeId id )
+{
+        Bu::MutexLocker ml( mGeneral );
+        Phenotype *pRet = hPhenotype.get( id );
+        hPhenotype.erase( id );
+        return pRet;
+}
+
+void Genetic::Population::clear()
+{
+        Bu::MutexLocker ml( mGeneral );
+        for( PhenotypeHash::iterator i = hPhenotype.begin(); i; i++ )
+                delete *i;
+        hPhenotype.clear();
+}
+
+Genetic::PhenotypeId Genetic::Population::timestep()
+{
+        Bu::MutexLocker ml( mGeneral );
+        return (++uTime);
+}
+
+Genetic::PhenotypeId Genetic::Population::time() const
+{
+        Bu::MutexLocker ml( mGeneral );
+        return uTime;
 }
 
diff --git a/src/population.h b/src/population.h
index 0fdad32..c2d1bec 100644
--- a/src/population.h
+++ b/src/population.h
@@ -1,13 +1,49 @@
 #ifndef GENETIC_POPULATION_H
 #define GENETIC_POPULATION_H
 
+#include <bu/string.h>
+#include <bu/hash.h>
+#include <bu/mutex.h>
+
+#include "genetic/config.h"
+
 namespace Genetic
 {
+        class Phenotype;
+
         class Population
         {
         public:
                 Population();
                 virtual ~Population();
+
+                typedef Bu::Hash<PhenotypeId, Phenotype *> PhenotypeHash;
+                typedef PhenotypeHash::iterator iterator;
+                typedef PhenotypeHash::const_iterator const_iterator;
+
+                iterator begin() { return hPhenotype.begin(); }
+                const_iterator begin() const { return hPhenotype.begin(); }
+                int getSize() const { return hPhenotype.getSize(); }
+
+                PhenotypeId addPhenotype( Phenotype *pNew );
+                Phenotype *getPhenotype( PhenotypeId id )
+                        { return hPhenotype.get( id ); }
+                bool hasProperty( PhenotypeId id, const Bu::String &sKey ) const;
+                Bu::Variant getProperty( PhenotypeId id, const Bu::String &sKey ) const;
+                void setProperty( PhenotypeId id, const Bu::String &sKey,
+                                Bu::Variant vValue );
+                void delPhenotype( PhenotypeId id );
+                Phenotype *takePhenotype( PhenotypeId id );
+                void clear();
+
+                PhenotypeId timestep();
+                PhenotypeId time() const;
+
+        private:
+                PhenotypeId uNextId;
+                PhenotypeId uTime;
+                PhenotypeHash hPhenotype;
+                mutable Bu::Mutex mGeneral;
         };
 };
 
diff --git a/src/tests/maxima/fitnessfunctioneq.cpp b/src/tests/maxima/fitnessfunctioneq.cpp
new file mode 100644
index 0000000..5694507
--- /dev/null
+++ b/src/tests/maxima/fitnessfunctioneq.cpp
@@ -0,0 +1,25 @@
+#include "fitnessfunctioneq.h"
+#include "genetic/phenotypebinary.h"
+
+FitnessFunctionEq::FitnessFunctionEq()
+{
+}
+
+FitnessFunctionEq::~FitnessFunctionEq()
+{
+}
+
+double FitnessFunctionEq::operator()( Genetic::Phenotype *pTest )
+{
+        Genetic::PhenotypeBinary *pbTest =
+                dynamic_cast<Genetic::PhenotypeBinary *>( pTest );
+        if( pbTest == NULL )
+                return 0.0;
+
+        uint32_t uRaw;
+        pbTest->extractBits( uRaw, 0, 32 );
+        double x = ((double)uRaw / (double)(0xfffffffful))*5.0 - 2.5;
+
+        return -1.8*(x*x*x*x) + 0.86*(x*x*x) + 4.0*(x*x);
+}
+
diff --git a/src/tests/maxima/fitnessfunctioneq.h b/src/tests/maxima/fitnessfunctioneq.h
new file mode 100644
index 0000000..7139532
--- /dev/null
+++ b/src/tests/maxima/fitnessfunctioneq.h
@@ -0,0 +1,15 @@
+#ifndef GENETIC_FITNESS_FUNCTION_EQ_H
+#define GENETIC_FITNESS_FUNCTION_EQ_H
+
+#include "genetic/fitnessfunction.h"
+
+class FitnessFunctionEq : public Genetic::FitnessFunction
+{
+public:
+        FitnessFunctionEq();
+        virtual ~FitnessFunctionEq();
+
+        virtual double operator()( Genetic::Phenotype *pTest );
+};
+
+#endif
diff --git a/src/tests/maxima/main.cpp b/src/tests/maxima/main.cpp
new file mode 100644
index 0000000..ab90e0b
--- /dev/null
+++ b/src/tests/maxima/main.cpp
@@ -0,0 +1,39 @@
+#include "genetic/population.h"
+#include "genetic/phenotypebinary.h"
+#include "genetic/operatorbasic.h"
+#include "genetic/explicitsimulation.h"
+#include "fitnessfunctioneq.h"
+
+#include <time.h>
+
+#include <bu/random.h>
+#include <bu/sio.h>
+using namespace Bu;
+
+int main( int argc, char *argv[] )
+{
+        Bu::Random::seed( time( NULL ) );
+        sio << "Global maxima equation test" << sio.nl
+                << "  - -1.8*x^4 + 0.86*x^3 + 4.0*x^2 == 3.53518 (approx)" << sio.nl
+                << sio.nl;
+
+        Genetic::ExplicitSimulation ex(
+                new Genetic::OperatorBasic(
+                        new Genetic::PhenotypeBinary( 32 ),
+                        0.01
+                        ),
+                new FitnessFunctionEq(),
+                100,
+                .1, .05
+                );
+
+        for( int j = 0; j < 100; j++ )
+        {
+                ex.timestep();
+                sio << ex.getMinFitness() << " - " << ex.getMaxFitness() <<
+                        sio.nl;
+        }
+
+        return 0;
+}
+