Commit 2fdd5d13 authored by jlopez's avatar jlopez
Browse files

Replace GSL lib

parent 2961dd1e
......@@ -7,22 +7,31 @@
#include <fstream>
#include <sstream>
#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
//#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);
/*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);
}
}*/
m_sizer = rows;
m_sizec = cols;
m2_matrix = new double*[rows];
for (unsigned long i = 0; i < rows; i++) {
m2_matrix[i] = new double[cols];
for (unsigned long j = 0; j < cols; j++) {
setValue(i, j, 0.0);
}
}
}
GMatrix::GMatrix(NumericMatrix mat) {
m_matrix = gsl_matrix_alloc(mat.rows(), mat.cols());
/*m_matrix = gsl_matrix_alloc(mat.rows(), mat.cols());
unsigned long x = 0;
......@@ -32,18 +41,38 @@ GMatrix::GMatrix(NumericMatrix mat) {
x += 1;
}
}*/
m_sizer = mat.rows();
m_sizec = mat.cols();
unsigned long x = 0;
m2_matrix = new double*[m_sizer];
for (unsigned long i = 0; i < m_sizer; i++) {
m2_matrix[i] = new double[m_sizec];
for (unsigned long j = 0; j < m_sizec; j++) {
setValue(i, j, mat[x]);
x += 1;
}
}
}
GMatrix::GMatrix(NumericVector vec) {
m_matrix = gsl_matrix_alloc(vec.size(), 1);
/*m_matrix = gsl_matrix_alloc(vec.size(), 1);
for (unsigned long i = 0; i < vec.size(); i++) {
setValue(i, 0, vec[i]);
}*/
m_sizer = vec.size();
m_sizec = 1;
m2_matrix = new double*[m_sizer];
for (unsigned long i = 0; i < m_sizer; i++) {
m2_matrix[i] = new double[m_sizec];
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) {
......@@ -80,84 +109,156 @@ GMatrix::GMatrix(std::string file, unsigned long rows, unsigned long cols, bool
}
infile.close();
}
}*/
m_sizer = rows;
m_sizec = cols;
m2_matrix = new double*[rows];
for (unsigned long i = 0; i < rows; i++) {
m2_matrix[i] = new double[cols];
for (unsigned long j = 0; j < cols; j++) {
setValue(i, j, 0.0);
}
}
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
if(reverse) {
setValue(j, i, value);
//gsl_matrix_set(m_matrix, j, i, value);
} else {
setValue(i, j, value);
//gsl_matrix_set(m_matrix, i, j, value);
}
j++;
}
i++;
}
infile.close();
}
GMatrix::~GMatrix() {
gsl_matrix_free(m_matrix);
//gsl_matrix_free(m_matrix);
for(int i = 0; i < m_sizer; ++i) {
delete [] m2_matrix[i];
}
delete [] m2_matrix;
}
double GMatrix::getValue(unsigned long i, unsigned long j) const {
return gsl_matrix_get(m_matrix, i, j);
//return gsl_matrix_get(m_matrix, i, j);
return m2_matrix[i][j];
}
void GMatrix::setValue(unsigned long i, unsigned long j, double value) {
gsl_matrix_set(m_matrix, i, j, value);
//gsl_matrix_set(m_matrix, i, j, value);
m2_matrix[i][j] = value;
}
gsl_matrix *GMatrix::getMatrix() const {
/*gsl_matrix *GMatrix::getMatrix() const {
return m_matrix;
}
}*/
void GMatrix::sub(const GMatrix *gMatrix) {
gsl_matrix_sub(m_matrix, gMatrix->getMatrix());
//gsl_matrix_sub(m_matrix, gMatrix->getMatrix());
for (unsigned long i = 0; i < m_sizer; i++) {
for (unsigned long j = 0; j < m_sizec; j++) {
m2_matrix[i][j] -= gMatrix->getValue(i, j);
}
}
}
double GMatrix::trace() const {
gsl_vector_view v = gsl_matrix_diagonal(m_matrix);
//gsl_vector_view v = gsl_matrix_diagonal(m_matrix);
unsigned long size = std::min(m_matrix->size1, m_matrix->size2);
//unsigned long size = std::min(m_matrix->size1, m_matrix->size2);
unsigned long size = std::min(m_sizer, m_sizec);
double sum = 0.0;
for (unsigned long i = 0; i < size; i++) {
sum += gsl_vector_get(&v.vector, i);
//sum += gsl_vector_get(&v.vector, i);
sum += getValue(i, i);
}
return sum;
}
void GMatrix::transpose() const {
gsl_matrix_transpose(m_matrix);
GMatrix GMatrix::getTranspose() const {
//gsl_matrix_transpose(m_matrix);
GMatrix matrix(m_sizec, m_sizer);
for (unsigned long i = 0; i < m_sizer; i++) {
for (unsigned long j = 0; j < m_sizec; j++) {
matrix.setValue(j, i, m2_matrix[i][j]);
}
}
return matrix;
}
GMatrix* GMatrix::mult(double val) {
for (unsigned long i = 0; i < m_matrix->size1; ++i) {
/*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);
}
}*/
for (unsigned long i = 0; i < m_sizer; ++i) {
for (unsigned long j = 0; j < m_sizec; ++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 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);
}
}*/
for (unsigned long i = 0; i < m_sizer; ++i) {
for (unsigned long j = 0; j < m_sizec; ++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 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);
}
}*/
for (unsigned long i = 0; i < m_sizer; ++i) {
for (unsigned long j = 0; j < m_sizec; ++j) {
double v = std::sqrt(getValue(i, j));
setValue(i, j, v);
}
}
return this;
}
unsigned long GMatrix::getRowSize() const {
return m_matrix->size1;
//return m_matrix->size1;
return m_sizer;
}
unsigned long GMatrix::getColSize() const {
return m_matrix->size2;
//return m_matrix->size2;
return m_sizec;
}
......
......@@ -13,7 +13,7 @@
#define POPULATION_INPUTMATRIX_H
#include <string>
#include<gsl/gsl_matrix.h>
//#include<gsl/gsl_matrix.h>
#include "Parameters.h"
#include <Rcpp.h>
......@@ -28,13 +28,13 @@ public:
GMatrix(NumericVector vec);
virtual ~GMatrix();
gsl_matrix *getMatrix() const;
//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 getTranspose() const;
GMatrix* mult(double val);
GMatrix* div(double val);
GMatrix* sqrt();
......@@ -42,7 +42,10 @@ public:
unsigned long getColSize() const;
private:
gsl_matrix* m_matrix;
//gsl_matrix* m_matrix;
double ** m2_matrix;
unsigned long m_sizer;
unsigned long m_sizec;
};
......
......@@ -7,27 +7,33 @@
#include <time.h>
GRandom::GRandom() {
const gsl_rng_type *Tgsl;
/*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));
gsl_rng_set (m_rgsl, (unsigned long) time (nullptr));*/
}
GRandom::~GRandom() {
delete m_rgsl;
//delete m_rgsl;
}
double GRandom::exponential(double mu) {
return gsl_ran_exponential(m_rgsl, 1.0 / mu);
//return gsl_ran_exponential(m_rgsl, 1.0 / mu);
std::exponential_distribution<double> distribution(mu);
return distribution(m_generator);
}
double GRandom::flat(double a, double b) {
return gsl_ran_flat(m_rgsl, a, b);
//return gsl_ran_flat(m_rgsl, a, b);
std::uniform_real_distribution<double> distribution(a,b);
return distribution(m_generator);
}
unsigned long GRandom::poisson(double mu) {
return gsl_ran_poisson(m_rgsl, mu);
//return gsl_ran_poisson(m_rgsl, mu);
std::poisson_distribution<int> distribution(mu);
return distribution(m_generator);
}
......@@ -12,8 +12,9 @@
#ifndef POPULATION_GRANDOM_H
#define POPULATION_GRANDOM_H
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
//#include <gsl/gsl_rng.h>
//#include <gsl/gsl_randist.h>
#include <random>
class GRandom {
......@@ -26,8 +27,8 @@ public:
virtual unsigned long poisson(double mu);
private:
gsl_rng* m_rgsl;
//gsl_rng* m_rgsl;
std::default_random_engine m_generator;
};
......
PKG_CXXFLAGS= -DCOMPATIBILITYRCPP
## This is a C++11 package
CXX_STD = CXX11
PKG_CPPFLAGS=-I$(LIB_GSL)/include -I../inst/include
PKG_LIBS=-L$(LIB_GSL)/lib -lgsl -lgslcblas $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
#PKG_CPPFLAGS=-I$(LIB_GSL)/include -I../inst/include
#PKG_LIBS=-L$(LIB_GSL)/lib -lgsl -lgslcblas $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
......@@ -4,5 +4,5 @@
PKG_CXXFLAGS= -DCOMPATIBILITYRCPP -fno-omit-frame-pointer
## This is a C++11 package
CXX_STD = CXX11
PKG_CPPFLAGS=-I$(LIB_GSL)/include -I../inst/include
PKG_LIBS=-L$(LIB_GSL)/lib -lgsl -lgslcblas $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
#PKG_CPPFLAGS=-I$(LIB_GSL)/include -I../inst/include
#PKG_LIBS=-L$(LIB_GSL)/lib -lgsl -lgslcblas $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
......@@ -14,7 +14,7 @@
#include <iostream>
#include <string>
#include <gsl/gsl_blas.h>
//#include <gsl/gsl_blas.h>
#include "GMatrix.h"
......@@ -27,11 +27,30 @@ public:
}
static void dgemmTN(GMatrix* A, GMatrix* B, GMatrix* C) {
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, A->getMatrix(), B->getMatrix(), 0.0, C->getMatrix());
//gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, A->getMatrix(), B->getMatrix(), 0.0, C->getMatrix());
GMatrix T = A->getTranspose();
for (unsigned long i = 0; i < C->getRowSize(); ++i) {
for (unsigned long j = 0; j < C->getColSize(); ++j) {
double val = 0.0;
for (unsigned long k = 0; k < T.getColSize(); ++k) {
val += T.getValue(i, k) * B->getValue(k, j);
}
C->setValue(i, j, val);
}
}
}
static void dgemmNN(GMatrix* A, GMatrix* B, GMatrix* C) {
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A->getMatrix(), B->getMatrix(), 0.0, C->getMatrix());
//gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A->getMatrix(), B->getMatrix(), 0.0, C->getMatrix());
for (unsigned long i = 0; i < C->getRowSize(); ++i) {
for (unsigned long j = 0; j < C->getColSize(); ++j) {
double val = 0.0;
for (unsigned long k = 0; k < A->getColSize(); ++k) {
val += A->getValue(i, k) * B->getValue(k, j);
}
C->setValue(i, j, val);
}
}
}
/**
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment