Commit 8ff32bd4 authored by jlopez's avatar jlopez
Browse files

Refactor C to C++

parents
cmake_minimum_required(VERSION 3.8)
project(Population)
set(CMAKE_CXX_STANDARD 11)
set(SOURCE_FILES src/main.cpp src/sources/Population.cpp src/headers/Population.h src/sources/Parameters.cpp src/headers/Parameters.h src/sources/GMatrix.cpp src/headers/GMatrix.h src/sources/GRandom.cpp src/headers/GRandom.h src/headers/Utils.h src/sources/DynamicSave.cpp src/headers/DynamicSave.h src/sources/Rescuers.cpp src/headers/Rescuers.h)
add_executable(Population ${SOURCE_FILES})
find_package(GSL REQUIRED)
target_link_libraries(Population GSL::gsl GSL::gslcblas)
\ No newline at end of file
//
// Created by jimmy on 03/10/17.
//
#ifndef POPULATION_DYNAMICSAVE_H
#define POPULATION_DYNAMICSAVE_H
#include <vector>
class DynamicSave {
public:
DynamicSave(unsigned long size);
virtual ~DynamicSave();
void print_output();
private:
std::vector<double> outtime;
std::vector<unsigned long> outNp;
std::vector<double> outEr;
std::vector<double> outVr;
std::vector<unsigned long> outne;
std::vector<unsigned long> outng;
std::vector<double> outmaxr;
std::vector<double> outmeanr;
std::vector<unsigned long> ntsr;
};
#endif //POPULATION_DYNAMICSAVE_H
//
// Created by jimmy on 02/10/17.
//
#ifndef POPULATION_INPUTMATRIX_H
#define POPULATION_INPUTMATRIX_H
#include <string>
#include<gsl/gsl_matrix.h>
#include "Parameters.h"
class GMatrix {
public:
GMatrix(unsigned long rows, unsigned long cols);
GMatrix(std::string file, unsigned long rows, unsigned long cols, bool reverse);
virtual ~GMatrix();
gsl_matrix *getMatrix() const;
double getValue(unsigned long i, unsigned long j) const;
void setValue(unsigned long i, unsigned long j, double value);
virtual void sub(const GMatrix* gMatrix);
private:
gsl_matrix* m_matrix;
};
#endif //POPULATION_INPUTMATRIX_H
//
// Created by jimmy on 03/10/17.
//
#ifndef POPULATION_GRANDOM_H
#define POPULATION_GRANDOM_H
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
class GRandom {
public:
GRandom();
virtual ~GRandom();
virtual double exponential(double mu);
virtual double flat(double a, double b);
virtual unsigned long poisson(double mu);
private:
gsl_rng* m_rgsl;
};
#endif //POPULATION_GRANDOM_H
//
// Created by jimmy on 02/10/17.
//
#ifndef POPULATION_PARAMETERS_H
#define POPULATION_PARAMETERS_H
#include <string>
#include <vector>
#include <ostream>
class Parameters {
public:
Parameters(bool dynamic, std::string file);
virtual ~Parameters() = default;
unsigned long getN() const;
unsigned long getNe() const;
unsigned long getNa() const;
unsigned long getN0() const;
double getBmax() const;
double getD() const;
unsigned long getN_rep() const;
double getMu() const;
double getT_max() const;
int getLog10Pextlim() const;
unsigned long getNbt() const;
double getSo() const;
unsigned long getNb_time_steps() const;
double getEs() const;
double getDtprint() const;
unsigned long getNatba() const;
private:
void init(bool dynamic, std::vector<std::string> params);
unsigned long m_n;
unsigned long m_ne;
unsigned long m_na;
unsigned long m_N0;
double m_Bmax;
double m_d;
unsigned long m_n_rep;
double m_mu;
double m_t_max;
int m_log10Pextlim;
unsigned long m_nbt;
double m_so;
unsigned long m_nb_time_steps;
double m_Es;
double m_dtprint;
unsigned long m_natba;
public:
friend std::ostream &operator<<(std::ostream &os, const Parameters &parameters);
};
#endif //POPULATION_PARAMETERS_H
//
// Created by jimmy on 02/10/17.
//
#ifndef POPULATION_POPULATION_H
#define POPULATION_POPULATION_H
#include "GMatrix.h"
#include <vector>
#include "GRandom.h"
class Population {
public:
Population(unsigned long rows, unsigned long cols);
virtual ~Population();
double getValue(unsigned long i, unsigned long j);
void setValue(unsigned long i, unsigned long j, double value);
double getBg(unsigned long index);
void setBg(unsigned long index, double value);
double getrg(unsigned long index);
void setrg(unsigned long index, double value);
double getMaxrg();
unsigned long getNg(unsigned long index);
void setNg(unsigned long index, unsigned long value);
double getTg(unsigned long index);
void setTg(unsigned long index, double value);
double getnm(unsigned long index);
void setnm(unsigned long index, unsigned long value);
double getBp() const;
void setBp(double Bp);
unsigned long getNp() const;
void setNp(unsigned long Np);
double getNpc() const;
void setNpc(double Npc);
unsigned long getng() const;
void setng(unsigned long ng);
virtual double compute_Er();
virtual void evalbirthrate(GMatrix* Sb, double BMax, const GMatrix* x0, unsigned long n, unsigned long nl);
virtual void birth(GRandom* rnd, double t, double u, GMatrix *A, double mu, unsigned long na, GMatrix *Sb, double Bmax, GMatrix *x0, unsigned long n, double d, double dt);
virtual void death(GRandom* rnd, double u, double d, unsigned long n, double dt);
virtual unsigned long compute_imaxr(unsigned long max);
private:
GMatrix* m_matrix;
std::vector<double> m_Bg; // Bg[i] = birth rate genotype i
std::vector<double> m_rg; // rg[i] = growth rate genotype i
std::vector<unsigned long> m_Ng; // Ng[i] = number of indiv with genotype i
std::vector<double> m_Tg; //Tg[i] = time at which genotype i appears
std::vector<unsigned long> m_nm; // nombre de mutations cumulées
double m_Bp; // population birth rate
unsigned long m_Np; // population size
double m_Npc; // cumulative population size
unsigned long m_ng; // nombre de génotypes différents
};
#endif //POPULATION_POPULATION_H
//
// Created by jimmy on 03/10/17.
//
#ifndef POPULATION_RESCUERS_H
#define POPULATION_RESCUERS_H
#include <vector>
class Rescuers {
public:
Rescuers(unsigned long size);
virtual ~Rescuers() = default;
unsigned long getP(unsigned long index) const;
double getR(unsigned long index) const;
double getT(unsigned long index) const;
unsigned long getNM(unsigned long index) const;
double getNPC(unsigned long index) const;
void setP(unsigned long index, unsigned long value);
void setR(unsigned long index, double value);
void setT(unsigned long index, double value);
void setNM(unsigned long index, unsigned long value);
void setNPC(unsigned long index, double value);
private:
std::vector<unsigned long> P;
std::vector<double> R;
std::vector<double> T;
std::vector<unsigned long> NM;
std::vector<double> NPC;
};
#endif //POPULATION_RESCUERS_H
//
// Created by jimmy on 03/10/17.
//
#ifndef POPULATION_UTILS_H
#define POPULATION_UTILS_H
#include <iostream>
#include <string>
#include <gsl/gsl_blas.h>
class Utils {
public:
static void debug(const std::string& message) {
std::cout << message << std::endl;
}
static void dgemmTN(GMatrix* A, GMatrix* B, GMatrix* C) {
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, A->getMatrix(), B->getMatrix(), 0.0, C->getMatrix());
}
static void dgemmNN(GMatrix* A, GMatrix* B, GMatrix* C) {
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A->getMatrix(), B->getMatrix(), 0.0, C->getMatrix());
}
static std::string separator() {
#ifdef _WIN32
return "\\";
#else
return "/";
#endif
}
};
#endif //POPULATION_UTILS_H
#include <iostream>
#include "headers/Population.h"
#include "headers/Parameters.h"
#include "headers/GMatrix.h"
#include "headers/GRandom.h"
#include "headers/Utils.h"
#include "headers/DynamicSave.h"
#include "headers/Rescuers.h"
#include <cmath>
#include <fstream>
int main() {
const bool SAVE_DYN = false;
const bool save_output_simulation = true;
const bool debug = false;
const unsigned long MAX_GENO = 500;
const int PARAM_INPUT = 16;
const std::string file = "/home/jimmy/jimmy/projets/guillaume/old_src/paramC.dat";
const std::string output = "/home/jimmy/Bureau/";
Parameters param(SAVE_DYN, file);
GMatrix A("/home/jimmy/jimmy/projets/guillaume/old_src/matA.dat", param.getNa(), param.getN(), true);
GMatrix Sb("/home/jimmy/jimmy/projets/guillaume/old_src/matSb.dat", param.getN(), param.getN(), false);
GMatrix x0("/home/jimmy/jimmy/projets/guillaume/old_src/vectxo.dat", param.getN(), 1, false);
GRandom grd;
const unsigned long nrnt = param.getNb_time_steps() * param.getN_rep();
Rescuers rescuers(param.getN_rep());
for(unsigned long i_rep = 0; i_rep < param.getN_rep(); i_rep++) {
std::cout << "i_rep: " << i_rep << std::endl;
Population pop(MAX_GENO, param.getN());
pop.setng(1);
pop.setNg(0, param.getN0());
pop.setTg(0, 0.0);
pop.evalbirthrate(&Sb, param.getBmax(), &x0, param.getN(), 0);
pop.setrg(0, pop.getBg(0) - param.getD());
pop.setNp(param.getN0());
pop.setNpc(0.0);
pop.setBp(pop.getBg(0)*param.getN0());
pop.setnm(0, 0);
double Er = pop.compute_Er();
// initialize times
double t = 0.0; // current time
unsigned int nb = 0; // number of birth since the last rescue test
rescuers.setP(i_rep, 1);
unsigned long imaxr = 0;
// simulation loop on time
while( t < param.getT_max()) {
// total rate of events and time at which next event occurs
double total_rate = pop.getBp() + param.getD() * pop.getNp();
double dt = grd.exponential(total_rate);
t += dt;
double u = grd.flat(0.0, total_rate);
if(u < pop.getBp()) {
if (pop.getng() >= MAX_GENO) {
std::cout << "ERROR: max number of genotypes reached\nIncrease MAX_GENO constant\n" << std::endl;
//TODO error
} else {
pop.birth(&grd, t, u, &A, param.getMu(), param.getNa(), &Sb, param.getBmax(), &x0, param.getN(), param.getD(), dt);
nb++;
if(nb > param.getNbt()) { // if nbt birth, test is there is a rescue
nb = 0;
Er = pop.compute_Er();
if(Er > 0) {
if(pop.getNp() * log(param.getD() / (Er+param.getD())) < param.getLog10Pextlim()) {
break;
}
}
}
}
} else {
pop.death(&grd, u, param.getD(), param.getN(), dt);
if(pop.getNp() == 0) {
rescuers.setP(i_rep, 0);
break;
}
}
} // end of simulation loop on time
// test if simul stops because t_max reached (in this case, failed simul)
if (t >= param.getT_max()) {
std::cout << "ERROR: t_max reached\nIncrease factsecuext parameter\n" << std::endl;
//TODO error
exit(1);
}
imaxr = pop.compute_imaxr(imaxr);
if(pop.getNp() != 0) {
rescuers.setT(i_rep, pop.getTg(imaxr));
rescuers.setR(i_rep, pop.getrg(imaxr));
rescuers.setNM(i_rep, static_cast<unsigned long>(pop.getnm(imaxr)));
} else {
rescuers.setT(i_rep, param.getT_max());
rescuers.setR(i_rep, 0.0);
rescuers.setNM(i_rep, 0);
}
rescuers.setNPC(i_rep, pop.getNpc());
}
if(save_output_simulation) {
std::ofstream outfile;
const std::string pathFile = output + Utils::separator()
+ "No" + std::to_string(param.getN0())
+ "_n" + std::to_string(param.getN())
+ "_ne" + std::to_string(param.getNe())
+ "_U" + std::to_string(param.getMu())
+ "_Es" + std::to_string(param.getEs())
+ "_so" + std::to_string(param.getSo())
+ "_rescuers";
outfile.open(pathFile);
outfile.clear();
outfile << "list(c(" << rescuers.getP(0);
for (unsigned long i = 1; i < param.getN_rep(); i++) {
outfile << ", " << rescuers.getP(i);
}
outfile << "), \nc(" << rescuers.getR(0);
for (unsigned long i = 1; i < param.getN_rep(); i++) {
outfile << ", " << rescuers.getR(i);
}
outfile << "), \nc(" << rescuers.getT(0);
for (unsigned long i = 1; i < param.getN_rep(); i++) {
outfile << ", " << rescuers.getT(i);
}
outfile << "), \nc(" << rescuers.getNM(0);
for (unsigned long i = 1; i < param.getN_rep(); i++) {
outfile << ", " << rescuers.getNM(i);
}
outfile << "), \nc(" << rescuers.getNPC(0);
for (unsigned long i = 1; i < param.getN_rep(); i++) {
outfile << ", " << rescuers.getNPC(i);
}
outfile << "))";
outfile.close();
}
return 0;
}
\ No newline at end of file
//
// Created by jimmy on 03/10/17.
//
#include "../headers/DynamicSave.h"
DynamicSave::DynamicSave(unsigned long size) {
outtime.resize(size);
outNp.resize(size);
outEr.resize(size);
outVr.resize(size);
outne.resize(size);
outng.resize(size);
outmaxr.resize(size);
outmeanr.resize(size);
ntsr.resize(size);
}
DynamicSave::~DynamicSave() {
}
void DynamicSave::print_output() {
}
//
// Created by jimmy on 02/10/17.
//
#include "../headers/GMatrix.h"
#include <string>
#include <fstream>
#include <sstream>
#include <iostream>
GMatrix::GMatrix(unsigned long rows, unsigned long cols) {
m_matrix = gsl_matrix_alloc(rows, cols);
for (unsigned long i = 0; i < rows; i++) {
for (unsigned long j = 0; j < cols; j++) {
setValue(i, j, 0);
}
}
}
GMatrix::GMatrix(std::string file, unsigned long rows, unsigned long cols, bool reverse) {
m_matrix = gsl_matrix_alloc(rows, cols);
if(reverse) {
std::ifstream infile(file);
std::string line;
unsigned long i = 0;
while (getline(infile,line)) // Read a line
{
std::istringstream iss(line);
double value;
unsigned long j = 0;
while (iss >> value) { // Read columns
gsl_matrix_set(m_matrix, j, i, value);
j++;
}
i++;
}
infile.close();
} else {
std::ifstream infile(file);
std::string line;
unsigned long i = 0;
while (getline(infile,line)) // Read a line
{
std::istringstream iss(line);
double value;
unsigned long j = 0;
while (iss >> value) { // Read columns
gsl_matrix_set(m_matrix, i, j, value);
j++;
}
i++;
}
infile.close();
}
}
GMatrix::~GMatrix() {
gsl_matrix_free(m_matrix);
}
double GMatrix::getValue(unsigned long i, unsigned long j) const {
return gsl_matrix_get(m_matrix, i, j);
}
void GMatrix::setValue(unsigned long i, unsigned long j, double value) {
gsl_matrix_set(m_matrix, i, j, value);
}
gsl_matrix *GMatrix::getMatrix() const {
return m_matrix;
}
void GMatrix::sub(const GMatrix *gMatrix) {
gsl_matrix_sub(m_matrix, gMatrix->getMatrix());
}
//
// Created by jimmy on 03/10/17.
//
#include "../headers/GRandom.h"
#include <time.h>
GRandom::GRandom() {
const gsl_rng_type *Tgsl;
gsl_rng_env_setup();
Tgsl = gsl_rng_default;
m_rgsl = gsl_rng_alloc(Tgsl);
gsl_rng_set (m_rgsl, (unsigned long) time (nullptr));
}
GRandom::~GRandom() {
delete m_rgsl;