Commit e282efd2 authored by jlopez's avatar jlopez
Browse files

Add project

parents
^.*\.Rproj$
^\.Rproj\.user$
# History files
.Rhistory
.Rapp.history
# Session Data files
.RData
# Example code in package build process
*-Ex.R
# Output files from R CMD build
/*.tar.gz
# Output files from R CMD check
/*.Rcheck/
# RStudio files
.Rproj.user/
# produced vignettes
vignettes/*.html
vignettes/*.pdf
# OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
.httr-oauth
# knitr and R markdown default cache directories
/*_cache/
/cache/
# Temporary files created by R markdown
*.utf8.md
*.knit.md
.Rproj.user
Package: PopulationR
Type: Package
Title: What the Package Does (Title Case)
Version: 0.1.0
Author: Who wrote it
Maintainer: The package maintainer <yourself@somewhere.net>
Description: More about what it does (maybe more than one line)
Use four spaces when indenting paragraphs within the Description.
License: What license is it under?
Encoding: UTF-8
LazyData: true
Imports: Rcpp (>= 0.12.13), MASS (>= 7.3-47), DistributionUtils (>= 0.5-1)
LinkingTo: Rcpp
exportPattern("^[^\\R]")
importFrom(Rcpp, evalCpp)
useDynLib(PopulationR, .registration = TRUE)
Version: 1.0
RestoreWorkspace: Default
SaveWorkspace: Default
AlwaysSaveHistory: Default
EnableCodeIndexing: Yes
UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8
RnwWeave: Sweave
LaTeX: pdfLaTeX
AutoAppendNewline: Yes
StripTrailingWhitespace: Yes
BuildType: Package
PackageUseDevtools: Yes
PackageInstallArgs: --no-multiarch --with-keep.source
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
Rrun <- function(max_geno, n, ne, na, N0, Bmax, d, n_rep, mu, t_max, log10Pextlim, nbt, so, nb_time_steps, Es, dtprint, natba, RSb, RA, Rxo, have_rate, have_rescuers, output_path, verbose) {
.Call(`_PopulationR_Rrun`, max_geno, n, ne, na, N0, Bmax, d, n_rep, mu, t_max, log10Pextlim, nbt, so, nb_time_steps, Es, dtprint, natba, RSb, RA, Rxo, have_rate, have_rescuers, output_path, verbose)
}
#' Run the simulation
#' @param max_geno integer : maximum genotype
#' @param range array : range
#' @param No integer : initial pop size
#' @param log10Pextlim integer (negative) : Log10 of the probability of extinction under which we consider that pop is rescued
#' @param rmax double : growth rate in exponential phase at the optimum (ro ~ Bmax - d)
#' @param factsecuext double :
#' @param nbt double : number of birth before testing for a rescue
#' @param ne integer :
#' @param n integer :
#' @param Es double :
#' @param Uo double :
#' @param na integer :
#' @param natba integer : Min nb of genotype copies to consider this genotype as advantageous
#' @param nbreps integer :
#' @param verbose logical : if you want a verbose program
#' @return array :
simulation <- function(max_geno = 500, range = c(-0.05, -0.1, -0.15),
No = 100, log10Pextlim = -10, rmax = 1, factsecuext = 10, nbt = 100,
ne = 10, n = 5, Es = 0.01, Uo = 0.01, na = 2000, natba = 10, nbreps = 50000,
have_rate = FALSE, have_rescuers = FALSE, output_path = "", verbose = FALSE)
{
library(MASS);
library(DistributionUtils);
print("Simulation Start")
rescuersL <- list()
d <- rmax # DOUBLE.
Bmax <- rmax + d # DOUBLE. bmax in order to get ro=Bmax - d
dt <- 50 / (No * (Bmax + d)) # DOUBLE.
rovals <- range * rmax
sovals <- rmax - rovals
nsovals <- length(sovals)
simname <- paste("No", No, "_n", n, "_ne", ne, "_U", Uo, "_Es", Es, "_so", sep="")
#Rprint(verbose, simname)
Ur <- numeric()
q <- numeric()
# loop on the values of so
for(iso in 1:nsovals) {
# set input
so <- sovals[iso]
tmax <- factsecuext * log(No+1) / abs(Bmax - so - d);
nbtimesteps = floor(tmax / dt) + 2; # nb time step printed (+ 2 because we include time 0 and t_max)
param <- c(n, ne, na, No, Bmax, d, nbreps, Uo, tmax, log10Pextlim, nbt, so, nbtimesteps, Es, dt, natba)
# build the landscape
Ini <- Rbuiltlandscape(n, ne, Es, so, na)
# Ur <- c(Ur, Uo / 2 * (1 - exp(- (B(Bmax, Ini$A, Ini$xo) - d) / d)))
L <- Ini$Sb %*% Ini$M
q <- c(q, (Ini$xo %*% L %*% L %*% Ini$xo) / (Ini$xo %*% L %*% Ini$xo) * Rtrace(L) / Rtrace(L %*% L))
# birth and death process replicated nbreps times
print(paste("running so=", so, "... Start", sep=""))
# ============================
# Start C++
# ============================
rescuersL[[iso]] <- Rrun(max_geno, n, ne, na, No, Bmax, d, nbreps, Uo, tmax, log10Pextlim, nbt, so, nbtimesteps, Es, dt, natba, Ini$Sb, Ini$A, Ini$xo, have_rate, have_rescuers, output_path, verbose)
# ============================
# End C++
# ============================
print(paste("running so=", so, "... End", sep=""))
}
pdf(paste(output_path,"distT1stepNo",No,"Uo",Uo,"n",n,"ne",ne,".pdf",sep=""), height=5, width=15)
layout(matrix(1:3, 1, 3))
for(iso in 1:nsovals) {
Prescuers <- rescuersL[[iso]][["P"]]
Trescuers <- rescuersL[[iso]][["T"]]
NMrescuers <- rescuersL[[iso]][["NM"]]
Trescuershist <- numeric()
for (i in 1:length(Trescuers))
if (Prescuers[i]==1) {
Trescuershist <- c(Trescuershist, Trescuers[i])
}
hist(Trescuershist, freq=FALSE, breaks=50, col="orange", border=FALSE, main=paste("so =", sovals[iso]))
xx <- density(Trescuershist)$x;
points(dexp(0:500, abs(rovals[iso])), col="black", type="l", lwd=2)
}
dev.off()
Prescuers <- rescuersL[[nsovals]][["P"]]
NMrescuers <- rescuersL[[nsovals]][["NM"]]
rescued <- NMrescuers[Prescuers == 1]
res <- c()
for (i in 0:max(rescued)) {
res <- c(res, length(rescued[rescued==i]))
}
print(res)
print("Simulation Stop")
return(res)
}
Rprint <- function (verbose, message)
{
if(verbose) {
print(message)
}
}
Rtrace <- function(matrice)
{
return(sum(diag(matrice)))
}
# building M & S as random matrices
# the distance to optimum xo
# the set of nba mutation phenotypic effects drawn into MVN(0,M)
# it is an K-allele-model with very large K = 2000
Rbuiltlandscape <- function(n, ne, Es, so, na)
{
if(ne >= n) m <- 1000 else m <- round(2 * n * ne / (n-ne) )
Xs <- matrix(rnorm(n*m),ncol=m)
S1 <- (1/m) * Xs %*% t(Xs)
Xm <- matrix(rnorm(n*m), ncol=m)
M1 <- (1/m) * Xm %*% t(Xm)
#scaling such that E(s)=bmax*tr(SM)/2-d
xx <- sqrt(2 * Es / Rtrace(S1 %*% M1) )
M <- xx * M1
S <- xx * S1
# matrix of allelic effects
A <- mvrnorm(na, numeric(n), M);
xo <- numeric(n)
if (so != 0) { #if so = 0, all the genotypic values stay at 0
x1 <- rnorm(n)
xo <- x1 * sqrt(2 * so / (t(x1) %*% S %*% x1) )
}
return(list(M=M, Sb=S, A=A, xo=xo))
}
\name{hello}
\alias{hello}
\title{Hello, World!}
\usage{
hello()
}
\description{
Prints 'Hello, world!'.
}
\examples{
hello()
}
\name{rcpp_hello}
\alias{rcpp_hello}
\title{Hello, Rcpp!}
\usage{
rcpp_hello()
}
\description{
Returns an \R \code{list} containing the character vector
\code{c("foo", "bar")} and the numeric vector \code{c(0, 1)}.
}
\examples{
rcpp_hello()
}
//
// Created by jimmy on 02/10/17.
//
#include "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(NumericMatrix mat) {
m_matrix = gsl_matrix_alloc(mat.rows(), mat.cols());
unsigned long x = 0;
for (unsigned long i = 0; i < mat.rows(); i++) {
for (unsigned long j = 0; j < mat.cols(); j++) {
setValue(i, j, mat[x]);
x += 1;
}
}
}
GMatrix::GMatrix(NumericVector vec) {
m_matrix = gsl_matrix_alloc(vec.size(), 1);
for (unsigned long i = 0; i < vec.size(); i++) {
setValue(i, 0, vec[i]);
}
}
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;
}
unsigned long GMatrix::getRowSize() const {
return m_matrix->size1;
}
unsigned long GMatrix::getColSize() const {
return m_matrix->size2;
}
//
// Created by jimmy on 02/10/17.
//
/**
* Class GMatrix
*
* Cette classe est utilisé pour encapsuler une matrice gsl
*
*/
#ifndef POPULATION_INPUTMATRIX_H
#define POPULATION_INPUTMATRIX_H
#include <string>
#include<gsl/gsl_matrix.h>
#include "Parameters.h"
#include <Rcpp.h>
using namespace Rcpp;
class GMatrix {
public:
GMatrix(unsigned long rows, unsigned long cols);
GMatrix(std::string file, unsigned long rows, unsigned long cols, bool reverse);
GMatrix(NumericMatrix mat);
GMatrix(NumericVector vec);
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.
//
#include "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;
}
double GRandom::exponential(double mu) {
return gsl_ran_exponential(m_rgsl, 1.0 / mu);
}
double GRandom::flat(double a, double b) {
return gsl_ran_flat(m_rgsl, a, b);
}
unsigned long GRandom::poisson(double mu) {
return gsl_ran_poisson(m_rgsl, mu);
}
//
// Created by jimmy on 03/10/17.
//
/**
* Class GRandom
*
* Cette classe est utilisé pour encapsuler le random de gsl
*
*/
#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.
//
#include "InputParser.h"
#include <algorithm>
InputParser::InputParser(int &argc, char **argv) {
for (int i=1; i < argc; ++i)
m_tokens.push_back(std::string(argv[i]));
}
const std::string &InputParser::getCmdOption(const std::string &option) const {
std::vector<std::string>::const_iterator itr;
itr = std::find(m_tokens.begin(), m_tokens.end(), option);
if (itr != m_tokens.end() && ++itr != m_tokens.end()){
return *itr;
}
static const std::string empty_string;
return empty_string;
}
bool InputParser::cmdOptionExists(const std::string &option) const {
return std::find(m_tokens.begin(), m_tokens.end(), option)
!= m_tokens.end();
}
//
// Created by jimmy on 05/10/17.
//
/**
* Class InputParser
*
* Cette classe est utilisé pour parser les paramètres d'entrées qui sont dans un fichier texte
*
*/
#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
PKG_CXXFLAGS= -DCOMPATIBILITYRCPP
## This is a C++11 package
CXX_STD = CXX11