diff --git a/TREES_code/bit_vector.cpp b/TREES_code/bit_vector.cpp
new file mode 100644
index 0000000..f9e9e87
--- /dev/null
+++ b/TREES_code/bit_vector.cpp
@@ -0,0 +1,96 @@
+//
+// TREES -
+// A TRait-based Eco-Evolutionary Simulation tool
+// Copyright (C) 2017 Jörgen Ripa
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see .
+//
+// Contact: jorgen.ripa@biol.lu.se
+//
+
+#include
+
+#include "bit_vector.h"
+
+bit_vector::bit_vector() { // default constructor
+ current_int = 0;
+ next_position = 0;
+}
+
+void bit_vector::reserve(size_type bit_count) {
+ size_type int_count = bit_count/(sizeof(store_type)*8) + 1;
+ std::vector:: reserve(int_count);
+}
+
+void bit_vector::push_back(bool bit) {
+ if (bit) {
+ store_type mask = store_type(1)<::push_back(current_int); // store in base class
+ current_int = 0;
+ next_position = 0;
+ }
+}
+
+bit_vector::size_type bit_vector::get_total_bit_count() {
+ return size()*sizeof(store_type)*8 + next_position;
+}
+
+bit_vector::size_type bit_vector::get_total_byte_count() {
+ return size()*sizeof(store_type);
+}
+
+void bit_vector::write_to_file(oSimfile& osf) {
+ // Write size of bitstream, followed by the bits themselves
+ osf.write((uint32_t)get_total_bit_count());
+ osf.write((uint32_t)sizeof(store_type));
+ if(size()>0) {
+ osf.writeArray(&at(0), (int)size());
+ }
+ if (next_position>0) {
+ osf.write(current_int);
+ }
+}
+
+void bit_vector::read_from_file(iSimfile &isf) {
+ clear();
+ uint32_t total_bits = isf.read();
+ uint32_t chunk_size = isf.read(); // size of sore_type (in bytes)
+ assert(chunk_size==sizeof(store_type));
+ uint32_t total_bytes = total_bits / (sizeof(store_type)*8);
+ next_position = total_bits % (sizeof(store_type)*8); // rest bits
+ if(total_bytes>0) {
+ isf.readArray(*this, total_bytes);
+ }
+ if (next_position>0) {
+ current_int = isf.read();
+ } else {
+ current_int = 0;
+ }
+}
+
+bool bit_vector::get_bit(size_type bit_pos) {
+ size_type int_pos = bit_pos/(sizeof(store_type)*8);
+ store_type the_int;
+ if (int_pos0;
+}
diff --git a/TREES_code/bit_vector.h b/TREES_code/bit_vector.h
new file mode 100644
index 0000000..6da89b7
--- /dev/null
+++ b/TREES_code/bit_vector.h
@@ -0,0 +1,47 @@
+//
+// TREES -
+// A TRait-based Eco-Evolutionary Simulation tool
+// Copyright (C) 2017 Jörgen Ripa
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see .
+//
+// Contact: jorgen.ripa@biol.lu.se
+//
+
+#ifndef bit_vector_hpp
+#define bit_vector_hpp
+
+#include
+#include
+#include
+#include "simfile.h"
+
+typedef uint64_t store_type;
+
+class bit_vector : protected std::vector {
+ protected:
+ store_type current_int;
+ int next_position;
+
+ public:
+ bit_vector(); // default constructor
+ void reserve(size_type bit_count);
+ void push_back(bool bit);
+ size_type get_total_bit_count();
+ size_type get_total_byte_count();
+ void write_to_file(oSimfile& osf);
+ void read_from_file(iSimfile& isf);
+ bool get_bit(size_type i);
+};
+#endif /* bit_vector_hpp */
diff --git a/TREES_code/build.sh b/TREES_code/build.sh
index c0799b6..82e36a9 100644
--- a/TREES_code/build.sh
+++ b/TREES_code/build.sh
@@ -2,25 +2,26 @@
# build.sh
# TREES
-# Created by Jörgen Ripa
+# Created by Jörgen Ripa on 26/6 -14.
g++ -c -std=c++0x main.cpp
g++ -c -std=c++0x simulation.cpp
-g++ -c -std=c++0x genetics.cpp
+g++ -c -std=c++0x population.cpp
g++ -c -std=c++0x trait.cpp
-g++ -c -std=c++0x space.cpp
+g++ -c -std=c++0x transform.cpp
g++ -c -std=c++0x fitness.cpp
g++ -c -std=c++0x mating.cpp
-g++ -c -std=c++0x population.cpp
+g++ -c -std=c++0x space.cpp
+g++ -c -std=c++0x genetics.cpp
g++ -c -std=c++0x geneTracking.cpp
+g++ -c -std=c++0x genotype_phenotype_map.cpp
g++ -c -std=c++0x parameterFile.cpp
-g++ -c -std=c++0x sample.cpp
g++ -c -std=c++0x simfile.cpp
+g++ -c -std=c++0x bit_vector.cpp
g++ -c -std=c++0x fastmath.cpp
-g++ -c -std=c++0x mytime.cpp
g++ -c -std=c++0x randgen.cpp
-g++ -c -std=c++0x transform.cpp
+g++ -c -std=c++0x mytime.cpp
-g++ -o TREES main.o simulation.o genetics.o trait.o space.o fitness.o population.o mating.o geneTracking.o sample.o parameterFile.o simfile.o fastmath.o randgen.o mytime.o transform.o
+g++ -o TREES main.o simulation.o population.o trait.o transform.o fitness.o mating.o space.o genetics.o geneTracking.o genotype_phenotype_map.o parameterFile.o simfile.o bit_vector.o fastmath.o randgen.o mytime.o
rm *.o
diff --git a/TREES_code/fitness.cpp b/TREES_code/fitness.cpp
index ebb7069..dea7f7c 100644
--- a/TREES_code/fitness.cpp
+++ b/TREES_code/fitness.cpp
@@ -19,6 +19,8 @@
// Contact: jorgen.ripa@biol.lu.se
//
+// Fitness calculations are generally done with double precision,
+// although trait values may be in single precision.
#include
#include
@@ -41,6 +43,7 @@ StabilizingSelection::StabilizingSelection(Population& p, ParameterFile& pf) :
Fitness(p) {
std::string zTraitname = pf.getString("trait");
zTrait = pop.findTrait(zTraitname);
+ optimum = pf.getDouble("optimal_value");
cost_coefficient = pf.getDouble("cost_coefficient");
cost_exponent = pf.getDouble("cost_exponent");
}
@@ -49,17 +52,18 @@ StabilizingSelection::~StabilizingSelection() {}
void StabilizingSelection::aggregateFitness(std::vector& fitness) {
//double* fi = &fitness[0];
+ Fastexp fexp;
for (int i=0; i0) {
double dzsum = 0;
- for (int d=0; dgetDims(); ++d) {
- double z = zTrait->traitValue(i,d);
- dzsum += z*z;
+ for (int d=0; dget_dims(); ++d) {
+ double dz = zTrait->traitValue(i,d) - optimum;
+ dzsum += dz*dz;
}
if(cost_exponent!=2) {
dzsum = fastabspow(dzsum,cost_exponent/2);
}
- fitness[i] *= std::max(0.0,(1-cost_coefficient*dzsum));
+ fitness[i] *= fexp.exp(-cost_coefficient*dzsum);
}
}
}
@@ -97,7 +101,7 @@ void DensityDependence::aggregateFitness(std::vector& fitness) {
// Within patch case:
} else if (pop.getSpace().isDiscrete() && s_space==0) {
- DiscreteSpace& theSpace = dynamic_cast(pop.getSpace());
+ Discrete_space& theSpace = dynamic_cast(pop.getSpace());
int n = pop.size();
// Vector of local patch densities (linear indexing):
@@ -152,13 +156,13 @@ ResourceLandscape::~ResourceLandscape() {
double ResourceLandscape::getTraitDist2( int ind1, int ind2) { // squared distance in trait space
double dx2 = 0;
- int dims = xTrait->getDims();
+ int dims = xTrait->get_dims();
if (dims==1) {
dx2 = getX(ind1) - getX(ind2);
dx2 *= dx2;
} else {
- double *xp1 = &getX(ind1,0);
- double *xp2 = &getX(ind2,0);
+ traitType *xp1 = &getX(ind1,0);
+ traitType *xp2 = &getX(ind2,0);
for (int d=0; d& fitness)
{
double twosa2 = 2*sa*sa;
int n = pop.size();
- int xdims = xTrait->getDims();
+ int xdims = xTrait->get_dims();
int sdims = pop.getSpace().getDims();
std::vector comp;
@@ -191,7 +195,7 @@ void ResourceLandscape::aggregateFitness(std::vector& fitness)
if (i& fitness)
return;
} else if (pop.getSpace().isDiscrete() && s_space==0) {
// Special case 2: within patch competition
- DiscreteSpace& theSpace = dynamic_cast(pop.getSpace());
+ Discrete_space& theSpace = dynamic_cast(pop.getSpace());
for (int i=0; i& fitness)
}
}
if (xdims==1) {
- double *xip = &getX(0,0);
+ traitType *xip = &getX(0,0);
for (int i=0; i& fitness)
++xip;
}
} else {
- double *xip = &getX(0,0);
+ traitType *xip = &getX(0,0);
for (int i=0; i& fitness) {
}
}
} else { // we can assume discrete space
- DiscreteSpace& theSpace = dynamic_cast(pop.getSpace());
+ Discrete_space& theSpace = dynamic_cast(pop.getSpace());
for (int patchi=0; patchitraitValue(individual); }
- inline double& getX(int individual, int dim) { return xTrait->traitValue(individual,dim); }
+ inline traitType& getX(int individual) { return xTrait->traitValue(individual); }
+ inline traitType& getX(int individual, int dim) { return xTrait->traitValue(individual,dim); }
double getTraitDist2( int ind1, int ind2); // squared distance in trait space
public:
ResourceLandscape(Population& pop, ParameterFile& pf);
@@ -100,7 +101,7 @@ class DiscreteResources : public Fitness {
double ta; // trade-off
double cmin; // minimal consumption level for status quo (standard = 1)
Trait* xTrait; // resource adaptation trait (only first dimension is used)
- inline double& getX(int individual) { return xTrait->traitValue(individual,0); }
+ inline traitType& getX(int individual) { return xTrait->traitValue(individual,0); }
public:
DiscreteResources(Population& pop, ParameterFile& pf);
virtual ~DiscreteResources();
diff --git a/TREES_code/geneTracking.cpp b/TREES_code/geneTracking.cpp
index bced18d..76c9e5b 100644
--- a/TREES_code/geneTracking.cpp
+++ b/TREES_code/geneTracking.cpp
@@ -27,7 +27,7 @@
#include "simfile.h"
#include "types.h"
-Gene::Gene(idType parent, timeType time, double effect) {
+Gene::Gene(id_type parent, time_type time, double effect) {
id = 0;
this->parent = parent;
birthTime = time;
@@ -47,11 +47,11 @@ Gene::Gene() {
sampled = false;
}
-void Gene::addChild(idType child) {
+void Gene::addChild(id_type child) {
children.push_back(child);
}
-void Gene::removeChild(idType child) {
+void Gene::removeChild(id_type child) {
for (childIter c=children.begin(); c!=children.end(); ++c) {
if (*c==child) {
children.erase(c);
@@ -62,19 +62,32 @@ void Gene::removeChild(idType child) {
exit(1);
}
-Simfile& operator<<(Simfile& sf, Gene& g) {
- sf.write(g.id);
- sf.write(g.parent);
- sf.write(g.birthTime);
- sf.write(g.deathTime);
- sf.write((float)g.geneticEffect);
- sf.write((int)g.children.size());
- if (g.children.size()>0) {
- sf.writeArray(&g.children[0], (int)g.children.size());
+void Gene::writeToFile(oSimfile& sf) {
+ sf.write(id);
+ sf.write(parent);
+ sf.write(birthTime);
+ sf.write(deathTime);
+ sf.write(geneticEffect);
+ sf.write((int)children.size());
+ if (children.size()>0) {
+ sf.writeArray(&children[0], (int)children.size());
}
- return sf;
}
+Gene::Gene(iSimfile& isf) {
+ id = isf.read();
+ parent = isf.read();
+ birthTime = isf.read();
+ deathTime = isf.read();
+ geneticEffect = isf.read();
+ sampled = true;
+ int childrenCount = isf.read();
+ children.clear();
+ children.reserve(childrenCount);
+ for (int ic=0; ic());
+ }
+}
GeneList::GeneList() {
clear();
@@ -82,10 +95,28 @@ GeneList::GeneList() {
nextId = 1; // 0 is reserved for invalid Genes
}
+GeneList::GeneList(iSimfile& isf) {
+ int count = isf.read();
+ clear();
+ reserve(count);
+ for (int gi=0; gi();
+}
+
GeneList::~GeneList() {
}
-idType GeneList::addGene(idType parent, timeType time, double effect) {
+void GeneList::writeGenes(oSimfile& sf) {
+ sf.write((int)size());
+ for( Gene& g: *this) {
+ g.writeToFile(sf);
+ }
+ sf.write(nextId);
+}
+
+id_type GeneList::addGene(id_type parent, time_type time, double effect) {
Gene newA( parent, time, effect);
newA.id = nextId++;
push_back(newA);
@@ -95,9 +126,9 @@ idType GeneList::addGene(idType parent, timeType time, double effect) {
return newA.id;
}
-int GeneList::getGeneIndex(idType id) {
+int GeneList::getGeneIndex(id_type id) {
for (int a=0; a=0) {
return at(getGeneIndex(id));
@@ -127,7 +158,7 @@ void GeneList::pruneGenes(std::vector& alive) {
}
++iw;
} else { // remove dead gene from parent's list of children
- idType parent = at(ir).parent;
+ id_type parent = at(ir).parent;
if (parent>0) { // zero is the parent of the root
int parent_i = getGeneIndex(parent);
// If parent still in the list and active:
diff --git a/TREES_code/geneTracking.h b/TREES_code/geneTracking.h
index 59600ac..1e16492 100644
--- a/TREES_code/geneTracking.h
+++ b/TREES_code/geneTracking.h
@@ -28,42 +28,48 @@
#include
#include "types.h"
-typedef uint64_t idType;
-class Simfile;
+class oSimfile;
+class iSimfile;
class Gene {
friend class GeneList;
friend class Genetics;
protected:
- idType id; // the id is assigned by the GeneList
- idType parent;
- timeType birthTime, deathTime;
- double geneticEffect;
- std::vector children;
+ id_type id; // the id is assigned by the GeneList
+ id_type parent;
+ time_type birthTime, deathTime;
+ float geneticEffect;
+ std::vector children;
bool sampled; // sampled genes should not be pruned
- typedef std::vector::iterator childIter;
+ typedef std::vector::iterator childIter;
public:
- Gene(); // somehow, this is needed
- Gene(idType parent, timeType time, double effect);
- void addChild(idType child);
- void removeChild(idType child);
- std::vector& getChildren() { return children;}
- friend Simfile& operator<<(Simfile& os, Gene& a);
+ Gene(); // needed for lists, etc.
+ Gene(id_type parent, time_type time, double effect);
+ Gene(iSimfile& isf);
+ void addChild(id_type child);
+ void removeChild(id_type child);
+ std::vector& getChildren() { return children;}
+ void writeToFile(oSimfile& osf);
};
-Simfile& operator<<(Simfile& os, Gene& a);
+//oSimfile& operator<<(oSimfile& os, Gene& a);
class GeneList : public std::vector {
protected:
- idType nextId;
+ id_type nextId;
public:
GeneList();
+ GeneList(iSimfile& isf);
~GeneList();
- idType addGene(idType parent, timeType time, double effect); // returns new Allele id
- Gene& getGene(idType id);
- int getGeneIndex(idType id);
+ id_type addGene(id_type parent, time_type time, double effect); // returns new Allele id
+ Gene& getGene(id_type id);
+ int getGeneIndex(id_type id);
void pruneGenes(std::vector& alive);
+ void writeGenes(oSimfile& osf);
};
+//oSimfile& operator<<(oSimfile& os, GeneList& list);
+
+
#endif /* defined(__Species__geneTracking__) */
diff --git a/TREES_code/genetics.cpp b/TREES_code/genetics.cpp
index 840e309..50e3c4f 100644
--- a/TREES_code/genetics.cpp
+++ b/TREES_code/genetics.cpp
@@ -27,7 +27,7 @@
#include "genetics.h"
#include "population.h"
-#include "sample.h"
+#include "genotype_phenotype_map.h"
#include "randgen.h"
#include "simfile.h"
@@ -39,10 +39,10 @@ pop(p) {
geneLists.clear();
if (geneTracking()) {
- G1id.reserve(10000);
- G2id.reserve(10000);
- newG1id.reserve(10000);
- newG2id.reserve(10000);
+ G1id.reserve(1,10000);
+ G2id.reserve(1,10000);
+ newG1id.reserve(1,10000);
+ newG2id.reserve(1,10000);
}
}
@@ -53,45 +53,47 @@ Genetics::~Genetics() {}
bool Genetics::geneTracking()
{ return pop.isTrackingGenes(); }
-void Genetics::initializeTracking(int n0) {
- idType firstId;
+void Genetics::initializeTracking() {
+ id_type firstId=0;
for (int l=0; lsampled) {
gp->sampled = true;
@@ -107,15 +109,58 @@ void Genetics::addGeneIDsToSample(Sample& s) {
}
}
+void Genetics::writeGeneLists(oSimfile &sf) {
+ //One list per locus:
+ for (GeneList& list : geneLists) {
+ list.writeGenes(sf);
+ }
+}
+
+void Genetics::readGeneLists(iSimfile &sf) {
+ //One list per locus
+ geneLists.clear();
+ geneLists.reserve(loci);
+ for (int li=0; li& glists) {
+ resumeGenomesFromCheckpoint(gsample);
+ if (pop.isTrackingGenes()) {
+ G1id = gsample->G1id;
+ G2id = gsample->G2id;
+ geneLists = glists; // deep copy?
+ }
+}
+
+//void Genetics::readGeneListsToCheckpoint(Checkpoint& cp, iSimfile& isf) {
+// std::vector& gls = cp.geneListsCopy;
+// gls.clear();
+// gls.reserve(loci);
+// for (int li=0; li alive;
alive.assign(list.size(), false);
// find all current alleles in population:
for (int i=0; i0) {
int parent_i = list.getGeneIndex(parent);
if (!alive.at(parent_i)) {
@@ -149,41 +194,95 @@ void Genetics::pruneGeneLists() {
void Genetics::compactData(std::vector& alive) {
if (geneTracking()) {
- // compact allele id information:
- int iwrite=0; // position for write
- int iread; // position for read
- for (iread=0; ireadiwrite) {
- for (int li=0; li G1id;
+ std::vector G2id;
+ Genetics_sample(Genetics& G);
+ public:
+ virtual void write_to_file(oSimfile& osf);
+};*/
+
+Genetics_sample::Genetics_sample() {
+ loci = 0;
+ pop_size = 0;
+ gene_tracking = false;
+ G1id.clear();
+ G2id.clear();
+}
+
+Genetics_sample::Genetics_sample(Genetics& G) {
+ loci = G.getLoci();
+ pop_size = G.pop.size();
+ gene_tracking = G.pop.isTrackingGenes();
+ if (gene_tracking) {
+ G1id = G.G1id;
+ G2id = G.G2id;
+ G.mark_all_living_genes_sampled();
+ }
+}
+
+void Genetics_sample::write_to_file(oSimfile& osf) {
+ osf.write(loci);
+ osf.write(pop_size);
+ osf.write((char)gene_tracking);
+ if (gene_tracking) {
+ G1id.write_to_file(osf);
+ G2id.write_to_file(osf);
+ }
+}
+
+Genetics_sample::Genetics_sample(iSimfile& isf) {
+ loci = isf.read();
+ pop_size = isf.read();
+ gene_tracking = (bool)isf.read();
+ if (gene_tracking) {
+ G1id.read_from_file(isf);
+ G2id.read_from_file(isf);
+ }
+}
+
+Genetics_sample::~Genetics_sample() {
+}
+
+Genetics_sample* Genetics_sample::create_from_file(iSimfile& isf) { // polymorphic read
+ Genetics_types type = (Genetics_types)isf.read();
+ switch (type) {
+ case Genetics_types::Continuous_type: {
+ ContinuousAlleles_sample* CAsam = new ContinuousAlleles_sample(isf);
+ return CAsam;
+ break;
+ }
+ case Genetics_types::Diallelic_type: {
+ Diallelic_sample* Dsam = new Diallelic_sample(isf);
+ return Dsam;
+ break;
}
+ default:
+ std::cout << "Wrong genetics type\n";
+ exit(1);
+ return NULL;
}
}
+//////////////////////////////////////////////////////////////////////////
+// ContinuousAlleles
+//////////////////////////////////////////////////////////////////////////
+
ContinuousAlleles::ContinuousAlleles(Population& pop, ParameterFile& pf) :
Genetics(pop, pf) {
- G1.reserve(10000);
- G2.reserve(10000);
- newG1.reserve(10000);
- newG2.reserve(10000);
+ G1.reserve(1,10000);
+ G2.reserve(1,10000);
+ newG1.reserve(1,10000);
+ newG2.reserve(1,10000);
}
ContinuousAlleles::~ContinuousAlleles() {
@@ -196,41 +295,38 @@ double ContinuousAlleles::getEffect2(int individual, int locus) {
return getGene2(individual, locus);
}
-void ContinuousAlleles::initialize(int n0) {
+Genotype_Phenotype_map* ContinuousAlleles::create_GP_map(Trait& T, int loci_per_dim, ParameterFile& pf) {
+ return new ContinuousAlleles_GP_map(T, *this, loci_per_dim, pf);
+}
+
+void ContinuousAlleles::initialize() {
double P_No_Mutations = pow(1-Pmut,loci);
P_At_least_one_mut = 1 - P_No_Mutations;
- G1.assign(loci*n0,0);
- G2.assign(loci*n0,0);
+ G1.assign(loci,pop.size(),0);
+ G2.assign(loci,pop.size(),0);
newG1.clear();
newG2.clear();
}
-void ContinuousAlleles::setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim ) {
- // sets first locus to initial value, the rest remain zero
- for (int ind=0; ind-1 && nG[mutlocus]<1);
}
} else { // no gene tracking
- double* g1 = &getGene1(parent,0);
- double* g2 = &getGene2(parent,0);
- double* target;
+ float* g1 = &getGene1(parent,0);
+ float* g2 = &getGene2(parent,0);
+ float* target;
if (targetHaplo==1) {
- target = &newG1.at(loci*newCount);
+ target = &newG1(0,newCount);
} else {
- target = &newG2.at(loci*newCount);
+ target = &newG2(0,newCount);
}
bitGenerator bG;
- int haplo = 1 + int(rand1()*2); // copy from parent haplotype 1 or 2
- double* target_it = target;
+ float* target_it = target;
for (int a=0; a-1 && nG[a]<1);
- if( bG.nextBit()){
- haplo = 3 - haplo; // changes 1 to 2 and vice versa
- }
}
// Possibly mutate:
while (rand1() < P_At_least_one_mut) {
int mutlocus = rand1()*loci;
- target[mutlocus] += randn(); // mutations are always of SD=1
- //assert(nG[mutlocus]>-1 && nG[mutlocus]<1);
+ target[mutlocus] += doubleExp(1); // mutations have variance 1
}
}
}
-double ContinuousAlleles::getGeneSum(int individual, int startlocus, int endlocus) {
- // no gene interactions:
- double sum = 0;
- double* g1 = &getGene1(individual,startlocus);
- double* g2 = &getGene2(individual,startlocus);
- for (int li=startlocus; li& alive) {
Genetics::compactData(alive);
- // compact genome:
- int iwrite=0; // position for write
- int iread; // position for read
- for (iread=0; ireadiwrite) {
- memcpy(&getGene1(iwrite, 0), &getGene1(iread, 0), loci*sizeof(G1[0]));
- memcpy(&getGene2(iwrite, 0), &getGene2(iread, 0), loci*sizeof(G1[0]));
-/* for (int a=0; a(gsample);
+ G1 = CA_sample->G1;
+ G2 = CA_sample->G2;
+}
+
+int ContinuousAlleles::add_loci(int L) {
+ int start_locus = loci;
+ loci += L;
+ return start_locus;
+}
+
+Genetics_sample* ContinuousAlleles::make_sample() {
+ return new ContinuousAlleles_sample(*this);
+}
+
+/*class ContinuousAlleles_sample : public Genetics_sample {
+ std::vector G1;
+ std::vector G2;
+
+ public:
+ ContinuousAlleles_sample(ContinuousAlleles& G);
+ virtual void write_to_file(oSimfile& osf);
+};*/
+
+ContinuousAlleles_sample::ContinuousAlleles_sample(ContinuousAlleles& G) : Genetics_sample(G) {
+ G1 = G.G1;
+ G2 = G.G2;
+}
+
+void ContinuousAlleles_sample::write_to_file(oSimfile& osf) {
+ // First specify type, convenient for reading:
+ osf.write((char)Genetics_types::Continuous_type);
+ // Let base class write basic data and gene tracking data, if applicable
+ Genetics_sample::write_to_file(osf);
+ G1.write_to_file(osf);
+ G2.write_to_file(osf);
}
-void ContinuousAlleles::addToSample(Sample& s) {
- s.addData(new FloatData(G1));
- s.addData(new FloatData(G2));
+ContinuousAlleles_sample::ContinuousAlleles_sample(iSimfile& isf) :
+Genetics_sample(isf) {
+ G1.read_from_file(isf);
+ G2.read_from_file(isf);
}
+// This virtual destructor is essential.
+ContinuousAlleles_sample::~ContinuousAlleles_sample() {
+}
-////////////////////////////////
+///////////////////////////////////////////////////////////////////
+// Diallelic
+///////////////////////////////////////////////////////////////////
Diallelic::Diallelic(Population& pop, ParameterFile& pf) :
Genetics(pop, pf) {
- G1.reserve(10000);
- G2.reserve(10000);
- newG1.reserve(10000);
- newG2.reserve(10000);
+ G1.reserve(1,10000);
+ G2.reserve(1,10000);
+ newG1.reserve(1,10000);
+ newG2.reserve(1,10000);
}
Diallelic::~Diallelic() {
}
double Diallelic::getEffect1(int individual, int locus) {
- return (double)getGene1(individual, locus);
+ return double(getGene1(individual, locus))/2.0;
}
double Diallelic::getEffect2(int individual, int locus) {
- return (double)getGene2(individual, locus);
+ return double(getGene2(individual, locus))/2.0;
}
+Genotype_Phenotype_map* Diallelic::create_GP_map(Trait& T, int loci_per_dim, ParameterFile& pf) {
+ return new Diallelic_GP_map(T, *this, loci_per_dim, pf);
+}
-void Diallelic::initialize(int n0) {
+void Diallelic::initialize() {
double P_No_Mutations = pow(1-Pmut,loci);
P_At_least_one_mut = 1 - P_No_Mutations;
- G1.assign(loci*n0,-1);
- G2.assign(loci*n0,-1);
+ G1.assign(loci,pop.size(),-1);
+ G2.assign(loci,pop.size(),-1);
newG1.clear();
newG2.clear();
-
-}
-
-void Diallelic::setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim ) {
- //int traitLoci = dims*lociPerDim;
- int plusCount = (int)((value+lociPerDim)/2.0);
- if (plusCount<0 || plusCount>lociPerDim) {
- std::cout << "Illegal initial trait value : " << value << '\n';
- exit(0);
- }
- // set the first loci to +1, the rest remain -1
- for (int ind=0; ind& alive) {
Genetics::compactData(alive);
// compact genome:
- int iwrite=0; // position for write
- int iread; // position for read
- for (iread=0; ireadiwrite) {
- memcpy(&getGene1(iwrite, 0), &getGene1(iread, 0), loci*sizeof(G1[0]));
- memcpy(&getGene2(iwrite, 0), &getGene2(iread, 0), loci*sizeof(G1[0]));
+ G1.compact_data(alive);
+ G2.compact_data(alive);
+}
+
+void Diallelic::resumeGenomesFromCheckpoint(Genetics_sample* gsample) {
+ Diallelic_sample* D_sample = dynamic_cast(gsample);
+
+ // copy stored bitvectors to genes
+ // (0/1) correspond to (-1/+1)
+ initialize(); // sets all genes to -1
+ for (int li=0; liG1.get_bit(li + ind*loci)) {
+ G1(li,ind) = +1;
+ }
+ if (D_sample->G2.get_bit(li + ind*loci)) {
+ G2(li,ind) = +1;
}
- ++iwrite;
}
}
- G1.resize(iwrite*loci);
- G2.resize(iwrite*loci);
}
-void Diallelic::addToSample(Sample& s) {
- s.addData(new IntData(G1));
- s.addData(new IntData(G2));
+int Diallelic::add_loci(int L) {
+ int start_locus = loci;
+ loci += L;
+ return start_locus;
}
+
+Genetics_sample* Diallelic::make_sample() {
+ return new Diallelic_sample(*this);
+}
+
+Diallelic_sample::Diallelic_sample(Diallelic& G) : Genetics_sample(G) {
+ G1.reserve(loci*pop_size);
+ G2.reserve(loci*pop_size);
+ for (int ind=0; ind0);
+ G2.push_back(G.getGene2(ind, locus)>0);
+ }
+ }
+}
+
+void Diallelic_sample::write_to_file(oSimfile& osf) {
+ osf.write((char)Genetics_types::Diallelic_type);
+ // First let base class write basic data and gene tracking data, if applicable
+ Genetics_sample::write_to_file(osf);
+ G1.write_to_file(osf);
+ G2.write_to_file(osf);
+}
+
+Diallelic_sample::Diallelic_sample(iSimfile& isf) :
+Genetics_sample(isf) {
+ G1.read_from_file(isf);
+ G2.read_from_file(isf);
+}
+
+// This virtual destructor is essential.
+Diallelic_sample::~Diallelic_sample() {
+}
+
+/////////////////////////////////////////
+// Omnigenic genetics model
+/////////////////////////////////////////
+Omnigenic::Omnigenic(Population& pop, ParameterFile& pf) :
+ContinuousAlleles(pop, pf) {
+ loci = pf.getPositiveInt("loci"); // This is only set once for this genetics module
+}
+
+Genotype_Phenotype_map* Omnigenic::create_GP_map(Trait& T, int loci_per_dim, ParameterFile& pf) {
+ return new Omnigenic_GP_map(T, *this, loci_per_dim, pf);
+}
+
diff --git a/TREES_code/genetics.h b/TREES_code/genetics.h
index b628b27..36e7750 100644
--- a/TREES_code/genetics.h
+++ b/TREES_code/genetics.h
@@ -25,20 +25,27 @@
#include
#include "parameterFile.h"
-#include "sample.h"
#include "geneTracking.h"
#include "types.h"
+#include "matrix.h"
+#include "bit_vector.h"
class Population;
+class Genotype_Phenotype_map;
+class Trait;
+class ParameterFile;
///////////////////////////////////////
// Genetics
// Base class for genetics modules
// Handles gene tacking, if activated
////////////////////////////////////////
+class Genetics_sample;
+
class Genetics {
- friend class Trait;
- friend class TraitConstant;
+// friend class Trait;
+// friend class TraitConstant;
+ friend class Genetics_sample;
protected:
Population& pop;
int loci; // the number of loci
@@ -49,14 +56,14 @@ class Genetics {
// Tracking:
std::vector geneLists; // Gene tables for tracking, one per locus
- void addGene(int locus, idType parent, double effect, timeType time);
- std::vector G1id, G2id; // mirrors of G1 and G2
- std::vector newG1id, newG2id; // mirrors of newG1 and newG2
- inline idType& getGene1id(int individual, int locus) {
- return G1id.at(individual*loci + locus);
+ void addGene(int locus, id_type parent, double effect, time_type time);
+ Matrix G1id, G2id; // mirrors of G1 and G2
+ Matrix newG1id, newG2id; // mirrors of newG1 and newG2
+ inline id_type& getGene1id(int individual, int locus) {
+ return G1id(locus,individual);
}
- inline idType& getGene2id(int individual, int locus) {
- return G2id.at(individual*loci + locus);
+ inline id_type& getGene2id(int individual, int locus) {
+ return G2id(locus,individual);
}
// protected constructor:
@@ -64,84 +71,164 @@ class Genetics {
public:
virtual ~Genetics();
//virtual bool isClonal(); // returns false by default
- void initializeTracking(int n0);
- virtual void initialize(int n0)=0;
- virtual void setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim )=0;
+ virtual Genotype_Phenotype_map* create_GP_map(Trait& T, int loci_per_dim, ParameterFile& pf)=0;
+ void initializeTracking();
+ virtual void initialize()=0;
+ //virtual void setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim )=0;
virtual void prepareNewGeneration(int size);
virtual void nextGeneration();
virtual void addChild(int mom, int dad)=0;
- virtual double getGeneSum(int individual, int startlocus, int endlocus)=0; // endlocus is exclusive!
+ //virtual double getGeneSum(int individual, int startlocus, int endlocus)=0; // endlocus is exclusive!
virtual void compactData(std::vector& alive);
- virtual void addToSample(Sample& s)=0;
+ virtual void resumeGenomesFromCheckpoint(Genetics_sample* gsample)=0;
+ virtual Genetics_sample* make_sample()=0;
+ void resumeFromCheckpoint(Genetics_sample* gsample, std::vector& glists);
- void addGeneIDsToSample(Sample& s);
int getLoci() { return loci;}
void pruneGeneLists();
bool geneTracking();
- void saveGeneLists(Simfile& sf);
+ void writeGeneLists(oSimfile& sf);
+ void readGeneLists(iSimfile& sf);
+ void mark_all_living_genes_sampled();
+ const std::vector& getGeneLists() { return geneLists; }
+};
+
+// Enumeration used for polymorphic file read/write
+enum Genetics_types {
+ Continuous_type,
+ Diallelic_type
+};
+
+class Genetics_sample {
+ friend class Genetics;
+ protected:
+ bool gene_tracking;
+ int loci;
+ int pop_size;
+ Matrix G1id;
+ Matrix G2id;
+ Genetics_sample();
+ Genetics_sample(Genetics& G);
+ Genetics_sample(iSimfile& isf); // base class read
+
+ public:
+ virtual ~Genetics_sample();
+ virtual void write_to_file(oSimfile& osf);
+ static Genetics_sample* create_from_file(iSimfile& isf); // polymorphic read
};
// Continuum of alleles genetics:
class ContinuousAlleles : public Genetics {
- friend class Trait;
+// friend class Trait;
+ friend class ContinuousAlleles_GP_map;
+ friend class ContinuousAlleles_sample;
+protected:
+ // The local coding is float, to save space
+ typedef float geneType;
public:
ContinuousAlleles(Population& pop, ParameterFile& pf);
virtual ~ContinuousAlleles();
- virtual void initialize(int n0);
- virtual void setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim );
+ virtual Genotype_Phenotype_map* create_GP_map(Trait& T, int loci_per_dim, ParameterFile& pf);
+ virtual void initialize();
+ //virtual void setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim );
virtual void prepareNewGeneration(int size);
virtual void nextGeneration();
virtual void addChild(int mom, int dad);
- virtual double getGeneSum(int individual, int startlocus, int endlocus); // endlocus is exclusive!
+ //virtual double getGeneSum(int individual, int startlocus, int endlocus); // endlocus is exclusive!
virtual void compactData(std::vector& alive);
- virtual void addToSample(Sample& s);
+ virtual void resumeGenomesFromCheckpoint(Genetics_sample* gsample);
+ virtual Genetics_sample* make_sample();
+ int add_loci(int L);
protected:
virtual double getEffect1(int individual, int locus);
virtual double getEffect2(int individual, int locus);
double P_At_least_one_mut;
- std::vector G1, G2;
+ Matrix G1, G2;
//int Gsize;
- std::vector newG1, newG2;
+ Matrix newG1, newG2;
//int newSize;
- inline double& getGene1(int individual, int locus) {
- return G1.at(individual*loci + locus);
+ inline geneType& getGene1(int individual, int locus) {
+ return G1(locus,individual);
}
- inline double& getGene2(int individual, int locus) {
- return G2.at(individual*loci + locus);
+ inline geneType& getGene2(int individual, int locus) {
+ return G2(locus,individual);
}
void produceGamete(int parent, int targetHaplo);
};
+class ContinuousAlleles_sample : public Genetics_sample {
+ friend class ContinuousAlleles;
+ protected:
+ Matrix G1;
+ Matrix G2;
+
+ public:
+ ContinuousAlleles_sample();
+ ContinuousAlleles_sample(ContinuousAlleles& G);
+ ContinuousAlleles_sample(iSimfile& isf);
+ virtual ~ContinuousAlleles_sample();
+ virtual void write_to_file(oSimfile& osf);
+};
-// Diallelic genetics (0/1 alleles):
+// Diallelic genetics (+-0.5 alleles):
class Diallelic : public Genetics {
- friend class Trait;
+ friend class Diallelic_GP_map;
+ friend class Diallelic_sample;
+protected:
+ // The local coding is +-1, although the effects are actually +-0.5.
+ // This is adjusted in 'getEffect1/2' and 'getGeneSum'.
+ typedef signed char geneType;
public:
Diallelic(Population& pop, ParameterFile& pf);
virtual ~Diallelic();
- virtual void initialize(int n0);
- virtual void setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim );
+ virtual Genotype_Phenotype_map* create_GP_map(Trait& T, int loci_per_dim, ParameterFile& pf);
+ virtual void initialize();
+ //virtual void setInitialValue(double value, int n, int startLocus, int dims, int lociPerDim );
virtual void prepareNewGeneration(int size);
virtual void nextGeneration();
virtual void addChild(int mom, int dad);
- virtual double getGeneSum(int individual, int startlocus, int endlocus); // endlocus is exclusive!
+ //virtual double getGeneSum(int individual, int startlocus, int endlocus); // endlocus is exclusive!
virtual void compactData(std::vector& alive);
- virtual void addToSample(Sample& s);
-
+ virtual void resumeGenomesFromCheckpoint(Genetics_sample* gsample);
+ virtual Genetics_sample* make_sample();
+ int add_loci(int L);
protected:
virtual double getEffect1(int individual, int locus);
virtual double getEffect2(int individual, int locus);
double P_At_least_one_mut;
- std::vector G1, G2;
- std::vector newG1, newG2;
- inline int& getGene1(int individual, int locus) {
- return G1.at(individual*loci + locus);
+
+ Matrix G1, G2;
+ Matrix newG1, newG2;
+ inline geneType& getGene1(int individual, int locus) {
+ return G1(locus,individual);
}
- inline int& getGene2(int individual, int locus) {
- return G2.at(individual*loci + locus);
+ inline geneType& getGene2(int individual, int locus) {
+ return G2(locus,individual);
}
void produceGamete(int parent, int targetHaplo);
};
+class Diallelic_sample : public Genetics_sample {
+ friend class Diallelic;
+ bit_vector G1; // compressed format of 0/1 bits
+ bit_vector G2;
+
+ public:
+ Diallelic_sample();
+ Diallelic_sample(Diallelic& G);
+ Diallelic_sample(iSimfile& isf);
+ virtual ~Diallelic_sample();
+ virtual void write_to_file(oSimfile& osf);
+};
+
+// Omnigenic model, complete pleitropy:
+class Omnigenic : public ContinuousAlleles {
+ // friend class Trait;
+ friend class Omnigenic_GP_map;
+public:
+ Omnigenic(Population& pop, ParameterFile& pf);
+ virtual Genotype_Phenotype_map* create_GP_map(Trait& T, int loci_per_dim, ParameterFile& pf);
+};
+
#endif /* defined(__Species__genetics__) */
diff --git a/TREES_code/genotype_phenotype_map.cpp b/TREES_code/genotype_phenotype_map.cpp
new file mode 100644
index 0000000..a013f35
--- /dev/null
+++ b/TREES_code/genotype_phenotype_map.cpp
@@ -0,0 +1,254 @@
+//
+// genotype_phenotype_map.cpp
+// TREES
+//
+// TREES -
+// A TRait-based Eco-Evolutionary Simulation tool
+// Copyright (C) 2017 Jörgen Ripa
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see .
+//
+// Contact: jorgen.ripa@biol.lu.se
+//
+
+#include
+#include "genotype_phenotype_map.h"
+#include "population.h"
+#include "genetics.h"
+#include "trait.h"
+#include "transform.h"
+#include "randgen.h"
+
+Genotype_Phenotype_map::Genotype_Phenotype_map(Trait& theTrait) : T(theTrait) {
+ transforms.clear();
+}
+
+Genotype_Phenotype_map::~Genotype_Phenotype_map() {
+ for (Transform* t : transforms) {
+ delete t;
+ }
+}
+
+double Genotype_Phenotype_map::do_all_transforms(double y) {
+ for (Transform* t : transforms) {
+ y = t->transform_single(y);
+ }
+ return y;
+}
+
+int Genotype_Phenotype_map::get_dims() {
+ return T.get_dims();
+}
+
+int Genotype_Phenotype_map::pop_size() {
+ return T.pop.size();
+}
+
+Trait& Genotype_Phenotype_map::get_trait() {
+ return T;
+}
+
+void Genotype_Phenotype_map::read_all_transforms(ParameterFile& pf) {
+ while (pf.get_next_name()=="transform") {
+ std::string tname = pf.get_next_value_string_lower();
+ if (tname=="linear") {
+ transforms.push_back(new LinearTransform(pf));
+ } else if (tname=="abs") {
+ transforms.push_back(new AbsTransform(pf));
+ } else if (tname=="logistic") {
+ transforms.push_back(new LogisticTransform(pf));
+ } else if (tname=="normal_deviate") {
+ transforms.push_back(new NormalDeviate(pf));
+ } else if (tname=="range") {
+ transforms.push_back(new Range(pf,*this));
+ } else {
+ std::cout << "Unknown transform : " << tname << '\n';
+ exit(1);
+ }
+ }
+}
+
+std::vector* Genotype_Phenotype_map::make_checkpoint_data() {
+ // Default is to store nothing
+ // All necessary information is normally taken from the parameter file.
+ return NULL;
+}
+
+void Genotype_Phenotype_map::resume_from_checkpoint_data(std::vector *data){
+ // default is to do nothing
+ // All necessary information is normally taken from the parameter file.
+}
+
+
+Diallelic_GP_map::Diallelic_GP_map(Trait& theT, Diallelic& theG, int loci_per_dim, ParameterFile& pf) : Genotype_Phenotype_map(theT), G(theG) {
+ this->loci_per_dim = loci_per_dim;
+ start_locus = G.add_loci(loci_per_dim * theT.get_dims());
+ read_all_transforms(pf);
+}
+
+Diallelic_GP_map::~Diallelic_GP_map() {
+}
+
+void Diallelic_GP_map::generate_phenotypes() {
+ T.assign(T.get_dims(), pop_size(), 0);
+ for (int ind=0; indloci_per_dim) {
+ std::cout << "Illegal initial trait value " << X_init << " of trait " << T.getName() << ".\n";
+ std::cout << "The allowed range is [-L, L], where L is the number of loci per dimension.\n";
+ std::cout << "Note that the initial value is an untransformed value.\n";
+ exit(1);
+ }
+ // set the first loci to +1, the rest remain -1
+ for (int ind=0; indloci_per_dim = loci_per_dim;
+ start_locus = G.add_loci(loci_per_dim * theT.get_dims());
+ read_all_transforms(pf);
+}
+
+ContinuousAlleles_GP_map::~ContinuousAlleles_GP_map() {
+}
+
+void ContinuousAlleles_GP_map::generate_phenotypes() {
+ T.assign(T.get_dims(), pop_size(), 0);
+ for (int ind=0; indX_init = X_init;
+}
+
+std::vector* Omnigenic_GP_map::make_checkpoint_data() {
+ // return a copy of weights:
+ return new std::vector(weights);
+}
+
+void Omnigenic_GP_map::resume_from_checkpoint_data(std::vector* data) {
+ weights = *data;
+}
+
+/*
+void Omnigenic_GP_map::add_to_checkpoint(Checkpoint &cp) {
+ // First add X_init as a vector with a single element:
+ std::vector X_init_vector(1, X_init);
+ cp.addData(new XData(X_init_vector));
+ cp.addData(new XData(weights));
+}
+
+void Omnigenic_GP_map::read_to_checkpoint(Checkpoint &cp, iSimfile &isf) {
+ cp.addData(new XData(isf,1)); // Read X_init parameter
+ cp.addData(new XData(isf,G.getLoci()*T.get_dims())); // Read weights
+}
+
+int Omnigenic_GP_map::resume_from_checkpoint(Checkpoint &cp, int dataIndex) {
+ X_init = (dynamic_cast&>(cp.getData(dataIndex++))).getData().at(0);
+ weights = (dynamic_cast&>(cp.getData(dataIndex++))).getData();
+ return dataIndex;
+}
+*/
+
diff --git a/TREES_code/genotype_phenotype_map.h b/TREES_code/genotype_phenotype_map.h
new file mode 100644
index 0000000..14fd909
--- /dev/null
+++ b/TREES_code/genotype_phenotype_map.h
@@ -0,0 +1,94 @@
+//
+// genotype_phenotype_map.hpp
+//
+// TREES -
+// A TRait-based Eco-Evolutionary Simulation tool
+// Copyright (C) 2017 Jörgen Ripa
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see .
+//
+// Contact: jorgen.ripa@biol.lu.se
+//
+
+#ifndef genotype_phenotype_map_h
+#define genotype_phenotype_map_h
+
+#include
+#include "types.h"
+#include "parameterFile.h"
+
+class Genetics;
+class Trait;
+class Transform;
+class Diallelic;
+class ContinuousAlleles;
+class Omnigenic;
+class Checkpoint;
+class iSimfile;
+class oSimfile;
+
+class Genotype_Phenotype_map {
+protected:
+ Trait& T;
+ Genotype_Phenotype_map(Trait& theTrait);
+ std::vector transforms;
+ void read_all_transforms(ParameterFile& pf);
+ double do_all_transforms(double y);
+public:
+ virtual ~Genotype_Phenotype_map();
+ virtual void generate_phenotypes()=0;
+ virtual void initialize(traitType X_init)=0;
+ virtual std::vector* make_checkpoint_data();
+ virtual void resume_from_checkpoint_data(std::vector* data);
+ int get_dims();
+ int pop_size();
+ Trait& get_trait();
+};
+
+class Diallelic_GP_map : public Genotype_Phenotype_map {
+ Diallelic& G;
+ int start_locus;
+ int loci_per_dim;
+public:
+ Diallelic_GP_map(Trait& theT, Diallelic& theG, int loci_per_dim, ParameterFile& pf);
+ virtual ~Diallelic_GP_map();
+ virtual void generate_phenotypes();
+ virtual void initialize(traitType X_init);
+ int get_loci_per_dim() { return loci_per_dim;}
+};
+
+class ContinuousAlleles_GP_map : public Genotype_Phenotype_map {
+ ContinuousAlleles& G;
+ int start_locus;
+ int loci_per_dim;
+public:
+ ContinuousAlleles_GP_map(Trait& theT, ContinuousAlleles& theG, int loci_per_dim, ParameterFile& pf);
+ virtual ~ContinuousAlleles_GP_map();
+ virtual void generate_phenotypes();
+ virtual void initialize(traitType X_init);
+};
+
+class Omnigenic_GP_map : public Genotype_Phenotype_map {
+ Omnigenic& G;
+ std::vector weights;
+ double X_init;
+public:
+ Omnigenic_GP_map(Trait& theT, Omnigenic& theG, int loci_per_dim, ParameterFile& pf);
+ virtual void generate_phenotypes();
+ virtual void initialize(traitType X_init);
+ virtual std::vector* make_checkpoint_data();
+ virtual void resume_from_checkpoint_data(std::vector* data);
+};
+
+#endif /* genotype_phenotype_map_hpp */
diff --git a/TREES_code/main.cpp b/TREES_code/main.cpp
index e4e2506..b9dc2be 100644
--- a/TREES_code/main.cpp
+++ b/TREES_code/main.cpp
@@ -38,18 +38,61 @@ int main(int argc, const char * argv[]) {
// Calculate table for interpolation of exponential function:
fillExpTable();// don't remove!!!!
+ int next_par_i = 1;
+ // check for command line options:
+ bool resume = false;
+ bool verbose = false;
+ while (next_par_i TREES parameterFile.txt [rep1 [rep2]]
std::string parameterFileName;
int rep1, rep2;
- if (argc<2) {
- std::cout << "Error: Missing parameter file. \n Syntax:\n Species parFile.txt [rep1] [rep2]\n";
- return 0;
+ if (argc-next_par_i < 1) {
+ std::cout << "Error: Missing parameter file. \n Syntax:\n TREES parameter_file [first_replicate [last_replicate]]\n";
+ return 1;
+ }
+ parameterFileName = argv[next_par_i++];
+ if (resume) {
+ // resume option, single replicate
+ // syntax: > TREES -resume parameter_file.txt results_file.txt [generation [new_results_file.txt]]
+ std::string resultsFile, newResultsFile;
+ time_type gen;
+ if (next_par_i=3) {
- rep1 = atoi(argv[2]);
- if(argc>=4) {
- rep2 = atoi(argv[3]);
+ // standard run option, one or more replicates
+ if(next_par_igetDims() != pref->getDims()) || (target->getDims() != strength->getDims()) ) {
+ if ((target->get_dims() != pref->get_dims()) || (target->get_dims() != strength->get_dims()) ) {
std::cout << "Target_selection error: Display, Preference and Strength traits must have the same number of dimensions\n";
exit(0);
}
@@ -63,13 +63,13 @@ Preference::Preference(Population& p, ParameterFile & pf) : pop(p)
dislimsq = disassortative_limit*disassortative_limit;
}
-double Preference::getPartnerWeight(int female, int male) {
+traitType Preference::getPartnerWeight(int female, int male) {
Fastexp fexp;
- double esum = 0;
- for (int d=0; dgetDims(); ++d) {
- double dx = getY(female,d) - getX(male,d);
- double cd = getC(female,d);
+ traitType esum = 0;
+ for (int d=0; dget_dims(); ++d) {
+ traitType dx = getY(female,d) - getX(male,d);
+ traitType cd = getC(female,d);
if (cd>=0) {
esum += cd * dx * dx;
} else {
diff --git a/TREES_code/mating.h b/TREES_code/mating.h
index 389705f..26f0e1e 100644
--- a/TREES_code/mating.h
+++ b/TREES_code/mating.h
@@ -28,6 +28,7 @@
#include "trait.h"
#include "parameterFile.h"
+#include "types.h" // defines traitType
class Population;
class Survival;
@@ -53,14 +54,14 @@ class Preference {
Trait* target;// X
Trait* pref; // Y
Trait* strength; // C
- double disassortative_limit;
- double dislimsq;
- double& getX(int individual,int dim) { return target->traitValue(individual,dim);}
- double& getY(int individual,int dim) { return pref->traitValue(individual,dim);}
- double& getC(int individual, int dim) { return strength->traitValue(individual,dim);}
+ traitType disassortative_limit;
+ traitType dislimsq;
+ traitType& getX(int individual,int dim) { return target->traitValue(individual,dim);}
+ traitType& getY(int individual,int dim) { return pref->traitValue(individual,dim);}
+ traitType& getC(int individual, int dim) { return strength->traitValue(individual,dim);}
public:
Preference(Population& p, ParameterFile & pf);
- double getPartnerWeight(int female, int male);
+ traitType getPartnerWeight(int female, int male);
};
diff --git a/TREES_code/matrix.h b/TREES_code/matrix.h
new file mode 100644
index 0000000..22df816
--- /dev/null
+++ b/TREES_code/matrix.h
@@ -0,0 +1,194 @@
+//
+// TREES -
+// A TRait-based Eco-Evolutionary Simulation tool
+// Copyright (C) 2017 Jörgen Ripa
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see .
+//
+// Contact: jorgen.ripa@biol.lu.se
+//
+
+#ifndef matrix_h
+#define matrix_h
+
+#include
+#include
+#include "simfile.h"
+
+// Matrix class without default constructing
+template class Matrix {
+protected:
+ T* X;
+ uint32_t M, N, X_size;
+ void initialize() {M=0; N=0; X=NULL; X_size=0;}
+public:
+ Matrix() {initialize();}
+
+ void reserve( uint32_t Mres, uint32_t Nres) {
+ if (Mres*Nres > X_size) {
+ T* newX = new T[Mres*Nres];
+ if (X && M*N>0) {
+ std::memcpy(newX, X, M*N*sizeof(T));
+ }
+ if (X) {
+ delete [] X;
+ }
+ X = newX;
+ X_size = Mres*Nres;
+ }
+ }
+
+ void assign(uint32_t newM, uint32_t newN, T val) {
+ reserve(newM,newN);
+ M = newM;
+ N = newN;
+ T* i = X;
+ for (int c=0; c& M2) {
+ initialize();
+ *this = M2;
+ }
+
+ // assignment operator:
+ Matrix& operator = (const Matrix &M2)
+ {
+ resize(M2.M, M2.N);
+ std::memcpy(X, M2.X, M2.M*M2.N*sizeof(T));
+ return *this;
+ }
+
+ Matrix(iSimfile& isf) {
+ initialize();
+ M = isf.read();
+ N = isf.read();
+ reserve(M,N);
+ isf.readCArray(X, M*N);
+ }
+ ~Matrix() { if (X) delete [] X; }
+
+ void write_to_file(oSimfile& osf) {
+ osf.write(M);
+ osf.write(N);
+ osf.writeCArray(X, M*N);
+ }
+ void read_from_file(iSimfile & isf) {
+ M = isf.read();
+ N = isf.read();
+ reserve(M,N);
+ isf.readCArray(X, M*N);
+ }
+ T& operator() (uint32_t row, uint32_t col) {
+ if (row0 && !my_isspace(instr[pos]) && !(instr[pos]==':') && bufpos<99) {
- buffer[bufpos++] = instr[pos++];
- }
- buffer[bufpos] = 0;
- if (bufpos>=99) {
- std::cout << "Buffer overflow : " << buffer << '\n';
- exit(0);
+std::string ParameterFile::get_next_value_string() {
+ std::string value_string("");
+ char c = skipWhiteSpace();
+ while (!eof() && !my_isspace(c) && c!=',' ) {
+ value_string.append(1,c);
+ c = pfile.get();
}
- S = buffer;
- return pos;
-}
-
-void stringToLower( std::string& s ) {
- std::transform(s.begin(), s.end(), s.begin(), ::tolower);
+ return value_string;
}
-void ParameterFile::getStringPair(std::string& s1, std::string& s2) {
- char line[500];
- getNextLine(line);
- int pos = 0;
- pos = ::skipWhiteSpace(line,pos);
- if (line[pos]==0) {
- s1 = "";
- s2 = "";
- } else {
- pos = ::getString(line,pos,s1);
- pos = ::skipWhiteSpace(line,pos);
- if (line[pos]==':') {
- ++pos;
- } else {
- std::cout << "Bad input : " << line << '\n';
- exit(0);
- }
- pos = ::skipWhiteSpace(line,pos);
- pos = ::getString(line,pos,s2);
- }
- stringToLower(s1);
- //stringToLower(s2);
+std::string ParameterFile::get_next_value_string_lower() {
+ std::string low_s = get_next_value_string();
+ std::transform(low_s.begin(), low_s.end(), low_s.begin(), ::tolower);
+ return low_s;
}
double ParameterFile::getDouble(const char* name) {
- std::string s1, s2;
- getStringPair(s1,s2);
- if (s1 != name ) {
+ if (get_next_name() != name ) {
std::cout << "Parameter error : \n";
std::cout << "Reading parameter : " << name << '\n';
- std::cout << "Found parameter : " << s1 << '\n';
+ std::cout << "Found parameter : " << get_next_name() << '\n';
exit(1);
}
+ std::string s2 = get_next_value_string();
return std::atof(s2.c_str());
}
@@ -126,7 +120,7 @@ int64_t ParameterFile::getPositiveLong(const char* name) {
int64_t i64 = getLong(name);
if (i64<=0) {
std::cout << "Parameter " << name << " has to be positive. Found value : " << i64 << '\n';
- exit(0);
+ exit(1);
}
return i64;
}
@@ -135,7 +129,7 @@ double ParameterFile::getPositiveDouble( const char* name) {
double d = getDouble(name);
if (d<=0) {
std::cout << "Parameter " << name << " has to be positive. Found value : " << d << '\n';
- exit(0);
+ exit(1);
}
return d;
@@ -145,7 +139,7 @@ int ParameterFile::getPositiveInt(const char* name) {
int i = (int)getLong(name);
if (i<=0) {
std::cout << "Parameter " << name << " has to be positive. Found value : " << i << '\n';
- exit(0);
+ exit(1);
}
return i;
}
@@ -168,20 +162,20 @@ bool ParameterFile::getBool(const char* name) {
std::string ParameterFile::getString(const char* name) {
- std::string s1, s2;
- getStringPair(s1,s2);
- if (s1!=name ) {
+ if (get_next_name() != name ) {
std::cout << "Parameter error : \n";
std::cout << "Reading parameter : " << name << '\n';
- std::cout << "Found parameter : " << s1 << '\n';
+ std::cout << "Found parameter : " << get_next_name() << '\n';
exit(1);
}
+ std::string s2 = get_next_value_string();
return s2;
}
std::string ParameterFile::getStringLower(const char* name) {
std::string s(getString(name));
- stringToLower(s);
+ std::transform(s.begin(), s.end(), s.begin(), ::tolower);
+ //stringToLower(s);
return s;
}
@@ -200,27 +194,78 @@ void ParameterFile::skipRestOfLine() {
}
}
-void ParameterFile::getNextLine(char *buffer) {
- int lastpos;
- // Skip comment lines starting with #, and blank lines
- buffer[0] = skipWhiteSpace();
- while (buffer[0]=='#') {
- skipRestOfLine();
- buffer[0] = skipWhiteSpace();
- }
- if (pfile.eof()) {
- buffer[0] = 0;
- return;
- }
- lastpos = 0;
- // read until comma, #, linefeed, or EOF:
- while (!(pfile.eof() || buffer[lastpos]=='\n' || buffer[lastpos]=='\r' || buffer[lastpos]==',' || buffer[lastpos]=='#' )) {
- char lastc = pfile.get();
- buffer[++lastpos] = lastc;
- }
- // skip comment at end of line:
- if (buffer[lastpos]=='#') {
- skipRestOfLine();
- }
- buffer[lastpos] = 0;
-}
+//void ParameterFile::getNextLine(char *buffer) {
+// int lastpos;
+// // Skip comment lines starting with #, and blank lines
+// buffer[0] = skipWhiteSpace();
+// while (buffer[0]=='#') {
+// skipRestOfLine();
+// buffer[0] = skipWhiteSpace();
+// }
+// if (pfile.eof()) {
+// buffer[0] = 0;
+// return;
+// }
+// lastpos = 0;
+// // read until comma, #, linefeed, or EOF:
+// while (!(pfile.eof() || buffer[lastpos]=='\n' || buffer[lastpos]=='\r' || buffer[lastpos]==',' || buffer[lastpos]=='#' )) {
+// char lastc = pfile.get();
+// buffer[++lastpos] = lastc;
+// }
+// // skip comment at end of line:
+// if (buffer[lastpos]=='#') {
+// skipRestOfLine();
+// }
+// buffer[lastpos] = 0;
+//}
+//int skipWhiteSpace(const char* str, int pos) {
+// while (str[pos]!=0 && my_isspace(str[pos])) {
+// ++pos;
+// }
+// return pos;
+//}
+//
+//int getString(const char* instr, int pos, std::string& S) {
+// char buffer[100];
+// buffer[0] = 0;
+// int bufpos = 0;
+// while ( instr[pos]>0 && !my_isspace(instr[pos]) && !(instr[pos]==':') && bufpos<99) {
+// buffer[bufpos++] = instr[pos++];
+// }
+// buffer[bufpos] = 0;
+// if (bufpos>=99) {
+// std::cout << "Buffer overflow : " << buffer << '\n';
+// exit(0);
+// }
+// S = buffer;
+// return pos;
+//}
+
+//void stringToLower( std::string& s ) {
+// std::transform(s.begin(), s.end(), s.begin(), ::tolower);
+//}
+
+//void ParameterFile::getStringPair(std::string& s1, std::string& s2) {
+// char line[500];
+// getNextLine(line);
+// int pos = 0;
+// pos = ::skipWhiteSpace(line,pos);
+// if (line[pos]==0) {
+// s1 = "";
+// s2 = "";
+// } else {
+// pos = ::getString(line,pos,s1);
+// pos = ::skipWhiteSpace(line,pos);
+// if (line[pos]==':') {
+// ++pos;
+// } else {
+// std::cout << "Bad input : " << line << '\n';
+// exit(0);
+// }
+// pos = ::skipWhiteSpace(line,pos);
+// pos = ::getString(line,pos,s2);
+// }
+// stringToLower(s1);
+// //stringToLower(s2);
+//}
+
diff --git a/TREES_code/parameterFile.h b/TREES_code/parameterFile.h
index d0321f9..75c43df 100644
--- a/TREES_code/parameterFile.h
+++ b/TREES_code/parameterFile.h
@@ -33,6 +33,9 @@ class ParameterFile {
public:
ParameterFile(const std::string& fileName);
~ParameterFile();
+ std::string get_next_name();
+ std::string get_next_value_string();
+ std::string get_next_value_string_lower();
double getDouble(const char* name);
double getPositiveDouble(const char* name);
int64_t getLong(const char* name);
@@ -42,14 +45,16 @@ class ParameterFile {
bool getBool(const char* name);
std::string getString(const char* name);
std::string getStringLower(const char* name);
- void getStringPair(std::string& s1, std::string& s2);
+ //void getStringPair(std::string& s1, std::string& s2);
void close();
bool eof() { return pfile.eof(); }
protected:
std::ifstream pfile;
+ std::string next_name;
+ std::ios::pos_type last_read_position;
char skipWhiteSpace();
void skipRestOfLine();
- void getNextLine(char* buffer);
+ //void getNextLine(char* buffer);
};
diff --git a/TREES_code/population.cpp b/TREES_code/population.cpp
index 77b50ed..48f47cc 100644
--- a/TREES_code/population.cpp
+++ b/TREES_code/population.cpp
@@ -28,14 +28,11 @@
#include "population.h"
#include "mating.h"
-#include "sample.h"
+#include "genotype_phenotype_map.h"
#include "randgen.h"
#include "fastmath.h"
Population::Population(ParameterFile& pf) {
- std::string modType;
- std::string modName;
-
// Read general parameters:
gene_tracking = pf.getBool("gene_tracking");
gene_sampling = pf.getBool("gene_sampling");
@@ -54,29 +51,27 @@ Population::Population(ParameterFile& pf) {
readGenetics(pf);
// Read traits:
- readTraits(pf, modType, modName);
+ readTraits(pf);
// read Space module, if specified:
- if (modType=="space") {
- addSpace(modName,pf);
- pf.getStringPair( modType, modName );
+ if (pf.get_next_name()=="space") {
+ addSpace(pf);
} else {
- space = new NullSpace(*this);
+ space = new Null_space(*this);
}
// Read fitness modules
fitnessList.clear();
- while (modType=="fitness") {
- addFitness(modName,pf);
- pf.getStringPair( modType, modName );
+ while (pf.get_next_name()=="fitness") {
+ addFitness(pf);
}
// Mating
- readMating(pf, modType, modName);
+ readMating(pf);
- if (modType.length()>0) {
- std::cout << "Unexpected module at end : " << modType << "\n";
+ if (pf.get_next_name().length()>0) {
+ std::cout << "Unexpected module at end : " << pf.get_next_name() << "\n";
exit(0);
}
@@ -96,35 +91,32 @@ void Population::readGenetics( ParameterFile& pf) {
genetics = new ContinuousAlleles(*this, pf);
} else if ( geneticsModel == "diallelic" ) {
genetics = new Diallelic(*this, pf);
+ } else if ( geneticsModel == "omnigenic" ) {
+ genetics = new Omnigenic(*this, pf);
} else {
std::cout << "Unknown Genetics model!\n";
exit(0);
}
}
-void Population::readTraits( ParameterFile& pf, std::string& modType, std::string& modName) {
+void Population::readTraits( ParameterFile& pf) {
traitList.clear();
- pf.getStringPair(modType, modName );
- while (modType.substr(0,5)=="trait") {
- if (modType=="trait") {
- addTrait(modName, pf);
- } else if (modType =="trait_constant") {
- addTraitConstant(modName, pf);
+ while (pf.get_next_name().substr(0,5)=="trait") {
+ if (pf.get_next_name()=="trait") {
+ addTrait(pf.get_next_value_string(), pf);
+ } else if (pf.get_next_name() =="trait_constant") {
+ addTraitConstant(pf.get_next_value_string(), pf);
} else {
- std::cout << "Unknown Trait type:" << modType << "\n";
- exit(0);
- }
- pf.getStringPair( modType, modName );
- while (modType=="transform") {
- addTransformToLastTrait(modName,pf);
- pf.getStringPair( modType, modName );
+ std::cout << "Unknown Trait type:" << pf.get_next_name() << "\n";
+ exit(1);
}
}
}
-void Population::readMating(ParameterFile& pf, std::string& modType, std::string& modName) {
- if (modType=="mating_pool") {
- switch (std::tolower(modName[0])) {
+void Population::readMating(ParameterFile& pf) {
+ if (pf.get_next_name()=="mating_pool") {
+ std::string mod_name = pf.get_next_value_string();
+ switch (std::tolower(mod_name[0])) {
case 's':
theMatingType = Selfing;
break;
@@ -136,23 +128,21 @@ void Population::readMating(ParameterFile& pf, std::string& modType, std::string
mate_s_space = pf.getDouble("s_space");
break;
default:
- std::cout << "Unknown Mating_pool type : " << modName << "\n";
+ std::cout << "Unknown Mating_pool type : " << mod_name << "\n";
exit(0);
break;
}
}else {
- std::cout << "Expected Mating_pool, found : " << modType << "\n";
- exit(0);
+ std::cout << "Expected Mating_pool, found : " << pf.get_next_name() << "\n";
+ exit(1);
}
mating_trials = pf.getPositiveInt("mating_trials");
// Read preference modules
matingPreferenceList.clear();
- pf.getStringPair( modType, modName );
- while (modType=="mating_preference") {
- addPreference(modName,pf);
- pf.getStringPair( modType, modName );
+ while (pf.get_next_name()=="mating_preference") {
+ addPreference(pf);
}
}
@@ -166,30 +156,13 @@ void Population::kill(int individual) {
somebodyDied = true;
}
-void Population::addTrait(std::string& traitName, ParameterFile& pf){
+void Population::addTrait(std::string traitName, ParameterFile& pf){
traitList.push_back(new Trait(traitName, *this, pf));
}
-void Population::addTraitConstant(std::string& traitName, ParameterFile& pf){
+void Population::addTraitConstant(std::string traitName, ParameterFile& pf){
traitList.push_back(new TraitConstant(traitName, *this, pf));
}
-
-void Population::addTransformToLastTrait(std::string& tname, ParameterFile& pf) {
- if (tname=="linear") {
- traitList.back()->addTransform(new LinearTransform(pf));
- } else if (tname=="abs") {
- traitList.back()->addTransform(new AbsTransform(pf));
- } else if (tname=="logistic") {
- traitList.back()->addTransform(new LogisticTransform(pf));
- } else if (tname=="normal_deviate") {
- traitList.back()->addTransform(new NormalDeviate(pf));
- } else {
- std::cout << "Unknown transform : " << tname << '\n';
- exit(0);
- }
-}
-
-
Trait* Population::findTrait(std::string& tname) {
for (Trait* t :traitList) {
if (t->getName() == tname) {
@@ -200,52 +173,46 @@ Trait* Population::findTrait(std::string& tname) {
exit(0);
}
-
-void Population::addFitness(std::string modName, ParameterFile &pf) {
- //Convert to lower case:
- std::transform(modName.begin(), modName.end(), modName.begin(), ::tolower);
-
- if (modName=="stabilizing_selection") {
+void Population::addFitness(ParameterFile &pf) {
+ std::string mod_name = pf.get_next_value_string_lower();
+ if (mod_name=="stabilizing_selection") {
fitnessList.push_back(new StabilizingSelection(*this, pf));
- } else if (modName=="density_dependence") {
+ } else if (mod_name=="density_dependence") {
fitnessList.push_back(new DensityDependence(*this, pf));
- } else if (modName == "resource_landscape") {
+ } else if (mod_name == "resource_landscape") {
fitnessList.push_back(new ResourceLandscape(*this, pf));
- } else if (modName=="discrete_resources") {
+ } else if (mod_name=="discrete_resources") {
fitnessList.push_back(new DiscreteResources(*this, pf));
- } else if (modName=="spatial_gradient") {
+ } else if (mod_name=="spatial_gradient") {
fitnessList.push_back(new SpatialGradient(*this, pf));
- } else if (modName=="catastrophes") {
+ } else if (mod_name=="catastrophes") {
fitnessList.push_back(new Catastrophe(*this, pf));
-// } else if (modName=="FinchModel") {
+// } else if (mod_name=="FinchModel") {
// fitnessList.push_back(new FinchModel(*this, pf));
} else {
- std::cout << "Unknown Fitness model : " << modName << "\n";
+ std::cout << "Unknown Fitness model : " << mod_name << "\n";
exit(0);
}
}
-void Population::addPreference(std::string modName, ParameterFile &pf) {
- //Convert to lower case:
- std::transform(modName.begin(), modName.end(), modName.begin(), ::tolower);
-
- if ( modName == "target_selection") {
+void Population::addPreference(ParameterFile &pf) {
+ std::string mod_name = pf.get_next_value_string_lower();
+ if ( mod_name == "target_selection") {
matingPreferenceList.push_back(new Preference(*this, pf));
} else {
- std::cout << "Unknown Mating_preference model : " << modName << "\n";
- exit(0);
+ std::cout << "Unknown Mating_preference model : " << mod_name << "\n";
+ exit(1);
}
}
-void Population::addSpace(std::string modName, ParameterFile &pf) {
- //Convert to lower case:
- std::transform(modName.begin(), modName.end(), modName.begin(), ::tolower);
- if ( modName == "none") {
- space = new NullSpace(*this);
- } else if ( modName == "discrete") {
- space = new DiscreteSpaceImp(*this, pf);
- } else if ( modName == "continuous") {
- space = new ContinuousSpace(*this, pf);
+void Population::addSpace(ParameterFile &pf) {
+ std::string mod_name = pf.get_next_value_string_lower();
+ if ( mod_name == "none") {
+ space = new Null_space(*this);
+ } else if ( mod_name == "discrete") {
+ space = new Discrete_space_imp(*this, pf);
+ } else if ( mod_name == "continuous") {
+ space = new Continuous_space(*this, pf);
} else {
std::cout << "Unknown Space model!\n";
exit(0);
@@ -302,7 +269,7 @@ std::vector Population::findDads() {
dads[mom] = mom;
}
} else if (withinPatchMating) {
- DiscreteSpace& theSpace = dynamic_cast(getSpace());
+ Discrete_space& theSpace = dynamic_cast(getSpace());
std::vector indList;
for (int patch=0; patch0) {
@@ -414,7 +381,6 @@ std::vector Population::findDads() {
++trialCount; // rejection, try another mate
}
}
-
}
}
return dads;
@@ -454,6 +420,7 @@ void Population::compactData(){
}
}
n = iw;
+ fitness.resize(n);
alive.assign(n, true);
somebodyDied = false;
}
@@ -474,38 +441,321 @@ Population::~Population() {
}
void Population::initialize() {
+ n = n0;
// initialize genomes:
- genetics->initialize(n0);
- // set genes to initial values:
+ genetics->initialize();// genome pre-allocation
+
for( Trait* tr : traitList) {
- tr->initialize(n0);
+ tr->initialize();// Normally sets genes to match initial trait values
}
if (gene_tracking) {
// create gene tables and assign initial values:
- genetics->initializeTracking(n0);
+ genetics->initializeTracking();
}
// assign initial position:
- space->initialize(n0);
- n = n0;
+ space->initialize();
alive.assign(n, true);
generatePhenotypes();
age = 0;
}
-Sample* Population::makeSample() {
- Sample* theSample = new Sample(age,n);
+int Population::calc_total_data_dimensions() {
+ int totd = 0;
+ for (Trait* tr : traitList) {
+ if (!tr->is_constant()) {
+ totd += tr->get_dims();
+ }
+ }
+ totd += space->getDims();
+ return totd;
+}
+
+void Population::resumeAtCheckpoint(Population_checkpoint& cp) {
+ n = cp.get_pop_size();
+ age = cp.get_generation();
+ genetics->resumeFromCheckpoint(cp.genetics, cp.geneListsCopy);
+ space->resumeFromCheckpoint(cp.space);
+ if(cp.GP_map_data.size()>0) {
+ int gp_map_index = 0;
+ // This only restores parameters (if any).
+ // Actual phenotypes are regenerated directly from genotypes (below).
+ for (Trait* t : traitList) {
+ if (!t->is_constant()) {
+ t->get_GP_map()->resume_from_checkpoint_data(cp.get_map_data(gp_map_index));
+ gp_map_index++;
+ }
+ }
+ }
+ generatePhenotypes();
+ alive.assign(n, true);
+ somebodyDied = false;
+ fitness.assign(n, 1.0);
+}
+
+Sample_base::~Sample_base(){
+}
+
+/*class Population_sample {
+ protected:
+ time_type generation;
+ double cputime; // in seconds
+ int pop_size;
+ Genetics_sample* genes_sample;
+ std::vector traits;
+ Space_sample* space_sample;
+ public:
+ Population_sample(Population& pop, double cputime);
+ void write_to_file(oSimfile& osf);
+ };
+*/
+
+Population_sample::Population_sample(Population& pop, double cputime) :
+ Sample_base(pop.getAge(), cputime, pop.size()) {
+
+ if (pop.gene_sampling) {
+ genes_sample = pop.getGenetics().make_sample();
+ } else {
+ genes_sample = NULL;
+ }
+ traits.clear();
+ for(Trait* tr : pop.traitList) {
+ if (!tr->is_constant()) {
+ traits.push_back(Trait_sample(*tr));
+ }
+ }
+ space_sample = pop.getSpace().make_sample();
+}
+
+Population_sample::~Population_sample() {
+ if (genes_sample) {
+ delete genes_sample;
+ }
+ delete space_sample;
+}
+
+void Population_sample::write_to_file(oSimfile &osf) {
+ osf.write(generation);
+ osf.write(cputime);
+ osf.write(pop_size);
+
+ if(genes_sample) {
+ osf.write((char)true); // sizeof(bool) is implementation-specific, but sizeof(char) is always 1
+ genes_sample->write_to_file(osf);
+ } else {
+ osf.write((char)false);
+ }
+
+ osf.write((int)traits.size()); // The number of sampled traits (not constants)
+ for(Trait_sample& trs : traits) {
+ trs.write_to_file(osf);
+ }
+
+ space_sample->write_to_file(osf);
+}
+
+Population_sample::Population_sample(iSimfile& isf) {
+ generation = isf.read();
+ cputime = isf.read();
+ pop_size = isf.read();
+ bool gene_sampling = (bool)isf.read();
if (gene_sampling) {
- genetics->addToSample(*theSample);
+ genes_sample = Genetics_sample::create_from_file(isf);
+ } else {
+ genes_sample = NULL; // Important!!!
+ }
+ int ntraits = isf.read();
+ traits.clear();
+ for (int ti=0; tiisConstant()) {
- tr->addToSample(*theSample);
+ space_sample = Space_sample::create_from_file(isf);
+}
+
+Population_checkpoint::Population_checkpoint(Population& pop, seed_type cpseed, double cputime) :
+Sample_base(pop.getAge(), cputime, pop.size()) {
+ seed = cpseed;
+ // copy all genomes, and possibly ID:s:
+ genetics = pop.getGenetics().make_sample();
+ // Store gene lists, if applicable:
+ if (pop.isTrackingGenes()) {
+ geneListsCopy = pop.getGenetics().getGeneLists(); // deep copy
+ } else {
+ geneListsCopy.clear();
+ }
+ GP_map_data.clear();
+ for(Trait* tr : pop.traitList) {
+ if (!tr->is_constant()) {
+ std::vector* data = tr->get_GP_map()->make_checkpoint_data();
+ if (data) {
+ GP_map_data.push_back(data);
+ }
}
}
- space->addToSample(*theSample);
- if (gene_tracking) {
- genetics->addGeneIDsToSample(*theSample);
+ space = pop.getSpace().make_sample();
+}
+
+Population_checkpoint::~Population_checkpoint() {
+ delete genetics;
+ delete space;
+}
+
+void Population_checkpoint::write_to_file(oSimfile &osf) {
+ osf.write(seed);
+ osf.write(generation);
+ osf.write(cputime);
+ osf.write(pop_size);
+ genetics->write_to_file(osf);
+ osf.write((size_type)geneListsCopy.size());
+ //One genelist per locus:
+ for (GeneList& list : geneListsCopy) {
+ list.writeGenes(osf);
+ }
+ osf.write((size_type)GP_map_data.size());
+ for (int ti=0; tiwrite_to_file(osf);
+}
+
+Population_checkpoint::Population_checkpoint(iSimfile& isf) {
+ seed = isf.read();
+ generation = isf.read();
+ cputime = isf.read();
+ pop_size = isf.read();
+ // Genetics:
+ genetics = Genetics_sample::create_from_file(isf);
+ size_type genelist_count = isf.read(); // number of loci, really
+ geneListsCopy.clear();
+ geneListsCopy.reserve(genelist_count);
+ for (int li=0; li();
+ GP_map_data.clear();
+ for (int mi=0; mi());
+ isf.readVector(*GP_map_data.back());
}
- return theSample;
+ // Space:
+ space = Space_sample::create_from_file(isf);
+}
+
+
+Microsample::Microsample(Population& pop, char option, double cputime) :
+Sample_base(pop.getAge(), cputime, pop.size()){
+ this->option = option;
+ means.clear();
+ covariances.clear();
+ switch (option) {
+ case 'n':
+ break;
+ case 'm':
+ calc_means(pop);
+ break;
+ case 'v':
+ calc_means(pop);
+ calc_variances(pop);
+ break;
+ case 'c':
+ calc_means(pop);
+ calc_covariances(pop);
+ break;
+ default:
+ std::cout << "Microsample error.'\n'";
+ exit(1);
+ break;
+ }
+}
+
+void Microsample::write_to_file(oSimfile &osf) {
+ osf.write(option);
+ osf.write(generation);
+ osf.write(cputime);
+ osf.write(pop_size);
+ osf.writeVector(means);
+ osf.writeVector(covariances);
}
+Microsample::Microsample(iSimfile& isf) {
+ option = isf.read();
+ generation = isf.read();
+ cputime = isf.read();
+ pop_size = isf.read();
+ isf.readVector(means);
+ isf.readVector(covariances);
+}
+
+void Microsample::calc_means(Population &pop) {
+ for (Trait* tr : pop.getTraits()) {
+ if (!tr->is_constant()) {
+ for (int dim=0; dimget_dims(); ++dim) {
+ means.push_back(tr->row_mean(dim));
+ }
+ }
+ }
+ // mean positions:
+ for (int dim=0; dimis_constant()) {
+ for (int dim=0; dimget_dims(); ++dim) {
+ covariances.push_back(tr->row_variance(dim,means.at(mean_i)));
+ ++mean_i;
+ }
+ }
+ }
+ // spatial variation:
+ for (int dim=0; dim bigM(pop_size,total_dims);
+ int Mdim=0;
+ for (Trait* tr : pop.getTraits()) {
+ if (!tr->is_constant()) {
+ for (int Tdim=0; Tdimget_dims(); ++Tdim) {
+ for (int ind=0; ind fitnessList;
MatingType theMatingType;
double mate_s_space;
+ int n0; // initial population size
+
+// Model structures:
+ Genetics* genetics;
+ Space* space;
std::vector matingPreferenceList;
- std::vector findDads();
std::vector traitList; // a list of traits
- Space* space;
+ std::vector fitnessList;
+
+// Dynamic parameters:
+ int n; //population size (sometimes including dead)
+ time_type age;
std::vector alive;
std::vector fitness;
bool somebodyDied;
+
+ std::vector findDads();
void allAlive();
void readGenetics( ParameterFile& pf);
- void readTraits( ParameterFile& pf, std::string& modType, std::string& modName);
- void addFitness(std::string modName, ParameterFile& pf);
- void readMating(ParameterFile& pf, std::string& modType, std::string& modName);
- void addPreference(std::string modName, ParameterFile& pf);
- void addSpace(std::string modName, ParameterFile& pf);
+ void readTraits( ParameterFile& pf);
+ void addFitness(ParameterFile& pf);
+ void readMating(ParameterFile& pf);
+ void addPreference(ParameterFile& pf);
+ void addSpace(ParameterFile& pf);
void compactData();
void generatePhenotypes();
- void addTrait(std::string& traitName, ParameterFile& pf);
- void addTraitConstant(std::string& traitName, ParameterFile& pf);
- void addTransformToLastTrait(std::string& tname, ParameterFile& pf);
- double& getTrait(int individual, int trait) { return traitList.at(trait)->traitValue(individual); }
+ void addTrait(std::string traitName, ParameterFile& pf);
+ void addTraitConstant(std::string traitName, ParameterFile& pf);
+ //void addTransformToLastTrait(std::string tname, ParameterFile& pf);
+ traitType& getTrait(int individual, int trait) { return traitList.at(trait)->traitValue(individual); }
std::vector& getTraits() {return traitList;}
+ int calc_total_data_dimensions();
+
+
public:
Population(ParameterFile& pf);
~Population();
- // Added as a quick fix, remove in Super4:
void initialize();
- Sample* makeSample();
+ void resumeAtCheckpoint(Population_checkpoint& cp);
+ // Running:
void makeNextGeneration(); //
void reproduce();
void survive();
void kill(int individual);
+ // Info:
int size() { return n;}
Genetics& getGenetics() {return *genetics;}
Space& getSpace() {return *space;}
int getF() { return F;}
Trait* findTrait(std::string& name);
- timeType getAge() { return age;}
+ time_type getAge() { return age;}
bool isTrackingGenes() { return gene_tracking; }
};
+class Sample_base {
+protected:
+ time_type generation;
+ double cputime; // in seconds
+ int pop_size;
+
+ void set(time_type gen, double cpu, int n)
+ { generation=gen; cputime=cpu; pop_size=n;}
+ Sample_base(time_type generation, double cputime, int n)
+ { set(generation,cputime,n);}
+ Sample_base()
+ {set(0,0,0);}
+
+public:
+ time_type get_generation() { return generation; }
+ double get_cputime() { return cputime;}
+ int get_pop_size() {return pop_size;}
+ virtual ~Sample_base();
+};
+
+
+class Population_sample : public Sample_base {
+ protected:
+ Genetics_sample* genes_sample;
+ std::vector traits;
+ Space_sample* space_sample;
+ public:
+ Population_sample(Population& pop, double cputime);
+ Population_sample(iSimfile& isf);
+ virtual ~Population_sample();
+ void write_to_file(oSimfile& osf);
+};
+
+class Population_checkpoint : public Sample_base {
+ friend class Population;
+ // stores genes, genelists and space, and possible trait parameters
+ // (trait values can be reconstructed from genes)
+ seed_type seed;
+ Genetics_sample* genetics;
+ std::vector geneListsCopy;
+ Space_sample* space;
+ // A vector of vectors (one per non-constant trait, or none):
+ std::vector*> GP_map_data; // possible trait parameters (see Omnigenic model)
+
+ public:
+ Population_checkpoint(Population& pop, seed_type cpseed, double cputime);
+ Population_checkpoint(iSimfile& isf);
+ ~Population_checkpoint();
+ void write_to_file(oSimfile& osf);
+ seed_type get_seed() { return seed;}
+ std::vector* get_map_data(int mi) { return GP_map_data.at(mi);}
+};
+
+class Microsample : public Sample_base {
+ // optionally stores means, variances and covariances of all traits and spatial positions
+ protected:
+ void calc_means(Population& pop);
+ void calc_variances(Population& pop);
+ void calc_covariances(Population& pop);
+
+ public:
+ char option;
+ std::vector means;
+ std::vector covariances; // either just variances or all covariances
+ Microsample(Population & pop, char option, double cputime);
+ Microsample(iSimfile& isf);
+ void write_to_file(oSimfile& osf);
+};
+
#endif /* defined(__Species__population__) */
diff --git a/TREES_code/randgen.cpp b/TREES_code/randgen.cpp
index 38617b0..2ab1ea0 100644
--- a/TREES_code/randgen.cpp
+++ b/TREES_code/randgen.cpp
@@ -27,37 +27,81 @@
#include
#include
#include
+#include
+
using namespace std;
#include "randgen.h"
+std::mt19937_64 the_generator; // Turns out the 64-bit version is faster than the 32-bit (Mac)
+std::uniform_real_distribution unif_dist(0.0,1.0);
+
int myround(double x) { return (int)floor(x+0.5);}
-unsigned randomize()
+seed_type randomize()
{
// set seed as microseconds since Epoch.
using namespace std::chrono;
system_clock::time_point now = system_clock::now();
- unsigned micros = (unsigned)duration_cast(now.time_since_epoch()).count();
+ seed_type micros = (seed_type)duration_cast(now.time_since_epoch()).count();
+
+ the_generator.seed(micros);
- srand(micros);
- std::cout << "Randomize seed : " << micros << '\n';
+ // old code:
+// srand(micros);
+// std::cout << "Randomize seed : " << micros << '\n';
+
return micros;
}
-void randseed(unsigned seed) {
- srand(seed);
- std::cout << "Randseed set : " << seed << '\n';
+void randseed(seed_type seed) {
+ // new code:
+ the_generator.seed(seed);
+// // old code:
+// srand(seed);
+// std::cout << "Randseed set : " << seed << '\n';
+}
+
+seed_type make_seed() {
+ return the_generator();
}
double rand1() // uniform [0,1)
{
- int r = rand();
- double C = int64_t(RAND_MAX)+1.0; // RAND_MAX is 2^31 in OS-X 10.9.5
- return r / C;//double(RAND_MAX+1);
+ // new code:
+ return unif_dist(the_generator);
+// // old code:
+// int r = rand();
+// double C = int64_t(RAND_MAX)+1.0; // RAND_MAX is 2^31 in OS-X 10.9.5
+// return r / C;//double(RAND_MAX+1);
+}
+
+bitGenerator::bitGenerator(){
+ regenerate();
+}
+
+void bitGenerator::regenerate() {
+ // new code:
+ bits = (uint64_t)the_generator();
+ bitsleft = 64;
+
+// // old code:
+// bits = (unsigned)(rand1()*65536);
+// bitsleft = 16;
}
+//bool bitGenerator::nextBit() {
+// //return rand1()<0.5;
+// if (bitsleft==0) {
+// regenerate();
+// }
+// bool reply = bits & 1;
+// bits >>= 1;
+// --bitsleft;
+// return reply;
+//}
+
double randn() {
static bool iset = false;
static double gset = 0;
@@ -81,15 +125,18 @@ double randn() {
// A double exponential distributed random number,
// i.e. exponentially distributed in both + and - direction
-// absmean is mean(abs(x))
-double doubleExp(double absmean) {
+double doubleExp(double variance) {
+ double a = sqrt(variance/2); // distribution parameter
double r = rand1(); // r is recycled
- if (r<0.5) {
- // now r is uniform [0, 0.5)
- return -absmean*log(2*r);
+ while (r==0) { // just in case, should be a rare event
+ r = rand1();
+ }
+ if (r<=0.5) {
+ // now r is uniform (0, 0.5]
+ return -a*log(2*r);
} else {
- // now r is uniform [0.5, 1)
- return absmean*log(2*r-1);
+ // now r is uniform (0.5, 1)
+ return a*log(2*r-1);
}
}
@@ -207,40 +254,4 @@ int weightedChoiceCumSum( vector& cumW ) {
return (int) (cwp - cumW.begin());
}
-// This is inline for speed:
-/* void bitGenerator::regenerate() {
- bits = (unsigned)rand1()*65536;
- bitsleft = 16;
-}*/
-
-bitGenerator::bitGenerator(){
- regenerate();
-}
-
-void bitGenerator::regenerate() {
- bits = (unsigned)(rand1()*65536);
- bitsleft = 16;
-}
-bool bitGenerator::nextBit() {
- //return rand1()<0.5;
- if (bitsleft==0) {
- regenerate();
- }
- bool reply = bits & 1;
- bits >>= 1;
- --bitsleft;
- return reply;
-}
-
-
-// This is inline for speed:
-/* bool bitGenerator::nextBit() {
- if (bitsleft==0) {
- regenerate();
- }
- bool reply = bits % 2;
- bits >>= 1;
- --bitsleft;
- return reply;
-} */
diff --git a/TREES_code/randgen.h b/TREES_code/randgen.h
index d1e59b2..2c2841e 100644
--- a/TREES_code/randgen.h
+++ b/TREES_code/randgen.h
@@ -25,29 +25,38 @@
#include
+#include "types.h"
+
+double rand1(); // uniform [0,1)
+double randn();
int poissrnd(double lambda);
int poissrndAppr(double lambda);
int binornd(int n, double p);
int binorndAppr(int n, double p);
int binoSmallP(int n, double p);
double doubleExp(double absmean);
-double randn();
-double rand1(); // uniform [0,1)
-unsigned randomize(); // also return the seed
-void randseed(unsigned seed);
+seed_type randomize(); // also return the seed
+void randseed(seed_type seed);
+seed_type make_seed();
int weightedChoice( std::vector& weights);
int weightedChoiceCumSum( std::vector& cumWeights);
// class to generate a stream of random bits (true/false)
class bitGenerator {
- unsigned bits;
+ uint64_t bits;
int bitsleft;
void regenerate();
public:
bitGenerator();
- bool nextBit();
-
+ inline bool nextBit(){
+ if (bitsleft-- == 0) {
+ regenerate();
+ }
+ bool reply = bits & 1;
+ bits >>= 1;
+ return reply;
+ }
};
#endif
diff --git a/TREES_code/simfile.cpp b/TREES_code/simfile.cpp
index f767f38..f2fcf11 100644
--- a/TREES_code/simfile.cpp
+++ b/TREES_code/simfile.cpp
@@ -21,23 +21,67 @@
#include
+#include
#include
#include "simfile.h"
-Simfile::Simfile(std::string filename) {
+void endian_test() {
// Test for Little Endian:
unsigned char SwapTest[2] = { 1, 0 };
short w = *(short *) SwapTest;
assert( w == 1 );
-
+}
+
+oSimfile::oSimfile(std::string filename) {
+ endian_test();
open(filename.c_str(), std::fstream::binary | std::fstream::out);
imbue(std::locale::classic());
+ write(current_file_version);
+}
+
+void oSimfile::close() {
+ std::ofstream::close();
}
-Simfile::~Simfile() {
+oSimfile::~oSimfile() {
close();
}
+void oSimfile::writeString(const std::string& s) {
+ std::ofstream::write(s.c_str(), s.length());
+ put(0);
+}
+
+iSimfile::iSimfile(std::string filename) {
+ endian_test();
+ open(filename.c_str(), std::fstream::binary | std::fstream::in);
+ if (!is_open()) {
+ std::cout << "Failed to open file " << filename << std::endl;
+ exit(1);
+ }
+ imbue(std::locale::classic());
+ int file_version = read();
+ if (file_version!=current_file_version ) {
+ std::cout << "File error. File " << filename << " is version " << file_version << " instead of " << current_file_version << std::endl;
+ exit(1);
+ }
+}
+
+iSimfile::~iSimfile() {
+ close();
+}
+
+std::string iSimfile::readString() {
+ std::stringstream buffer;
+ std::ifstream::get(*buffer.rdbuf(), char(0));
+ // get the null character:
+ std::ifstream::get();
+ return buffer.str();
+}
+
+
+
+/*
void Simfile::write(float x) {
std::ofstream::write((char*)&x, sizeof(float));
}
@@ -50,8 +94,15 @@ void Simfile::write(int64_t i) {
std::ofstream::write((char*)&i, sizeof(int64_t));
}
-void Simfile::write(idType id) {
- std::ofstream::write((char*)&id, sizeof(idType));
+void Simfile::write(id_type id) {
+ std::ofstream::write((char*)&id, sizeof(id_type));
+}
+*/
+
+
+/*
+void Simfile::writeArray( signed char* A, int size) {
+ std::ofstream::write((char*)A, sizeof(signed char)*size);
}
void Simfile::writeArray( float* A, int size) {
@@ -72,12 +123,7 @@ void Simfile::writeArray( int* A, int size) {
std::ofstream::write((char*)A, sizeof(int)*size);
}
-void Simfile::writeArray( idType* A, int size) {
- std::ofstream::write((char*)A, sizeof(idType)*size);
-}
-
-void Simfile::writeString(std::string s) {
- std::ofstream::write(s.c_str(), s.length());
- put('\n');
- put(0);
+void Simfile::writeArray( id_type* A, int size) {
+ std::ofstream::write((char*)A, sizeof(id_type)*size);
}
+*/
diff --git a/TREES_code/simfile.h b/TREES_code/simfile.h
index 7a96dcc..789e92c 100644
--- a/TREES_code/simfile.h
+++ b/TREES_code/simfile.h
@@ -24,27 +24,97 @@
#define __Species__simfile__
#include
+#include
#include
#include
#include "geneTracking.h"
+static const uint32_t current_file_version = 4; // as of TREES 1.3
+
// class for saving data in binary format
-// protected inheritance to prevent other output than
-// of supported types
-class Simfile : protected std::ofstream {
+// protected inheritance to prevent non-supported output
+class oSimfile : protected std::ofstream {
public:
- Simfile(std::string filename);
- ~Simfile();
- void write(float x);
+ oSimfile(std::string filename);
+ ~oSimfile();
+ template void write(const WriteType data) {
+ std::ofstream::write((char*)&data, sizeof(WriteType));
+
+ }
+ templatevoid writeArray( const WriteType * A, int size) {
+ std::ofstream::write((char*)A, sizeof(WriteType)*size);
+ }
+ templatevoid writeCArray( const WriteType * A, int size) {
+ std::ofstream::write((char*)A, sizeof(WriteType)*size);
+ }
+ // This also stores vector length:
+ templatevoid writeVector( const std::vector