From 70c021dfe740b33abccde1df09a3cd4e653d17f6 Mon Sep 17 00:00:00 2001 From: jorgenripa Date: Thu, 23 Jan 2020 11:57:08 +0100 Subject: [PATCH] Add files via upload --- TREES_code/bit_vector.cpp | 96 +++++ TREES_code/bit_vector.h | 47 ++ TREES_code/build.sh | 17 +- TREES_code/fitness.cpp | 32 +- TREES_code/fitness.h | 7 +- TREES_code/geneTracking.cpp | 67 ++- TREES_code/geneTracking.h | 44 +- TREES_code/genetics.cpp | 597 ++++++++++++++++---------- TREES_code/genetics.h | 169 ++++++-- TREES_code/genotype_phenotype_map.cpp | 254 +++++++++++ TREES_code/genotype_phenotype_map.h | 94 ++++ TREES_code/main.cpp | 83 +++- TREES_code/mating.cpp | 12 +- TREES_code/mating.h | 13 +- TREES_code/matrix.h | 194 +++++++++ TREES_code/parameterFile.cpp | 205 +++++---- TREES_code/parameterFile.h | 9 +- TREES_code/population.cpp | 460 +++++++++++++++----- TREES_code/population.h | 123 +++++- TREES_code/randgen.cpp | 117 ++--- TREES_code/randgen.h | 23 +- TREES_code/simfile.cpp | 72 +++- TREES_code/simfile.h | 90 +++- TREES_code/simulation.cpp | 185 ++++---- TREES_code/simulation.h | 23 +- TREES_code/space.cpp | 337 +++++++++------ TREES_code/space.h | 129 ++++-- TREES_code/trait.cpp | 143 +++--- TREES_code/trait.h | 76 ++-- TREES_code/transform.cpp | 58 ++- TREES_code/transform.h | 42 +- TREES_code/types.h | 9 +- 32 files changed, 2800 insertions(+), 1027 deletions(-) create mode 100644 TREES_code/bit_vector.cpp create mode 100644 TREES_code/bit_vector.h create mode 100644 TREES_code/genotype_phenotype_map.cpp create mode 100644 TREES_code/genotype_phenotype_map.h create mode 100644 TREES_code/matrix.h 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& v) { + write((size_type)v.size()); + std::ofstream::write((char*)&v[0], sizeof(WriteType)*v.size()); + } + + void writeString(const std::string& s); + + void close(); + +/* void write(float x); void write(int i); void write(int64_t i); - void write(idType id); + void write(id_type id); + */ +/* + void writeArray( signed char* A, int size); + void writeArray( int* A, int size); void writeArray( float* A, int size); void writeArray( double* A, int size); - void writeArray( int* A, int size); - void writeArray( idType* A, int size); - void writeString(std::string s); + void writeArray( id_type* A, int size); + */ +}; + +/* +template oSimfile& operator <<(oSimfile& osf, std::vector& v) { + osf.write((int)v.size()); + for (int i=0; i void read(ReadType& data) { +// std::ifstream::read((char*)&data, sizeof(ReadType)); +// +// } + template ReadType read() { + ReadType data; + std::ifstream::read((char*)&data, sizeof(ReadType)); + return data; + } + templatevoid readArray( std::vector& A, int size) { + A.assign(size,(ReadType)0); + std::ifstream::read((char*)&A[0], sizeof(ReadType)*size); + } + templatevoid readCArray( ReadType* Ap, int count) { + std::ifstream::read((char*)Ap, sizeof(ReadType)*count); + } + templatevoid readVector(std::vector& v) { + size_type vlength = read(); + readArray(v, vlength); + } + std::string readString(); }; +template iSimfile& operator >>(iSimfile& isf, std::vector& v) { + int count = isf.read(); + v.assign(count,T()); + for (int i=0; i> v[i]; // only works for some supported types + } +} + #endif /* defined(__Species__simfile__) */ diff --git a/TREES_code/simulation.cpp b/TREES_code/simulation.cpp index 9ff6acb..a2d0e46 100644 --- a/TREES_code/simulation.cpp +++ b/TREES_code/simulation.cpp @@ -32,27 +32,27 @@ #include "population.h" #include "randgen.h" #include "simfile.h" -#include "sample.h" #include "mytime.h" -Simulation::Simulation(std::string& parameterFileName){ +Simulation::Simulation(std::string& parameterFileName, bool verbose){ ParameterFile pf(parameterFileName); + this->verbose = verbose; simName = parameterFileName; // skip file extension: simName.erase(simName.find_last_of('.')); std::cout << "name: " << simName << '\n'; - verbose = pf.getBool("verbose"); tmax = pf.getPositiveLong("t_max"); sampleInterval = pf.getPositiveLong("sample_interval"); + std::string option_string = pf.getStringLower("microsamples"); + microsample_option = option_string[0]; checkpointInterval = pf.getLong("checkpoint_interval"); // checkPoints==0 turns this feature off keepOldCheckpoints = pf.getBool("keep_old_checkpoints"); std::string seedS = pf.getStringLower("seed"); if (seedS[0]=='r') { seed = randomize(); } else { - seed = (unsigned)std::atol(seedS.c_str()); - randseed(seed); + seed = (seed_type)std::atoll(seedS.c_str()); } //Parse population parameters (including gene tracking/sampling): @@ -71,18 +71,25 @@ Simulation::Simulation(std::string& parameterFileName){ Simulation::~Simulation() { //std::cout << "Killing simulation\n"; delete thePop; - for (int si=0; siinitialize(); - sampleList.clear(); + sample_list.clear(); + microsample_list.clear(); checkpointList.clear(); makeFileName(replicate); // Sample generation 0: - sampleList.push_back(thePop->makeSample()); + sample_list.push_back(new Population_sample(*thePop,0)); + if (microsample_option!='n') { + microsample_list.push_back(new Microsample(*thePop,microsample_option,0)); + } } void Simulation::makeFileName(int replicate) { @@ -91,35 +98,49 @@ void Simulation::makeFileName(int replicate) { resultsFileName = fileName.str(); } -void Simulation::runFromGeneration(timeType gen) { - double t1 = getNow(); - double prevSec = 0; - - for(timeType t=gen+1; t<=tmax; ++t) { +void Simulation::runFromGeneration(time_type gen, seed_type gen_seed, double time_so_far) { + double start_time = getNow(); + double prevSec = time_so_far; + randseed(gen_seed); + for(time_type t=gen+1; t<=tmax; ++t) { + if (verbose) { + std::cout << t << std::flush; + } //std::cout << "make gen " << t << '\n'; thePop->makeNextGeneration(); + // Seed for next generation (may be saved in checkpoint): + seed_type next_seed = make_seed(); + randseed(next_seed); + + double cputime = getNow()-start_time+time_so_far; + + if (verbose) { + int backsteps = int(log10(t))+1; + std::cout << std::string(backsteps,'\b'); + } if (thePop->size()==0) { std::cout<< "Extinction at t = " << t << "!\n"; break; } + if (microsample_option!='n') { + microsample_list.push_back(new Microsample(*thePop, microsample_option, cputime)); + } if( t % sampleInterval == 0 ) { - sampleList.push_back(thePop->makeSample()); + sample_list.push_back(new Population_sample(*thePop,cputime)); if (verbose) { std::cout << "Generation " << t << ", population size " << thePop->size() << '\n'; - double nowSec = getNow()-t1; - int timeLeft = int( (nowSec-prevSec)*(tmax-t)/sampleInterval ); - std::cout << "Time " << time2str(int(nowSec)) << ", time left : " << time2str(timeLeft) << "\n"; - prevSec = nowSec; + int timeLeft = int( (cputime-prevSec)*(tmax-t)/sampleInterval ); + std::cout << "Time " << time2str(int(cputime)) << ", time left : " << time2str(timeLeft) << "\n"; + prevSec = cputime; } } if (checkpointInterval>0 && t%checkpointInterval==0) { - unsigned seed = rand(); + Population_checkpoint* cp = new Population_checkpoint(*thePop,next_seed,cputime); if (keepOldCheckpoints) { - checkpointList.push_back(thePop->makeCheckpoint(seed)); + checkpointList.push_back(cp); } else { - checkpointList.assign(1,thePop->makeCheckpoint(seed)); + checkpointList.assign(1,cp); } - randseed(seed); save(); } } @@ -128,98 +149,114 @@ void Simulation::runFromGeneration(timeType gen) { } } +void truncate_sample_list(std::vector& list, time_type last_gen) { + if (list.size()>0) { + int si=0; + time_type gen = list.at(si)->get_generation(); + while (genget_generation(); + } + if (gen>last_gen) { + --si; + } + // si is now index of last generation, si+1 is size of valid list + if (list.size() > si+1) { + for (int del_i=si+1; del_igetGeneration(); + if (resumeGen==-1) { // use last checkpoint: + resumeGen = checkpointList.back()->get_generation(); } else { // Find the right checkpoint (they may not be saved at equal intervals): - int cpi = 0; - timeType gen = checkpointList[cpi]->getGeneration(); - while (gengetGeneration(); - } - if (gen!=resumeGen) { + truncate_sample_list(checkpointList, resumeGen); + if (checkpointList.size()==0 || checkpointList.back()->get_generation() !=resumeGen) { std::cout << "Error. Can't find checkpoint for generation " << resumeGen << " in file " << resultsFile << '\n'; exit(0); } - // truncate lists: - checkpointList.resize(cpi+1); - } - int si = 0; - timeType gen = sampleList[si]->getGeneration(); - while (gengetGeneration(); } - if (gen>resumeGen) { - --si; - } - sampleList.resize(si+1); - thePop->resumeFromCheckpoint(*checkpointList.back()); - randseed(checkpointList.back()->seed); + // truncate microsample list: + truncate_sample_list(microsample_list, resumeGen); + // truncate sample list: + truncate_sample_list(sample_list, resumeGen); + + Population_checkpoint& cp = *dynamic_cast(checkpointList.back()); + thePop->resumeAtCheckpoint(cp); resultsFileName = newResultsFile; - runFromGeneration(resumeGen); + runFromGeneration(resumeGen, cp.get_seed(), cp.get_cputime()); } void Simulation::save() { - int fileVersion = 2; // As of v. 1.1 - std::cout << "Saving " << resultsFileName << '\n'; + std::cout << "Saving " << resultsFileName << "..." << std::flush; oSimfile simfile(resultsFileName); - simfile.write(fileVersion); - simfile.write((int)seed); + simfile.write(seed); simfile.writeString(parameterFileCopy.str()); - simfile.write(thePop->getGenetics().getLoci()); - // save all samples: - simfile.write((int)sampleList.size()); - for(Sample* s : sampleList) { - thePop->writeSample(*s,simfile); + simfile.write(thePop->getGenetics().getLoci()); + // save microsamples + simfile.write((size_type)microsample_list.size()); + for(Sample_base* ms : microsample_list) { + dynamic_cast(ms)->write_to_file(simfile); + } + // save samples: + simfile.write((size_type)sample_list.size()); + for(Sample_base* s : sample_list) { + dynamic_cast(s)->write_to_file(simfile); } if (thePop->isTrackingGenes()) { thePop->getGenetics().writeGeneLists(simfile); } if (checkpointInterval>0) { - simfile.write((int)checkpointList.size()); - for(Checkpoint* cp : checkpointList) { - thePop->writeCheckpoint(*cp, simfile); + simfile.write((size_type)checkpointList.size()); + for(Sample_base* cp : checkpointList) { + dynamic_cast(cp)->write_to_file(simfile); } } + simfile.close(); + std::cout << "Done.\n" << std::flush;; } // Load samples and checkpoints from results file. Don't change simulation parameters. void Simulation::load(std::string& resultsFileName) { iSimfile resFile(resultsFileName); - int fileVersion = resFile.read(); - if( fileVersion!=2 ) { - std::cout << "Wrong file version of results file " << resultsFileName << '\n'; - exit(0); - } - // read seed (not used): - resFile.read(); + // read seed (this should be kept): + seed = resFile.read(); // read (old) parameter-file copy (not used): resFile.readString(); // read loci (not used): - resFile.read(); + resFile.read(); - int sampleCount = resFile.read(); - sampleList.clear(); - sampleList.reserve(sampleCount); - for (int si=0; sireadSample(resFile)); + // read microsamples + size_type microsample_count = resFile.read(); + microsample_list.clear(); + for (size_type si=0; si(); + sample_list.clear(); + sample_list.reserve(sampleCount); + for (size_type si=0; siisTrackingGenes()) { thePop->getGenetics().readGeneLists(resFile); } if (checkpointInterval>0) { - int count = resFile.read(); + size_type count = resFile.read(); checkpointList.clear(); checkpointList.reserve(count); - for (int cpi=0; cpireadCheckpoint(resFile)); + for (size_type cpi=0; cpi sampleList; - std::vector checkpointList; + std::vector sample_list; + std::vector microsample_list; + // CheckpointList is implemented as a list of Sample_base* to be able to administrate it as such. The pointers are typecast to Population_checkpoint* (a subclass of Sample_base) when necessary. + std::vector checkpointList; void load(std::string& resultsFile); void save(); void makeFileName(int replicate); public: - Simulation(std::string& parameterFileName); + Simulation(std::string& parameterFileName, bool verbose); ~Simulation(); + seed_type get_seed() { return seed;} void initialize(int replicate); - void runFromGeneration(timeType gen); - void resume(std::string resultsFile, timeType gen, std::string newResultsFile); + void runFromGeneration(time_type gen, seed_type gen_seed, double time_so_far); + void resume(std::string resultsFile, time_type gen, std::string newResultsFile); }; + #endif /* defined(__Species__simulation__) */ diff --git a/TREES_code/space.cpp b/TREES_code/space.cpp index d70222e..95173fa 100644 --- a/TREES_code/space.cpp +++ b/TREES_code/space.cpp @@ -23,7 +23,7 @@ #include #include #include -#include "stdlib.h" +#include #include "space.h" #include "population.h" @@ -34,43 +34,124 @@ Space::Space(Population& p) : pop(p), Dims(0) {} Space::~Space() {} -DiscreteSpace::DiscreteSpace(Population& p) : +/*class Space_sample { + protected: + int pop_size; + int dims; + Space_sample(Space& S); + public: + virtual void write_to_file(oSimfile& osf)=0; +};*/ + +Space_sample::Space_sample(Space& S) { + pop_size = S.pop.size(); + dims = S.getDims(); +} + +Space_sample::~Space_sample(){ +} + +void Space_sample::write_to_file(oSimfile &osf) { + osf.write(pop_size); + osf.write(dims); +} + +Space_sample::Space_sample(iSimfile& isf) { + pop_size = isf.read(); + dims = isf.read(); +} + +Space_sample* Space_sample::create_from_file(iSimfile &isf) { + Space_types type = (Space_types)isf.read(); + switch (type) { + case Space_types::Continuous_space_type: + return new Continuous_space_sample(isf); + break; + + case Space_types::Discrete_space_type: + return new Discrete_space_sample(isf); + break; + + default: + std::cout << "Wrong space type!\n"; + exit(1); + return NULL; + break; + } +} + +Discrete_space::Discrete_space(Population& p) : Space(p) {} -bool DiscreteSpace::isDiscrete() { return true; } - -NullSpace::NullSpace(Population& p) : DiscreteSpace(p) {} -NullSpace::~NullSpace() {} -void NullSpace::initialize(int n0) {} -void NullSpace::disperse() {} -void NullSpace::prepareNewGeneration(int size) {} -void NullSpace::nextGeneration() {} -void NullSpace::addChild(int mom, int dad) {} -void NullSpace::compactData(std::vector& alive) {} -void NullSpace::addToSample(Sample& theSample) {} -void NullSpace::readToSample(Sample& theSample, iSimfile& isf, int n) {}; -int NullSpace::resumeFromCheckpoint(Checkpoint& cp, int dataIndex) {return dataIndex; } -double NullSpace::getPosition(int individual, int dim) { return 0.0; } -double NullSpace::getDist2(int ind1, int ind2) {return 0.0;} - -int NullSpace::getLinearPatch(int individual) { return 0; } -int NullSpace::popSizeInPatch(int linearPatch) { +bool Discrete_space::isDiscrete() { return true; } + +Space_sample* Discrete_space::make_sample() { + return new Discrete_space_sample(*this); +} + +/*class DiscreteSpace_sample { + protected: + std::vector linear_patches; + public: + DiscreteSpace_sample(DiscreteSpace& S); +}*/ + +Discrete_space_sample::Discrete_space_sample(Discrete_space& S) : Space_sample(S) { + linear_patches.clear(); + if (dims>0) { // S may be a NullSpace + linear_patches.reserve(pop_size); + for (int ind=0; ind((char)Discrete_space_type); + Space_sample::write_to_file(osf); // writes basic parameters + osf.writeVector(linear_patches); +} + +Discrete_space_sample::Discrete_space_sample(iSimfile& isf) : +Space_sample(isf){ + isf.readVector(linear_patches); +} + +Null_space::Null_space(Population& p) : Discrete_space(p) {} +Null_space::~Null_space() {} +void Null_space::initialize() {} +void Null_space::disperse() {} +void Null_space::prepareNewGeneration(int size) {} +void Null_space::nextGeneration() {} +void Null_space::addChild(int mom, int dad) {} +void Null_space::compactData(std::vector& alive) {} +void Null_space::resumeFromCheckpoint(Space_sample* S_s) {} +double Null_space::getPosition(int individual, int dim) { return 0.0; } +double Null_space::getDist2(int ind1, int ind2) {return 0.0;} + +int Null_space::getLinearPatch(int individual) { return 0; } +int Null_space::popSizeInPatch(int linearPatch) { if (linearPatch==0) { return pop.size(); } else { return 0; } } -int NullSpace::getPatchCount() { return 1; } +int Null_space::getPatchCount() { return 1; } +double Null_space::get_mean(int dim) { return 0.0; } +double Null_space::get_variance(int dim, double mean) { return 0.0; }; -DiscreteSpaceImp::DiscreteSpaceImp(Population& p, ParameterFile& pf) : -DiscreteSpace(p) { +Discrete_space_imp::Discrete_space_imp(Population& p, ParameterFile& pf) : +Discrete_space(p) { length = pf.getPositiveInt("size"); Dims = pf.getPositiveInt("dimensions"); std::string PdTraitname = pf.getString("p_disperse"); PDisp = pop.findTrait(PdTraitname); - if (PDisp->getDims() > 1) { + if (PDisp->get_dims() > 1) { std::cout << "Dispersal probability can not be multidimensional. Trait : " << PdTraitname << '\n'; exit(0); } @@ -105,25 +186,25 @@ DiscreteSpace(p) { exit(0); } - patches.reserve(10000*Dims); - newPatches.reserve(10000*Dims); + patches.reserve(Dims,10000); + newPatches.reserve(Dims,10000); patchPopSizes.assign(getPatchCount(), 0); generationLastCount = -1; } -DiscreteSpaceImp::~DiscreteSpaceImp() { +Discrete_space_imp::~Discrete_space_imp() { } -int DiscreteSpaceImp::getLinearPatch(int individual) { +int Discrete_space_imp::getLinearPatch(int individual) { return linearPatches[individual]; } -double DiscreteSpaceImp::getPosition(int indiviudal, int dim) { - return patches[indiviudal*Dims + dim]; +double Discrete_space_imp::getPosition(int indiviudal, int dim) { + return patches(dim,indiviudal); } -double DiscreteSpaceImp::getDist2(int ind1, int ind2) { // squared distance +double Discrete_space_imp::getDist2(int ind1, int ind2) { // squared distance double dist2=0.0; int *p1 = &getPatch(ind1, 0); int *p2 = &getPatch(ind2, 0); @@ -136,7 +217,7 @@ double DiscreteSpaceImp::getDist2(int ind1, int ind2) { // squared distance return dist2; } -int DiscreteSpaceImp::popSizeInPatch(int linearPatch) { +int Discrete_space_imp::popSizeInPatch(int linearPatch) { if (pop.getAge()!=generationLastCount) { patchPopSizes.assign(getPatchCount(), 0); for (int i=0; i1) { for (int i=0; itraitValue(i) ) { + if (rand1()traitValue(i)) { double z; int dim, dist; switch (dispType) { @@ -204,7 +285,7 @@ void DiscreteSpaceImp::disperse() { } // if (length>1) } -int DiscreteSpaceImp::treatBoundaries(int pos, int individual) { +int Discrete_space_imp::treatBoundaries(int pos, int individual) { if (pos>=length || pos<0) { switch (boundaryCondition) { case Reflective: @@ -234,17 +315,16 @@ int DiscreteSpaceImp::treatBoundaries(int pos, int individual) { return pos; } - -void DiscreteSpaceImp::prepareNewGeneration(int size) { - newPatches.clear(); +void Discrete_space_imp::prepareNewGeneration(int size) { + newPatches.assign(getDims(), 0, 0); } -void DiscreteSpaceImp::nextGeneration() { +void Discrete_space_imp::nextGeneration() { patches = newPatches; assignLinearPatches(); } -int DiscreteSpaceImp::calcLinearPatch(int individual) { +int Discrete_space_imp::calcLinearPatch(int individual) { int p = 0; for (int d=0; d Discrete_space_imp::calc_coords_from_linear(int linear) { + std::vector coords(Dims,0); + for (int d=Dims-1; d>=0; --d) { + coords.at(d) = linear % length; + linear = linear/length; + } + return coords; +} + +void Discrete_space_imp::assignLinearPatches() { if (Dims==1) { - linearPatches = patches; + linearPatches.resize(patches.get_N()); + std::memcpy(&linearPatches.at(0), patches.getX(), patches.get_N()*sizeof(int)); } else { - linearPatches.assign(patches.size(), 0); - for (int i=0; i& alive) { +void Discrete_space_imp::compactData(std::vector& alive) { + patches.compact_data(alive); int iw=0; for (int ir=0; ir(patches)); -} -void DiscreteSpaceImp::readToSample(Sample& theSample, iSimfile& isf, int n) { - theSample.addData(new XData(isf,n)); -} - -int DiscreteSpaceImp::resumeFromCheckpoint(Checkpoint &cp, int dataIndex) { - patches = (dynamic_cast&>(cp.getData(dataIndex++))).getData(); +void Discrete_space_imp::resumeFromCheckpoint(Space_sample* S_Sample) { + Discrete_space_sample* DS_sample = dynamic_cast(S_Sample); + linearPatches = DS_sample->linear_patches; + assert(linearPatches.size()==pop.size()); + patches.assign(Dims, pop.size(), 0); + for (int ind=0; ind coords = calc_coords_from_linear(linearPatches.at(ind)); + for (int d=0; dgetDims() > 1) { + if (PDisp->get_dims() > 1) { std::cout << "Dispersal probability can not be multidimensional. Trait : " << PdTraitname << '\n'; exit(0); } @@ -327,31 +418,27 @@ Space(p) { std::cout << "ContinuousSpace: Initial position outside range : " << initialPosition << '\n'; exit(0); } - pos.reserve(10000*Dims); - newPos.reserve(10000*Dims); + pos.reserve(Dims,10000); + newPos.reserve(Dims,10000); } -ContinuousSpace::~ContinuousSpace() { +Continuous_space::~Continuous_space() { } -bool ContinuousSpace::isDiscrete() { return false; } +bool Continuous_space::isDiscrete() { return false; } -ContinuousSpace::positionType& ContinuousSpace::getCoord(int indiviudal, int dim) { - return pos[indiviudal*Dims + dim]; +double Continuous_space::getPosition(int indiviudal, int dim) { + return pos(dim,indiviudal); } -double ContinuousSpace::getPosition(int indiviudal, int dim) { - return pos[indiviudal*Dims + dim]; -} - -double ContinuousSpace::getDist2(int ind1, int ind2) { // squared distance +double Continuous_space::getDist2(int ind1, int ind2) { // squared distance if (Dims==1) { - double dx = getCoord(ind1, 0)-getCoord(ind2, 0); + double dx = pos(0,ind1)-pos(0,ind2); return dx*dx; } else { double dist2 = 0.0; - positionType *p1 = &getCoord(ind1, 0); - positionType *p2 = &getCoord(ind2, 0); + positionType *p1 = &pos(0,ind1); + positionType *p2 = &pos(0,ind2); for (int d=0; dmaxPos/2) { - return maxPos-dx; - } else { - return dx; - } -}*/ - -void ContinuousSpace::initialize(int n0) { +void Continuous_space::initialize() { // Put everyone in initialPosition - pos.assign(n0*Dims, initialPosition); + pos.assign(Dims,pop.size(), initialPosition); } -ContinuousSpace::positionType ContinuousSpace::treatBoundaries(ContinuousSpace::positionType pos, int individual) { +Continuous_space::positionType Continuous_space::treatBoundaries(Continuous_space::positionType pos, int individual) { if (pos>maxPos || pos<0) { switch (boundaryCondition) { case Reflective: @@ -407,7 +484,7 @@ ContinuousSpace::positionType ContinuousSpace::treatBoundaries(ContinuousSpace:: return pos; } -void ContinuousSpace::disperse() { +void Continuous_space::disperse() { for (int i=0; itraitValue(i); @@ -419,7 +496,7 @@ void ContinuousSpace::disperse() { if (z dir(Dims,0); @@ -430,52 +507,62 @@ void ContinuousSpace::disperse() { } double C = dist/sqrt(r2); for (int d=0; d& alive) { - int iw=0; // writing position - for (int ir=0; ir& alive) { + pos.compact_data(alive); } -void ContinuousSpace::addToSample(Sample& theSample) { - // add sample spatial position: - theSample.addData(new XData(pos)); +void Continuous_space::resumeFromCheckpoint(Space_sample* S_sample) { + Continuous_space_sample* CS_sample = dynamic_cast(S_sample); + pos = CS_sample->positions; } -void ContinuousSpace::readToSample(Sample& theSample, iSimfile& isf, int n) { - theSample.addData(new XData(isf,n*Dims)); +double Continuous_space::get_mean(int dim) { + return pos.row_mean(dim); +} +double Continuous_space::get_variance(int dim, double mean) { + return pos.row_variance(dim, mean); +} + +Continuous_space_sample::Continuous_space_sample(Continuous_space& S) : Space_sample(S) { + positions = S.pos; +} + +void Continuous_space_sample::write_to_file(oSimfile& osf) { + osf.write((char)Continuous_space_type); + Space_sample::write_to_file(osf); + assert(positions.get_M()==dims); + assert(positions.get_N()==pop_size); + positions.write_to_file(osf); +} + +Continuous_space_sample::Continuous_space_sample(iSimfile& isf) : +Space_sample(isf){ + positions.read_from_file(isf); } -int ContinuousSpace::resumeFromCheckpoint(Checkpoint &cp, int dataIndex) { - pos = (dynamic_cast&>(cp.getData(dataIndex++))).getData(); - return dataIndex; +Continuous_space_sample::~Continuous_space_sample() { } +Space_sample* Continuous_space::make_sample() { + return new Continuous_space_sample(*this); +} diff --git a/TREES_code/space.h b/TREES_code/space.h index ac07302..f8c4b64 100644 --- a/TREES_code/space.h +++ b/TREES_code/space.h @@ -28,8 +28,8 @@ #include #include "parameterFile.h" -#include "sample.h" #include "trait.h" +#include "matrix.h" ///////////////////////////// // Space modules @@ -38,9 +38,11 @@ ///////////////////////////// class Population; +class Space_sample; // The space interface: class Space { + friend class Space_sample; protected: Population& pop; int Dims; @@ -48,52 +50,83 @@ class Space { Space(Population& p); virtual ~Space(); int getDims() { return Dims; } - virtual void initialize(int n0)=0; + virtual void initialize()=0; virtual void disperse()=0; virtual void prepareNewGeneration(int size)=0; virtual void nextGeneration()=0; virtual void addChild(int mom, int dad)=0; virtual void compactData(std::vector& alive)=0; - virtual void addToSample(Sample& theSample)=0; - virtual void readToSample(Sample& theSample, iSimfile& isf, int n)=0; - virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex)=0; + virtual Space_sample* make_sample()=0; + virtual void resumeFromCheckpoint(Space_sample* S_sample)=0; virtual double getPosition(int individual, int dimension)=0; virtual double getDist2(int ind1, int ind2)=0; // squared distance + virtual double get_mean(int dim)=0; + virtual double get_variance(int dim, double mean)=0; virtual bool isDiscrete()=0; }; +enum Space_types { + Discrete_space_type, + Continuous_space_type +}; + +class Space_sample { +protected: + int pop_size; + int dims; + Space_sample(Space& S); + Space_sample(iSimfile& isf); +public: + virtual ~Space_sample(); + virtual void write_to_file(oSimfile& osf); + static Space_sample* create_from_file(iSimfile& isf); +}; + // The DiscreteSpace interface -class DiscreteSpace : public Space { +class Discrete_space : public Space { + friend class Discrete_space_sample; public: - DiscreteSpace(Population& p); + Discrete_space(Population& p); //virtual ~DiscreteSpace()=0; // new functions: virtual int getLinearPatch(int individual)=0; virtual int popSizeInPatch(int linearPatch)=0; virtual int getPatchCount()=0; // From Space interface: - virtual void initialize(int n0)=0; + virtual void initialize()=0; virtual void disperse()=0; virtual void prepareNewGeneration(int size)=0; virtual void nextGeneration()=0; virtual void addChild(int mom, int dad)=0; virtual void compactData(std::vector& alive)=0; - virtual void addToSample(Sample& theSample)=0; - virtual void readToSample(Sample& theSample, iSimfile& isf, int n)=0; - virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex)=0; + virtual void resumeFromCheckpoint(Space_sample* S_sample)=0; virtual double getPosition(int individual, int dim)=0; virtual double getDist2(int ind1, int ind2)=0; // squared distance + virtual double get_mean(int dim)=0; + virtual double get_variance(int dim, double mean)=0; virtual bool isDiscrete(); + virtual Space_sample* make_sample(); }; -// A discrete space implementation, with global dispersal. -class DiscreteSpaceImp : public DiscreteSpace { +class Discrete_space_sample : public Space_sample { + friend class Discrete_space_imp; protected: - std::vector patches; - std::vector newPatches; + std::vector linear_patches; +public: + Discrete_space_sample(Discrete_space& S); + Discrete_space_sample(iSimfile& isf); + virtual ~Discrete_space_sample(); + virtual void write_to_file(oSimfile& osf); +}; + +// A discrete space implementation +class Discrete_space_imp : public Discrete_space { +protected: + Matrix patches; + Matrix newPatches; std::vector patchPopSizes; - std::vector linearPatches; - timeType generationLastCount; + std::vector linearPatches; // one per individual + time_type generationLastCount; int length; // size of space = length^Dims int initialPatch; // starting patch for all individuals (in all dimensions) Trait* PDisp; @@ -119,60 +152,62 @@ class DiscreteSpaceImp : public DiscreteSpace { } } int calcLinearPatch(int individual); + std::vector calc_coords_from_linear(int linear); void assignLinearPatches(); public: - DiscreteSpaceImp(Population& p, ParameterFile& pf); - virtual ~DiscreteSpaceImp(); - int& getPatch(int individual, int dim) {return patches[individual*Dims + dim];} + Discrete_space_imp(Population& p, ParameterFile& pf); + virtual ~Discrete_space_imp(); + int& getPatch(int individual, int dim) {return patches(dim,individual);} // DiscreteSpace interface: virtual int getLinearPatch(int individual); virtual int popSizeInPatch(int linearPatch); virtual int getPatchCount() { return std::pow(length,Dims); } // Space interface: - virtual void initialize(int n0); + virtual void initialize(); virtual void disperse(); virtual void prepareNewGeneration(int size); virtual void nextGeneration(); virtual void addChild(int mom, int dad); virtual void compactData(std::vector& alive); - virtual void addToSample(Sample& theSample); - virtual void readToSample(Sample& theSample, iSimfile& isf, int n); - virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex); + virtual void resumeFromCheckpoint(Space_sample* S_sample); virtual double getPosition(int individual, int dim); virtual double getDist2(int ind1, int ind2); // squared distance + virtual double get_mean(int dim); + virtual double get_variance(int dim, double mean); }; // The NullSpace class is a DiscreteSpace implementation // with no stored positions. // All individuals are regarded as positioned in patch 0. -class NullSpace : public DiscreteSpace { +class Null_space : public Discrete_space { public: - NullSpace(Population& p); - virtual ~NullSpace(); + Null_space(Population& p); + virtual ~Null_space(); virtual int getLinearPatch(int individual); virtual int popSizeInPatch(int linearPatch); virtual int getPatchCount(); // From Space interface: - virtual void initialize(int n0); + virtual void initialize(); virtual void disperse(); virtual void prepareNewGeneration(int size); virtual void nextGeneration(); virtual void addChild(int mom, int dad); virtual void compactData(std::vector& alive); - virtual void addToSample(Sample& theSample); - virtual void readToSample(Sample& theSample, iSimfile& isf, int n); - virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex); + virtual void resumeFromCheckpoint(Space_sample* S_s); virtual double getPosition(int individual, int dimension); virtual double getDist2(int ind1, int ind2); // squared distance + virtual double get_mean(int dim); + virtual double get_variance(int dim, double mean); }; -class ContinuousSpace : public Space { +class Continuous_space : public Space { + friend class Continuous_space_sample; protected: typedef float positionType; - std::vector pos; - std::vector newPos; + Matrix pos; + Matrix newPos; double initialPosition; // starting position for all individuals (in all dimensions) Trait* PDisp; positionType dispdist; // mean dispersal distance @@ -195,22 +230,32 @@ class ContinuousSpace : public Space { } } public: - ContinuousSpace(Population& p, ParameterFile& pf); - virtual ~ContinuousSpace(); + Continuous_space(Population& p, ParameterFile& pf); + virtual ~Continuous_space(); virtual double getPosition(int individual, int dim); - positionType& getCoord(int individual, int dim); // position with reference - virtual void initialize(int n0); + //positionType& getCoord(int individual, int dim); // position with reference + virtual void initialize(); virtual void disperse(); virtual void prepareNewGeneration(int size); virtual void nextGeneration(); virtual void addChild(int mom, int dad); virtual void compactData(std::vector& alive); - virtual void addToSample(Sample& theSample); - virtual void readToSample(Sample& theSample, iSimfile& isf, int n); - virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex); + virtual void resumeFromCheckpoint(Space_sample* S_sample); virtual double getDist2(int ind1, int ind2); // squared distance + virtual double get_mean(int dim); + virtual double get_variance(int dim, double mean); virtual bool isDiscrete(); + virtual Space_sample* make_sample(); }; - +class Continuous_space_sample : public Space_sample { + friend class Continuous_space; + protected: + Matrix positions; + public: + Continuous_space_sample(Continuous_space& S); + Continuous_space_sample(iSimfile& isf); + virtual ~Continuous_space_sample(); + virtual void write_to_file(oSimfile& osf); +}; #endif /* defined(__Species__space__) */ diff --git a/TREES_code/trait.cpp b/TREES_code/trait.cpp index 9d74d96..129c985 100644 --- a/TREES_code/trait.cpp +++ b/TREES_code/trait.cpp @@ -19,113 +19,136 @@ // Contact: jorgen.ripa@biol.lu.se // +#include +#include #include "trait.h" #include "population.h" -#include "sample.h" +#include "genotype_phenotype_map.h" #include "randgen.h" // protected constructor, used by subclasses: Trait::Trait(std::string& tname, Population& p) : -name(tname), pop(p), genes(p.getGenetics()) { - dims = 0; - lociPerDim = 0; +Matrix(0,0), name(tname), pop(p), genes(p.getGenetics()) { + GP_map = NULL; Xinit = 0; - startLocus = 0; } Trait::Trait(std::string& tname, Population& p, ParameterFile& pf) : -name(tname), pop(p), genes(p.getGenetics()) { +Matrix(0,0), name(tname), pop(p), genes(p.getGenetics()) { - // Make genomic space for this trait: - dims = pf.getPositiveInt("dimensions"); - lociPerDim = pf.getPositiveInt("loci_per_dim"); + set_dims(pf.getPositiveInt("dimensions")); + int lociPerDim = pf.getPositiveInt("loci_per_dim"); Xinit =pf.getDouble("initial_value"); - startLocus = genes.loci; - genes.loci += lociPerDim*dims; + GP_map = genes.create_GP_map(*this,lociPerDim, pf); // Reserve memory space for phenotypes (this will be expanded when necessary): - X.reserve(10000*dims); - // clear transform list: - transforms.clear(); + reserve(get_dims(),10000); } -bool Trait::isConstant() { return false; } +bool Trait::is_constant() { return false; } -void Trait::initialize(int n0) { - genes.setInitialValue(Xinit,n0,startLocus,dims,lociPerDim); +void Trait::initialize() { + assign(get_dims(), pop.size(), 0); + GP_map->initialize(Xinit); // this initializes the genes, not trait values } Trait::~Trait() { - for (int ti=0; tigenerate_phenotypes(); } -void Trait::generatePhenotypes() { - X.assign(pop.size()*dims, 0); - if (lociPerDim>0) { - for (int i=0; itransform(&X[0], X.size()); - } +traitType& Trait::traitValue(int individual, int dim) { + assert(dimadd_to_checkpoint(cp); } -void Trait::compactData(std::vector& alive) { - // compact arrays: - int iw = 0; // write index - for (int ir=0; iriw) { - for (int d=0; dread_to_checkpoint(cp,isf); } -void Trait::addToSample(Sample& s) { - s.addData(new FloatData(X)); +int Trait::resume_from_checkpoint(Checkpoint &cp, int dataIndex) { + dataIndex = GP_map->resume_from_checkpoint(cp,dataIndex); + return dataIndex; } +*/ + +////////////////////////////////////////////////////////// +// TraitConstant +/////////////////////////////////////////////////// TraitConstant::TraitConstant(std::string& name, Population& p, ParameterFile& pf): Trait(name,p) { // A constant trait is implemented as a trait with zero loci: - dims = pf.getPositiveInt("dimensions"); - lociPerDim = 0; + set_dims(pf.getPositiveInt("dimensions")); Xinit =pf.getDouble("initial_value"); - startLocus = genes.loci; // Reserve memory space for phenotypes (this will be expanded when necessary): - X.reserve(10000*dims); - // clear transform list: - transforms.clear(); + reserve(get_dims(),10000); } -bool TraitConstant::isConstant() { return true; } +bool TraitConstant::is_constant() { return true; } -void TraitConstant::initialize(int n0) { +void TraitConstant::initialize() { // do nothing } void TraitConstant::generatePhenotypes() { - X.assign(pop.size()*dims, Xinit); + // We don't need a GP map for this + assign(get_dims(), pop.size(), Xinit); +} + +/* +void TraitConstant::add_to_checkpoint(Checkpoint &cp) { + // do nothing +} + +void TraitConstant::read_to_checkpoint(Checkpoint &cp, iSimfile &isf) { + // niente +} + +int TraitConstant::resume_from_checkpoint(Checkpoint &cp, int dataIndex){ + return dataIndex; +} +*/ + +/* class Trait_sample { + protected: + int dims; + int size; + std::vector trait_values; + public: + Trait_sample(Trait& T); + int get_dims() {return dims;} + int get_size() {return size;} + void write_to_file(oSimfile &osf); +};*/ + +Trait_sample::Trait_sample(Trait & T) : Matrix(T){ + //assert(dims*size == trait_values.size()); } +//void Trait_sample::write_to_file(oSimfile &osf) { +// osf.write(dims); +// osf.write(size); +// osf.writeArray(&trait_values.at(0), dims*size); +//} + +Trait_sample::Trait_sample(iSimfile& isf) : Matrix(isf){ +// dims = isf.read(); +// size = isf.read(); +// isf.readArray(trait_values, dims*size); +} diff --git a/TREES_code/trait.h b/TREES_code/trait.h index e523f77..b2a8648 100644 --- a/TREES_code/trait.h +++ b/TREES_code/trait.h @@ -29,72 +29,70 @@ #include "transform.h" #include "parameterFile.h" +#include "types.h" // defines traitType +#include "matrix.h" class Population; class Genetics; -class Sample; +class Genotype_Phenotype_map; +class Checkpoint; +class iSimfile; ///////////////////////////// // Trait class -// A trait is really a genotype-phenotype map -// The Trait class is friends with the Genetics class +// Contains all trait values +// Has a genotype-phenotype map ///////////////////////////// -class Trait { +class Trait : public Matrix { + friend class Population; + friend class Genotype_Phenotype_map; + friend class Diallelic_GP_map; + friend class ContinuousAlleles_GP_map; + friend class Omnigenic_GP_map; + friend class Trait_sample; protected: Trait(std::string& name, Population& p); std::string name; Population& pop; Genetics& genes; - int dims; // number of dimensions - std::vector X; // phenotypic data of entire population - double Xinit; // initial value - typedef std::vector::iterator Xiter; - int startLocus; - int lociPerDim; // loci per dimension - std::vector transforms; + Genotype_Phenotype_map* GP_map; + void set_dims(int d) { M = d;} + traitType Xinit; // initial value public: Trait(std::string& name, Population& p, ParameterFile& pf); ~Trait(); - void addTransform( Transform* t); - virtual void initialize(int n0); + virtual bool is_constant(); + virtual void initialize(); virtual void generatePhenotypes(); - void compactData(std::vector& alive); - double& traitValue(int individual, int dim=0); - // {return X[individual];} // not worth it! - //std::vector& getX() { return X;} + void compactData(std::vector& alive) { compact_data(alive);} + traitType& traitValue(int individual, int dim=0); std::string& getName() { return name; } - int getDims() { return dims; } - void addToSample(Sample& s); - virtual bool isConstant(); + int get_dims() { return M; } + Genotype_Phenotype_map* get_GP_map() { return GP_map;} +// virtual void add_to_checkpoint(Checkpoint& cp); +// virtual void read_to_checkpoint(Checkpoint& cp, iSimfile& isf); +// virtual int resume_from_checkpoint(Checkpoint& cp, int dataIndex); }; class TraitConstant : public Trait { public: TraitConstant(std::string& name, Population& p, ParameterFile& pf); - virtual void initialize(int n0); + virtual bool is_constant(); + virtual void initialize(); virtual void generatePhenotypes(); - virtual bool isConstant(); +// virtual void add_to_checkpoint(Checkpoint& cp); +// virtual void read_to_checkpoint(Checkpoint& cp, iSimfile& isf); +// virtual int resume_from_checkpoint(Checkpoint& cp, int dataIndex); }; -/* -class MultiTraitIterator { +class Trait_sample : public Matrix{ protected: - int traitCount; // height of virtual array - double** traitp; // list of pointers to trait values public: - MultiTraitIterator(std::vector& traitList, int startpos=0); - MultiTraitIterator(std::vector& traitList, int startpos=0); - ~MultiTraitIterator(); - inline double& operator [](int t) {return *traitp[t];} - inline MultiTraitIterator& operator ++() { - for (int t=0; t +#include + #include "transform.h" +#include "genotype_phenotype_map.h" #include "fastmath.h" #include "parameterFile.h" #include "randgen.h" +#include "trait.h" +#include "genetics.h" Transform::~Transform() {} @@ -32,40 +37,69 @@ LinearTransform::LinearTransform(ParameterFile& pf) { offset = pf.getDouble("offset"); scale = pf.getDouble("scale"); } + +LinearTransform::LinearTransform(traitType offset, traitType scale) { + this->offset = offset; + this->scale = scale; +} + LinearTransform::~LinearTransform() {} -void LinearTransform::transform(double* x, size_t count) { +void LinearTransform::transform(traitType* x, size_t count) { for (size_t i=0; i(&GP_map); + if (the_map) { + traitType min_value = pf.getDouble("min"); + traitType max_value = pf.getDouble("max"); + int L = the_map->get_loci_per_dim(); // untransformed values [-L,L] + offset = L + min_value; + scale = (max_value-min_value)/2.0/L; + } else { + std::cout << "Error! \n" << "Trait " << GP_map.get_trait().getName() << ": A Range transform requires Diallelic Genetics.\n"; + exit(1); + } +} + +Range::~Range() {}; + AbsTransform::AbsTransform(ParameterFile& pf) {} AbsTransform::~AbsTransform() {} -void AbsTransform::transform(double* x, size_t count) { +void AbsTransform::transform(traitType* x, size_t count) { for (size_t i=0; i #include "parameterFile.h" +#include "types.h" // defines traitType +#include "fastmath.h" + +class Genotype_Phenotype_map; // neede here to prevent circular header files. class Transform { public: - virtual void transform(double* x, size_t count)=0; + virtual void transform(traitType* x, size_t count)=0; + virtual double transform_single(double x)=0; virtual ~Transform(); }; class LinearTransform : public Transform { - double offset; - double scale; +protected: + traitType offset; + traitType scale; public: LinearTransform(ParameterFile& pf); + LinearTransform(traitType offset, traitType scale); virtual ~LinearTransform(); - virtual void transform(double* x, size_t count); + virtual void transform(traitType* x, size_t count); + virtual double transform_single(double x); +}; + + +class Range : public LinearTransform { +public: + Range(ParameterFile& pf, Genotype_Phenotype_map& GP_map); + virtual ~Range(); }; + class AbsTransform : public Transform { public: AbsTransform(ParameterFile& pf); virtual ~AbsTransform(); - virtual void transform(double* x, size_t count); + virtual void transform(traitType* x, size_t count); + virtual double transform_single(double x); }; class LogisticTransform : public Transform { - double xmin; - double xmax; + traitType xmin; + traitType scale; + Fastexp fexp; public: LogisticTransform(ParameterFile& pf);; virtual ~LogisticTransform(); - virtual void transform(double* x, size_t count); + virtual void transform(traitType* x, size_t count); + virtual double transform_single(double x); }; class NormalDeviate : public Transform { - double SD; + traitType SD; public: NormalDeviate(ParameterFile& pf); virtual ~NormalDeviate(); - virtual void transform(double* x, size_t count); + virtual void transform(traitType* x, size_t count); + virtual double transform_single(double x); }; - - #endif /* defined(__Species__transform__) */ diff --git a/TREES_code/types.h b/TREES_code/types.h index 5461b24..489afed 100644 --- a/TREES_code/types.h +++ b/TREES_code/types.h @@ -25,7 +25,10 @@ #include -//typedef uint64_t idType; -typedef int64_t timeType; - +typedef uint64_t id_type; +typedef int64_t time_type; +typedef uint64_t seed_type; +typedef uint32_t size_type; +typedef float traitType; +typedef float geneStorage; // used for file storage of gene effects #endif