DNASTORE.H:
#ifndef DNASTORE_H
#define DNASTORE_H
#include <iostream>
#include <fstream>
using namespace std;
class DNAStore
{
public:
DNAStore(int aTraits, int aMutations, int aGoodMutations,
int aBadMutations, double aInc, double aFitness);
// For file load operations
DNAStore();
// Reset the store's contents
void reset(int aTraits, int aMutations, int aGoodMutations,
int aBadMutations, double aInc, double aFitness);
friend ostream& operator<<( ostream& os, DNAStore& dna );
friend istream& operator>>( istream& os, DNAStore& dna );
// Compilation constants
enum {
// Max number of selectable traits to have
TRAITS = 23,
// How many copies of gene in the DNA
// Note that you can't just change COPIES and expect
// everything to work!
COPIES = 2,
// How many children per critter
// Note that if you change CHILDREN, you need to fix up
// the ages that Critters reproduce at.
CHILDREN = 2,
// Base count of critters in a generation
CRITTERS = 1000000,
// How old a critter has to be before reproducing
MATURE = 17,
// Years between children
CHILD_GAP = 6,
};
void summary();
//Run a trial for a setup
void runTrial();
// Print some mutation info
void printStats();
protected:
// Constants derived from the public ones; leave these alone.
enum {
// Variance in reproduction age
CHILD_VAR = CHILD_GAP / 2,
// Average length of a generation
GEN = MATURE + CHILD_GAP * (CHILDREN - 1) / 2,
// Maximum age of a critter
MAX_AGE = MATURE + CHILD_GAP * (CHILDREN - 1),
// Age at which a critter dies
DEAD = MAX_AGE + 1,
// Size of the critter array; needs to hold everyone at once
MAX_CRIT = CRITTERS * 3,
// Normal beneficial mutation level
BENEFICIAL = 1,
};
// Reproduce a kid from two parents
void reproduce(int aKid, int aP1, int aP2);
// Net beneficial mutations for this critter.
int fitness(int aPos);
// Advance a year with all the reproduction involved
void nextYear();
// Print stats and see if we should stop yet
bool analyze();
// Check dates and times; pause if necessary
bool pauseExecution();
// Pick a number between 0 and max
int choice(int aMax);
// Generate a random number between 0 and 1
double survivalChance();
// Rate of increase, must be <= CHILDREN
double mInc;
// Fitness improvement per beneficial mutation
double mFitness;
// Number of bad mutations per good one
int mBadMutations;
// Number of good mutations per generation
int mGoodMutations;
// How many years have gone by
int mYears;
// How many traits to actually use
int mTraits;
// How many mutations have been substituted
int mMutations;
// How many mutations need to be substituted
int mMaxMutations;
// What day was the last save done?
int mLastSave;
// The array of critter DNA
signed char mDNA[MAX_CRIT][COPIES][TRAITS];
// Which critter is the first one younger than this;
// Note that the youngest are at the end, but the storage
// is circular.
int mAges[DEAD + 1];
// Debugging variables
double mFitSum;
double mFitCount;
};
ostream& operator<<( ostream& os, DNAStore& dna );
istream& operator>>( istream& os, DNAStore& dna );
#endif
DNASTORE.CPP
#include <stdlib.h>
#include <iostream>
#include "dnastore.h"
#include "math.h"
#include "time.h"
#include "unistd.h"
#include "assert.h"
#define MAX(a, b) (a <= b ? b : a)
DNAStore::DNAStore(int aTraits, int aMutations, int aGoodMutations,
int aBadMutations, double aInc, double aFitness) :
mTraits(aTraits), mMaxMutations(aMutations),
mGoodMutations(aGoodMutations), mBadMutations(aBadMutations),
mInc(aInc), mFitness(aFitness), mLastSave(-1)
{
reset(aTraits, aMutations, aGoodMutations, aBadMutations,
aInc, aFitness);
}
// Reset the store's contents
void DNAStore::reset(int aTraits, int aMutations, int aGoodMutations,
int aBadMutations, double aInc, double aFitness)
{
mYears = 0;
mMutations = 0;
mTraits = aTraits;
mMaxMutations = aMutations;
mBadMutations = aBadMutations;
mInc = aInc;
mFitness = aFitness;
mLastSave = -1;
mGoodMutations = aGoodMutations;
// Initialize the array
for (int i = 0; i < MAX_CRIT; i++)
{
for (int j = 0; j < COPIES; j++)
{
for (int k = 0; k < TRAITS; k++)
{
mDNA[i][j][k] = 0;
}
}
}
// Start off a number of critters
for (int i = 0; i < DEAD + 1; i++)
{
mAges[i] = CRITTERS - ((i * CRITTERS) / DEAD);
}
// One beneficial mutation of each type to start
for (int i = 0; i < mTraits; i++)
{
mDNA[CRITTERS - 1 - i][0][i] = 1;
}
mFitSum = 0;
mFitCount = 0;
}
DNAStore::DNAStore() : mLastSave(-1)
{
// Expect a file read to fill all this in.
}
ostream& operator<<( ostream& os, DNAStore& dna )
{
os << dna.mInc << endl;
os << dna.mFitness << endl;
os << dna.mGoodMutations << endl;
os << dna.mBadMutations << endl;
os << dna.mYears << endl;
os << dna.mTraits << endl;
os << dna.mMutations << endl;
os << dna.mMaxMutations << endl;
for (int i = 0; i < dna.DEAD + 1; i++)
{
os << dna.mAges[i] << " ";
}
os << endl;
for (int i = 0; i < dna.MAX_CRIT; i++)
{
for (int j = 0; j < dna.COPIES; j++)
{
for (int k = 0; k < dna.TRAITS; k++)
{
os.put(dna.mDNA[i][j][k]);
}
}
os << endl;
}
os << endl;
return os;
}
istream& operator>>( istream& is, DNAStore& dna )
{
is >> dna.mInc;
is >> dna.mFitness;
is >> dna.mGoodMutations;
is >> dna.mBadMutations;
is >> dna.mYears;
is >> dna.mTraits;
is >> dna.mMutations;
is >> dna.mMaxMutations;
for (int i = 0; i < dna.DEAD + 1; i++)
{
is >> dna.mAges[i];
}
is.get(); // space
is.get(); // newline
for (int i = 0; i < dna.MAX_CRIT; i++)
{
for (int j = 0; j < dna.COPIES; j++)
{
for (int k = 0; k < dna.TRAITS; k++)
{
dna.mDNA[i][j][k] = is.get();
}
}
is.get(); // newline
}
return is;
}
// Reproduce a kid from two parents
void DNAStore::reproduce(int aKid, int aP1, int aP2)
{
int lBits = 0;
int c1 = 0;
int c2 = 0;
for (int i = 0; i < TRAITS; i++)
{
if (lBits <= 0)
{
c1 = rand() % 0x1000;
c2 = rand() % 0x1000;
lBits = 12;
}
// Choose one chromosome randomly from each parent.
mDNA[aKid][0][i] = mDNA[aP1][c1%COPIES][i];
mDNA[aKid][1][i] = mDNA[aP2][c2%COPIES][i];
c1 >>= 1;
c2 >>= 1;
lBits--;
}
}
// Net beneficial mutations for this critter.
int DNAStore::fitness(int aPos)
{
int retval = 0;
for (int i = 0; i < mTraits; i++)
{
retval += MAX(mDNA[aPos][0][i], mDNA[aPos][1][i]);
}
return retval;
}
// Pick a number between 0 and max
int DNAStore::choice(int aMax)
{
// Build a random number that's much bigger than aMax
// NB: Assumes RAND_MAX is 0x7fff as per g++
int myrand = rand();
int lNeed = aMax << 4;
while (lNeed > RAND_MAX)
{
lNeed >>= 15;
myrand = myrand << 15 + rand();
}
return myrand % aMax;
}
// Generate a random number between 0 and 1
double DNAStore::survivalChance()
{
double retval = (((rand() % 0x1000) << 12) + rand() % 0x1000);
retval = retval / 0x1000000;
mFitSum += retval;
mFitCount++;
return retval;
}
// Advance a year with all the reproduction involved
void DNAStore::nextYear()
{
// What fraction over or under the standard population are we?
int pop = mAges[0] - mAges[DEAD];
double lCrowding = (pop - CRITTERS) * mTraits;
lCrowding /= CRITTERS;
// Count all the kids that are going to be born this year
int lKids = mAges[MATURE] - mAges[MATURE + CHILD_VAR];
lKids += mAges[DEAD - CHILD_VAR] - mAges[DEAD];
lKids /= CHILD_VAR;
// Compare to average for this population; get a second crowding
// factor that will even out yearly reproduction
double lKidCrowd = lKids * DEAD / 2;
lKidCrowd /= pop;
lKidCrowd -= 1;
// Range of adults (possible partners for reproducers)
int lAdRange = mAges[MATURE] - mAges[DEAD];
// Where to start putting kids in the array
int lKid = mAges[0];
// Debugging stats
int lKidCountGood = 0;
int lKidCountBad = 0;
double lKidFitGood = 0;
double lKidFitBad = 0;
for (int age = MATURE; age < DEAD; age++)
{
int lStart = mAges[age + 1];
int lFinish = mAges[age];
int lNum = 0;
for (int i = lStart; i < lFinish; i++)
{
int aPos = i % MAX_CRIT;
int lOffset = aPos % (CHILD_VAR);
// The offset is to spread out each generation's reproduction
// as per real life and to smooth out the population.
// It is controlled by array position - cheap to calculate.
if ((age == MATURE + lOffset) ||
(age == MATURE + CHILD_GAP - lOffset))
{
lNum++;
int i2 = choice(lAdRange);
reproduce(lKid % MAX_CRIT, i % MAX_CRIT, i2 % MAX_CRIT);
double lFit = fitness(lKid % MAX_CRIT) - lCrowding - lKidCrowd;
double lFitChance = pow(mFitness, lFit);
// No exceeding Haldane limit!
if (lFitChance > mInc)
{
lFitChance = mInc;
}
// Chance this kid survives:
lFitChance /= CHILDREN;
if (lFitChance >= survivalChance())
{
lKid++;
lKidFitGood += lFitChance;
lKidCountGood++;
}
else
{
lKidFitBad += lFitChance;
lKidCountBad++;
}
}
}
}
// We keep increasing our indices and getting back inside
// the array with mod operations; but there's a limit to that.
int lShift = 0;
if (mAges[DEAD] > MAX_CRIT)
{
lShift = MAX_CRIT;
}
// Crowding should stop the population getting this big.
if (mAges[0] - mAges[DEAD] > MAX_CRIT)
{
cout << "Population exceeded limits!!!" << endl;
for (int i = 0; i < DEAD + 1; i++)
{
cout << "Array " << i << " is " << mAges[i] << endl;
}
cout << "Fitsum: " << mFitSum << " Count: " << mFitCount <<
" Avg: " << (mFitSum / mFitCount) << endl;
assert(mAges[0] - mAges[DEAD] <= MAX_CRIT);
}
// Everyone shifts up an age, and the oldest fall off the end.
for (int i = DEAD; i > 0; i--)
{
mAges[i] = mAges[i-1] - lShift;
}
mAges[0] = lKid - lShift;
mYears++;
}
// Print stats and see if we should stop yet
bool DNAStore::analyze()
{
bool bStop = false;
cout << ".";
int maxval[TRAITS];
for (int i = 0; i < TRAITS; i++)
{
maxval[i] = 0;
}
int pop = mAges[0] - mAges[DEAD];
int lBase = mAges[DEAD];
for (int t = 0; t < mTraits; t++)
{
int tc = 0;
for (int i = 0; i < pop; i++)
{
int cr = (i + lBase) % MAX_CRIT;
int loc = MAX(mDNA[cr][0][t], mDNA[cr][1][t]);
tc += mDNA[cr][0][t] + mDNA[cr][1][t];
maxval[t] = MAX(maxval[t], loc);
}
// Avg critter has a beneficial mutation in this trait
if (tc >= pop * 2)
{
// Mutation substituted!
// Start a new one!
for (int i = 0; i < MAX_CRIT; i++)
{
mDNA[i][0][t] -= 1;
mDNA[i][1][t] -= 1;
}
maxval[t] -= 1;
mMutations++;
cout << "!" << mMutations;
if (mMutations >= mMaxMutations)
{
bStop = true;
}
}
}
// Every generation some newborn gets a beneficial mutation
int lYoung = mAges[0] - mAges[1];
// Complain if there are no newborns - something's wrong!
if (lYoung <= 0)
{
cout << "0:" << mAges[0] << " 1:" << mAges[1] <<
" d:" << lYoung << endl;
for (int i = 0; i < DEAD + 1; i++)
{
cout << i << ":" << mAges[i] << " ";
}
cout << endl << "Gap: ";
for (int i = 1; i < DEAD + 1; i++)
{
cout << i << ":" << (mAges[i-1] - mAges[i]) << " ";
}
cout << endl;
assert(lYoung > 0);
}
// Beneficial mutation
for (int j = 0; j < mGoodMutations; j++)
{
int lChoose = (choice(lYoung) + mAges[1]) % MAX_CRIT;
int lWhat = choice(mTraits * 2);
mDNA[lChoose][lWhat % 2][lWhat / 2] += 1;
}
// Detrimental mutation is much more common
for (int j = 0; j < mBadMutations; j++)
{
int lChoose = (choice(lYoung) + mAges[1]) % MAX_CRIT;
int lWhat = choice(mTraits * 2);
mDNA[lChoose][lWhat % 2][lWhat / 2] -= 1;
}
cout.flush();
return bStop;
}
// Check dates and times; pause if necessary
// Return true for pause, false for keep going.
bool DNAStore::pauseExecution()
{
time_t lNow = time(0);
tm* lTime = localtime(&lNow);
bool lAfterHours = ((lTime->tm_wday == 5) || (lTime->tm_wday == 6) ||
(lTime->tm_hour <= 7) || (lTime->tm_hour >= 20));
if ((mLastSave != lTime->tm_yday) && (lTime->tm_hour >= 8))
{
// Save a file at the end of each session
ofstream fsave( "big.out", ios::trunc | ios::out );
fsave << *this;
mLastSave = lTime->tm_yday;
}
if (!lAfterHours)
{
ifstream ifcheck( ".after_hours", ios::in );
if ( !ifcheck.fail() )
{
lAfterHours = true;
}
ifcheck.close();
}
// After hours is when we run.
if (!lAfterHours)
{
// Now sleep until the next time check
sleep(10);
}
return !lAfterHours;
}
// Run a trial for this setup
void DNAStore::runTrial()
{
time_t lStart = time(0);
cout << "Increase: " << mInc << " Fitness: " << mFitness <<
" Gen: " << GEN << " Pop: " << CRITTERS <<
" Good: " << mGoodMutations << " Bad: " << mBadMutations <<
" Traits: " << mTraits << " Expect: " << mMaxMutations << endl;
bool bStop = false;
while (!bStop)
{
if (!pauseExecution())
{
nextYear();
if (mYears % GEN == 0)
{
bStop = analyze();
}
if (mYears % 10000 == 0)
{
cout << "@" << (mYears / GEN);
}
if ((mYears / GEN) >= 500000)
{
bStop = true;
}
}
}
int lGen = mYears / GEN;
time_t lFinish = time(0);
double duration = lFinish - lStart;
cout << endl << "Time taken " << duration <<
"(sec) Per generation: " << (duration / lGen) << endl;
cout << "Increase: " << mInc << " Fitness: " << mFitness <<
" Gen: " << GEN << " Pop: " << CRITTERS <<
" Good: " << mGoodMutations <<
" Bad: " << mBadMutations << endl;
cout << "Generations: " << lGen << " Mutations: " << mMutations <<
" Rate:" << ((lGen + mMutations - 1) / mMutations) <<
" Traits: " << mTraits << endl;
ofstream report( "result.out", ios::out | ios::app);
report << "Increase: " << mInc << " Fitness: " << mFitness <<
" Gen: " << GEN << " Pop: " << CRITTERS <<
" Good: " << mGoodMutations <<
" Bad: " << mBadMutations << endl;
report << "Generations: " << lGen << " Mutations: " << mMutations <<
" Rate:" << ((lGen + mMutations - 1) / mMutations) <<
" Traits: " << mTraits << endl;
}
void DNAStore::summary()
{
cout << "Increment:" << mInc << endl;
cout << "Fitness:" << mFitness << endl;
cout << "Good Mutations:" << mGoodMutations << endl;
cout << "Bad Mutations:" << mBadMutations << endl;
cout << "Years:" << mYears << endl;
cout << "Traits:" << mTraits << endl;
cout << "Mutations:" << mMutations << endl;
cout << "MaxMutations:" << mMaxMutations << endl;
cout << "Age Markers: ";
for (int i = 0; i < DEAD + 1; i++)
{
cout << mAges[i] << " ";
}
cout << endl;
}
// Print some mutation info
void DNAStore::printStats()
{
summary();
int count[260];
cout << "Increase: " << mInc << " Fitness: " << mFitness <<
" Gen: " << GEN << " Pop: " << CRITTERS <<
" Good: " << mGoodMutations << " Bad: " << mBadMutations <<
" Traits: " << mTraits << " Expect: " << mMaxMutations << endl;
cout << "Pop size limits:" << mAges[0] <<
" and " << mAges[DEAD] << ", total " <<
(mAges[0] - mAges[DEAD])<< endl;
cout << "Counting DNA types" << endl;
for (int i = 0; i < 260; i++)
{
count[i] = 0;
}
int lS = mAges[DEAD];
int lF = mAges[0];
for (int p = lS; p < lF; p++)
{
int i = p % MAX_CRIT;
if (p % 10000 == 0)
{
cout << ".";
}
for (int j = 0; j < COPIES; j++)
{
for (int k = 0; k < TRAITS; k++)
{
int index = mDNA[i][j][k] + 128;
count[index]++;
}
}
}
cout << endl;
int sum = 0;
for (int i = 0; i < 260; i++)
{
if (count[i] != 0)
{
cout << "(" << (i-128) << ") " << count[i] << endl;
sum += count[i];
}
}
cout << "Sum is " << sum << " expected " <<
((lF - lS) * COPIES * TRAITS) << endl;
cout << "Population is " << (mAges[0] - mAges[DEAD]) << endl;
}
MAIN.CPP
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include "DNASTORE.H"
#include <time.h>
#include <assert.h>
void runTrial(DNAStore* aWorld, int aTraits, int aCount,
int aGood, double aInc, double aFitness)
{
// Don't allow excessive bad mutations, but require at least
// one to match each beneficial.
// Each generation is about CRITTERS / 20; max out the harmful
// mutations at 1 in 50.
int bad = DNAStore::CRITTERS / 1000;
if (bad < aGood)
{
bad = aGood;
}
aWorld->reset(aTraits, aCount, aGood, bad, aInc, aFitness);
aWorld->runTrial();
}
void runTrials(DNAStore* aWorld, int aGood, double aInc, double aFitness)
{
runTrial(aWorld, 1, 10, aGood, aInc, aFitness);
runTrial(aWorld, 3, 30, aGood, aInc, aFitness);
runTrial(aWorld, 7, 70, aGood, aInc, aFitness);
runTrial(aWorld, 14, 140, aGood, aInc, aFitness);
runTrial(aWorld, 23, 300, aGood, aInc, aFitness);
}
// Not really beautiful yet; if you have to stop half way through a
// run, only the trial you were running gets restarted.
// If that happens, you need to rewrite and recompile.
int main()
{
srand( (unsigned)time( NULL ) );
DNAStore* world = new DNAStore();
ifstream ifsave( "big.out", ios::in );
if ( !ifsave.fail() )
{
ifsave >> *world;
world->printStats();
world->runTrial();
}
else
{
runTrials(world, 1, 2, 2);
runTrials(world, 1, 2, 1.1);
runTrials(world, 1, 1.1, 2);
runTrials(world, 1, 1.1, 1.1);
}
return 0;
}
No comments:
Post a Comment