Commit 44f326c3 authored by jlopez's avatar jlopez
Browse files

Change project to R package

parents e282efd2 2e6eb9c7
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 src/sources/InputParser.cpp src/headers/InputParser.h)
add_executable(Population ${SOURCE_FILES})
find_package(GSL REQUIRED)
target_link_libraries(Population libgsl.a libgslcblas.a)
\ No newline at end of file
File added
#!/bin/bash
#$ -cwd
#$ -j y
#$ -S /bin/bash
#$ -N popu
#$ -V
module load gcc5.3
/home/jlopez/work_pop/Population -i /home/jlopez/work_pop/paramC.dat -o $TMPDIR/ -g 500 -m /home/jlopez/work_pop/
cp $TMPDIR/*_rescuers ./
This diff is collapsed.
0.0595168352221398 -0.00200580999852227 -0.00161243484323404 -0.00138749847828824 0.00253150110707747
-0.00200580999852227 0.0589221483895156 0.00322204140865314 -0.00109475248990945 -0.00203834308770295
-0.00161243484323404 0.00322204140865314 0.0644277222040999 0.00179266470707301 -0.000226647497984715
-0.00138749847828824 -0.00109475248990945 0.00179266470707301 0.0625259614913305 -0.000121287903527481
0.00253150110707747 -0.00203834308770295 -0.000226647497984715 -0.000121287903527481 0.0642562268177967
0.0638600572001953 -0.000893600318243619 0.00352203303962418 -0.00165509383509936 -0.00115572623507183
-0.000893600318243619 0.0643683920748695 -0.00090814298906903 0.00436001515282374 0.00194350598048018
0.00352203303962418 -0.00090814298906903 0.0640014380005563 -0.00286556487521142 0.00133350102848826
-0.00165509383509936 0.00436001515282374 -0.00286556487521142 0.0690277135634814 -0.00098449127692254
-0.00115572623507183 0.00194350598048018 0.00133350102848826 -0.00098449127692254 0.0624068628542177
#Comment_with_no_space
#n
5
#ne
10
#na
2000
#N0
100
#Bmax
2
#d
1
#n_rep
50000
#mu
0.01
#t_max
307.674701122751
#log10Pextlim
-10
#nbt
100
#so
1.15
#nb_time_steps
1848
#Es
0.01
#dtprint
0.166666666666667
#natba
10
-0.225374562894369
-2.0520651180874
5.40592540892945
1.79027686017605
0.363782519858795
//
// 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);
double trace() const;
void transpose() const;
GMatrix* mult(double val);
GMatrix* div(double val);
GMatrix* sqrt();
unsigned long getRowSize() const;
unsigned long getColSize() const;
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 05/10/17.
//
#ifndef POPULATION_INPUTPARSER_H
#define POPULATION_INPUTPARSER_H
#include <string>
#include <vector>
class InputParser {
public:
InputParser (int &argc, char **argv);
virtual ~InputParser() = default;
const std::string& getCmdOption(const std::string &option) const;
bool cmdOptionExists(const std::string &option) const;
private:
std::vector <std::string> m_tokens;
};
#endif //POPULATION_INPUTPARSER_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> m_P;
std::vector<double> m_R;
std::vector<double> m_T;
std::vector<unsigned long> m_NM;
std::vector<double> m_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>
#include "GMatrix.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
}
static void showUsage() {
std::cerr << "Usage: \n"
<< "\t-h\t\tShow this help message\n"
<< "\t-i INPUT\tSpecify the input path\n"
<< "\t-o OUTPUT\tSpecify the output path\n"
<< "\t-g MAX_GENO\tSpecify the max genotype\toptional & default 500\n"
<< "\t-m MATRIX\tSpecify the input matrix path\tpath to matA.dat, matSb.dat & vectxo.dat\n"
<< "\t-r RATE\tIf you want the vector of positive rate\toptional\n"
<< std::endl;
exit(-1);
}
static GMatrix* multiplyMatrix(GMatrix* ma, GMatrix* mb) {
GMatrix* mc = new GMatrix(ma->getRowSize(), mb->getColSize());
if(ma->getColSize() != mb->getRowSize()) {
return mc;
}
for (unsigned long i = 0; i < ma->getRowSize(); i++) {
for (unsigned long j = 0; j < mb->getColSize(); j++) {
double sum = 0;
for (unsigned long k = 0; k < mb->getRowSize(); k++) {
sum += ma->getValue(i, k) * mb->getValue(k, j);
}
mc->setValue(i, j, sum);
}
}
}
};
#endif //POPULATION_UTILS_H
#include <iostream>
#include "Population.h"
#include "Utils.h"
#include "Rescuers.h"
......@@ -471,4 +470,5 @@ int main(int argc, char *argv[]) {
std::cout << "Finish" << std::endl;
return 0;
}
//
// 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>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <cmath>
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());
}
double GMatrix::trace() const {
gsl_vector_view v = gsl_matrix_diagonal(m_matrix);
unsigned long size = std::min(m_matrix->size1, m_matrix->size2);
double sum = 0.0;
for (unsigned long i = 0; i < size; i++) {
sum += gsl_vector_get(&v.vector, i);
}
return sum;
}
void GMatrix::transpose() const {
gsl_matrix_transpose(m_matrix);
}
GMatrix* GMatrix::mult(double val) {
for (unsigned long i = 0; i < m_matrix->size1; ++i) {
for (unsigned long j = 0; j <m_matrix->size2 ; ++j) {
double v = getValue(i, j) * val;
setValue(i, j, v);
}
}
return this;
}
GMatrix* GMatrix::div(double val) {
for (unsigned long i = 0; i < m_matrix->size1; ++i) {
for (unsigned long j = 0; j <m_matrix->size2 ; ++j) {
double v = val / getValue(i, j);
setValue(i, j, v);
}
}
return this;
}
GMatrix* GMatrix::sqrt() {
for (unsigned long i = 0; i < m_matrix->size1; ++i) {
for (unsigned long j = 0; j <m_matrix->size2 ; ++j) {
double v = std::sqrt(getValue(i, j));
setValue(i, j, v);
}
}
return this;