From c1132606169875be6d07b54b30e8ae9446341bc2 Mon Sep 17 00:00:00 2001
From: Peter Carbonetto
Date: Sun, 4 Jun 2017 12:06:36 -0500
Subject: Removed FORCE_FLOAT from prdt.h/prdt.cpp.
---
src/eigenlib.cpp | 94 +++++++++----------
src/eigenlib.h | 9 +-
src/logistic.cpp | 61 ++++++-------
src/logistic.h | 93 +++++++++----------
src/prdt.cpp | 274 +++++++++++++++++++++++++++----------------------------
src/prdt.h | 23 ++---
6 files changed, 274 insertions(+), 280 deletions(-)
(limited to 'src')
diff --git a/src/eigenlib.cpp b/src/eigenlib.cpp
index 14ffbf1..7ad250f 100644
--- a/src/eigenlib.cpp
+++ b/src/eigenlib.cpp
@@ -1,6 +1,6 @@
/*
- Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011 Xiang Zhou
+ Genome-wide Efficient Mixed Model Association (GEMMA)
+ Copyright (C) 2011-2017, Xiang Zhou
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
@@ -13,7 +13,7 @@
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 .
+ along with this program. If not, see .
*/
#include
@@ -27,19 +27,23 @@
using namespace std;
using namespace Eigen;
-
-
-//on two different clusters, compare eigen vs lapack/gsl
-//dgemm, 5x or 0.5x faster or slower than lapack, 5x or 4x faster than gsl
-//dgemv, 20x or 4x faster than gsl,
-//eigen, 1x or 0.3x slower than lapack
-//invert, 20x or 10x faster than lapack
-
-void eigenlib_dgemm (const char *TransA, const char *TransB, const double alpha, const gsl_matrix *A, const gsl_matrix *B, const double beta, gsl_matrix *C)
-{
- Map, 0, OuterStride > A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
- Map, 0, OuterStride > B_mat(B->data, B->size1, B->size2, OuterStride(B->tda) );
- Map, 0, OuterStride > C_mat(C->data, C->size1, C->size2, OuterStride(C->tda) );
+// On two different clusters, compare eigen vs lapack/gsl:
+//
+// dgemm, 5x or 0.5x faster or slower than lapack, 5x or 4x faster than gsl
+// dgemv, 20x or 4x faster than gsl,
+// eigen, 1x or 0.3x slower than lapack
+// invert, 20x or 10x faster than lapack
+//
+void eigenlib_dgemm (const char *TransA, const char *TransB,
+ const double alpha, const gsl_matrix *A,
+ const gsl_matrix *B, const double beta,
+ gsl_matrix *C) {
+ Map, 0, OuterStride >
+ A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
+ Map, 0, OuterStride >
+ B_mat(B->data, B->size1, B->size2, OuterStride(B->tda) );
+ Map, 0, OuterStride >
+ C_mat(C->data, C->size1, C->size2, OuterStride(C->tda) );
if (*TransA=='N' || *TransA=='n') {
if (*TransB=='N' || *TransB=='n') {
@@ -55,19 +59,18 @@ void eigenlib_dgemm (const char *TransA, const char *TransB, const double alpha,
}
}
- //gsl_matrix_view C_view = gsl_matrix_view_array (C_mat.data(), C->size1, C->size2);
- //gsl_matrix_memcpy (C, &C_view.matrix);
-
return;
}
-
-
-void eigenlib_dgemv (const char *TransA, const double alpha, const gsl_matrix *A, const gsl_vector *x, const double beta, gsl_vector *y)
-{
- Map, 0, OuterStride > A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
- Map, 0, InnerStride > x_vec(x->data, x->size, InnerStride(x->stride) );
- Map, 0, InnerStride > y_vec(y->data, y->size, InnerStride(y->stride) );
+void eigenlib_dgemv (const char *TransA, const double alpha,
+ const gsl_matrix *A, const gsl_vector *x,
+ const double beta, gsl_vector *y) {
+ Map, 0, OuterStride >
+ A_mat(A->data, A->size1, A->size2, OuterStride(A->tda) );
+ Map, 0, InnerStride >
+ x_vec(x->data, x->size, InnerStride(x->stride) );
+ Map, 0, InnerStride >
+ y_vec(y->data, y->size, InnerStride(y->stride) );
if (*TransA=='N' || *TransA=='n') {
y_vec=alpha*A_mat*x_vec+beta*y_vec;
@@ -78,38 +81,35 @@ void eigenlib_dgemv (const char *TransA, const double alpha, const gsl_matrix *A
return;
}
-
-
-void eigenlib_invert(gsl_matrix *A)
-{
- Map > A_mat(A->data, A->size1, A->size2);
+void eigenlib_invert(gsl_matrix *A) {
+ Map >
+ A_mat(A->data, A->size1, A->size2);
A_mat=A_mat.inverse();
return;
}
-
-void eigenlib_dsyr (const double alpha, const gsl_vector *b, gsl_matrix *A)
-{
- Map > A_mat(A->data, A->size1, A->size2);
- Map, 0, OuterStride > b_vec(b->data, b->size, OuterStride(b->stride) );
-
+void eigenlib_dsyr (const double alpha, const gsl_vector *b, gsl_matrix *A) {
+ Map >
+ A_mat(A->data, A->size1, A->size2);
+ Map, 0, OuterStride >
+ b_vec(b->data, b->size, OuterStride(b->stride) );
A_mat=alpha*b_vec*b_vec.transpose()+A_mat;
-
return;
}
-
-void eigenlib_eigensymm (const gsl_matrix *G, gsl_matrix *U, gsl_vector *eval)
-{
- Map, 0, OuterStride > G_mat(G->data, G->size1, G->size2, OuterStride(G->tda) );
- Map, 0, OuterStride > U_mat(U->data, U->size1, U->size2, OuterStride(U->tda) );
- Map, 0, OuterStride > eval_vec(eval->data, eval->size, OuterStride(eval->stride) );
+void eigenlib_eigensymm (const gsl_matrix *G, gsl_matrix *U,
+ gsl_vector *eval) {
+ Map, 0, OuterStride >
+ G_mat(G->data, G->size1, G->size2, OuterStride(G->tda) );
+ Map, 0, OuterStride >
+ U_mat(U->data, U->size1, U->size2, OuterStride(U->tda) );
+ Map, 0, OuterStride >
+ eval_vec(eval->data, eval->size, OuterStride(eval->stride) );
SelfAdjointEigenSolver es(G_mat);
- if (es.info() != Success) abort();
-
+ if (es.info() != Success)
+ abort();
eval_vec=es.eigenvalues();
U_mat=es.eigenvectors();
-
return;
}
diff --git a/src/eigenlib.h b/src/eigenlib.h
index f869786..8cb8880 100644
--- a/src/eigenlib.h
+++ b/src/eigenlib.h
@@ -23,8 +23,13 @@
using namespace std;
-void eigenlib_dgemm (const char *TransA, const char *TransB, const double alpha, const gsl_matrix *A, const gsl_matrix *B, const double beta, gsl_matrix *C);
-void eigenlib_dgemv (const char *TransA, const double alpha, const gsl_matrix *A, const gsl_vector *x, const double beta, gsl_vector *y);
+void eigenlib_dgemm (const char *TransA, const char *TransB,
+ const double alpha, const gsl_matrix *A,
+ const gsl_matrix *B, const double beta,
+ gsl_matrix *C);
+void eigenlib_dgemv (const char *TransA, const double alpha,
+ const gsl_matrix *A, const gsl_vector *x,
+ const double beta, gsl_vector *y);
void eigenlib_invert(gsl_matrix *A);
void eigenlib_dsyr (const double alpha, const gsl_vector *b, gsl_matrix *A);
void eigenlib_eigensymm (const gsl_matrix *G, gsl_matrix *U, gsl_vector *eval);
diff --git a/src/logistic.cpp b/src/logistic.cpp
index 1b47946..002ce98 100644
--- a/src/logistic.cpp
+++ b/src/logistic.cpp
@@ -7,45 +7,40 @@
#include
#include "logistic.h"
-// I need to bundle all the data that goes to the function to optimze together.
+// I need to bundle all the data that goes to the function to optimze
+// together.
typedef struct{
gsl_matrix_int *X;
gsl_vector_int *nlev;
gsl_vector *y;
- gsl_matrix *Xc; // continuous covariates Matrix Nobs x Kc (NULL if not used)
+ gsl_matrix *Xc; // continuous covariates matrix Nobs x Kc (NULL if not used)
double lambdaL1;
double lambdaL2;
-}fix_parm_mixed_T;
-
-
-
-
-
-
-double fLogit_mixed(gsl_vector *beta
- ,gsl_matrix_int *X
- ,gsl_vector_int *nlev
- ,gsl_matrix *Xc
- ,gsl_vector *y
- ,double lambdaL1
- ,double lambdaL2)
-{
+} fix_parm_mixed_T;
+
+double fLogit_mixed(gsl_vector *beta,
+ gsl_matrix_int *X,
+ gsl_vector_int *nlev,
+ gsl_matrix *Xc,
+ gsl_vector *y,
+ double lambdaL1,
+ double lambdaL2) {
int n = y->size;
- // int k = X->size2;
int npar = beta->size;
double total = 0;
double aux = 0;
- /* omp_set_num_threads(ompthr); */
- /* /\* Changed loop start at 1 instead of 0 to avoid regularization of beta_0*\/ */
- /* /\*#pragma omp parallel for reduction (+:total)*\/ */
+ // Changed loop start at 1 instead of 0 to avoid regularization of
+ // beta_0*\/ */
+ // #pragma omp parallel for reduction (+:total)
for(int i = 1; i < npar; ++i)
total += beta->data[i]*beta->data[i];
total = (-total*lambdaL2/2);
- /* /\*#pragma omp parallel for reduction (+:aux)*\/ */
+ // #pragma omp parallel for reduction (+:aux)
for(int i = 1; i < npar; ++i)
aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]);
total = total-aux*lambdaL1;
- /* #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) reduction (+:total) */
+ // #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y)
+ // #reduction (+:total)
for(int i = 0; i < n; ++i) {
double Xbetai=beta->data[0];
int iParm=1;
@@ -94,11 +89,12 @@ wgsl_mixed_optim_df (const gsl_vector *beta, void *params,
int n = p->y->size;
int K = p->X->size2;
int Kc = p->Xc->size2;
- int npar = beta->size;
+ int npar = beta->size;
+
// Intitialize gradient out necessary?
for(int i = 0; i < npar; ++i)
out->data[i]= 0;
- /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0 */
+ // Changed loop start at 1 instead of 0 to avoid regularization of beta 0.
for(int i = 1; i < npar; ++i)
out->data[i]= p->lambdaL2*beta->data[i];
for(int i = 1; i < npar; ++i)
@@ -113,7 +109,8 @@ wgsl_mixed_optim_df (const gsl_vector *beta, void *params,
Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm];
iParm+=p->nlev->data[k]-1;
}
- // Adding the continuous
+
+ // Adding the continuous.
for(int k = 0; k < Kc; ++k)
Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm++];
@@ -126,7 +123,8 @@ wgsl_mixed_optim_df (const gsl_vector *beta, void *params,
out->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]+=pn;
iParm+=p->nlev->data[k]-1;
}
- // Adding the continuous
+
+ // Adding the continuous.
for(int k = 0; k < Kc; ++k) {
out->data[iParm++] += gsl_matrix_get(p->Xc,i,k)*pn;
}
@@ -134,12 +132,9 @@ wgsl_mixed_optim_df (const gsl_vector *beta, void *params,
}
-
-/* The Hessian of f */
-void
-wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params,
- gsl_matrix *out)
-{
+// The Hessian of f.
+void wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params,
+ gsl_matrix *out) {
fix_parm_mixed_T *p = (fix_parm_mixed_T *)params;
int n = p->y->size;
int K = p->X->size2;
diff --git a/src/logistic.h b/src/logistic.h
index a68ee09..e951935 100644
--- a/src/logistic.h
+++ b/src/logistic.h
@@ -1,52 +1,54 @@
#ifndef LOGISTIC_H_ /* Include guard */
#define LOGISTIC_H_
-/* Mixed interface */
-void logistic_mixed_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc
- ,gsl_matrix_int *X //Matrix Nobs x K
- ,gsl_vector_int *nlev // Vector with number categories
- ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc
- ,gsl_vector *yhat //Vector of prob. predicted by the logistic
- );
+// Mixed interface.
+void logistic_mixed_pred(gsl_vector *beta, // Vector of parameters
+ // length = 1+Sum_k(C_k-1)+Kc.
+ gsl_matrix_int *X, // Matrix Nobs x K.
+ gsl_vector_int *nlev, // Vector with num. categories.
+ gsl_matrix *Xc, // Continuous covariates matrix
+ // Nobs x Kc
+ gsl_vector *yhat); // Vector of prob. predicted by
+ // the logistic.
-int logistic_mixed_fit(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc
- ,gsl_matrix_int *X //Matrix Nobs x K
- ,gsl_vector_int *nlev // Vector with number categories
- ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc
- ,gsl_vector *y //Vector of prob. to predict
- ,double lambdaL1 // Regularization L1 0.0 if not used
- ,double lambdaL2); // Regularization L2 0.0 if not used
+int logistic_mixed_fit(gsl_vector *beta, // Vector of parameters
+ // length = 1+Sum_k(C_k-1)+Kc
+ gsl_matrix_int *X, // Matrix Nobs x K.
+ gsl_vector_int *nlev, // Vector with number categories.
+ gsl_matrix *Xc, // Continuous covariates
+ // matrix Nobs x Kc
+ gsl_vector *y, // Vector of prob. to predict.
+ double lambdaL1, // Reg. L1 0.0 if not used.
+ double lambdaL2); // Reg. L2 0.0 if not used.
-double fLogit_mixed(gsl_vector *beta
- ,gsl_matrix_int *X
- ,gsl_vector_int *nlev
- ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc
- ,gsl_vector *y
- ,double lambdaL1
- ,double lambdaL2);
+double fLogit_mixed(gsl_vector *beta,
+ gsl_matrix_int *X,
+ gsl_vector_int *nlev,
+ gsl_matrix *Xc, // continuous covariates matrix Nobs x Kc
+ gsl_vector *y,
+ double lambdaL1,
+ double lambdaL2);
-/* Categorical only interface */
-void logistic_cat_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc
- ,gsl_matrix_int *X //Matrix Nobs x K
- ,gsl_vector_int *nlev // Vector with number categories
- ,gsl_vector *yhat //Vector of prob. predicted by the logistic
- );
+// Categorical-only interface.
+void logistic_cat_pred(gsl_vector *beta, // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc
+ gsl_matrix_int *X, //Matrix Nobs x K
+ gsl_vector_int *nlev, // Vector with number categories
+ gsl_vector *yhat); //Vector of prob. predicted by the logistic
-int logistic_cat_fit(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc
- ,gsl_matrix_int *X //Matrix Nobs x K
- ,gsl_vector_int *nlev // Vector with number categories
- ,gsl_vector *y //Vector of prob. to predict
- ,double lambdaL1 // Regularization L1 0.0 if not used
- ,double lambdaL2); // Regularization L2 0.0 if not used
-
-double fLogit_cat(gsl_vector *beta
- ,gsl_matrix_int *X
- ,gsl_vector_int *nlev
- ,gsl_vector *y
- ,double lambdaL1
- ,double lambdaL2);
+int logistic_cat_fit(gsl_vector *beta, // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc
+ gsl_matrix_int *X, //Matrix Nobs x K
+ gsl_vector_int *nlev, // Vector with number categories
+ gsl_vector *y, //Vector of prob. to predict
+ double lambdaL1, // Regularization L1 0.0 if not used
+ double lambdaL2); // Regularization L2 0.0 if not used
+double fLogit_cat(gsl_vector *beta,
+ gsl_matrix_int *X,
+ gsl_vector_int *nlev,
+ gsl_vector *y,
+ double lambdaL1,
+ double lambdaL2);
/* Continuous only interface */
void logistic_cont_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) + Kc
@@ -60,11 +62,10 @@ int logistic_cont_fit(gsl_vector *beta // Vector of parameters length = 1 + Sum
,double lambdaL1 // Regularization L1 0.0 if not used
,double lambdaL2); // Regularization L2 0.0 if not used
-double fLogit_cont(gsl_vector *beta
- ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc
- ,gsl_vector *y
- ,double lambdaL1
- ,double lambdaL2);
-
+double fLogit_cont(gsl_vector *beta,
+ gsl_matrix *Xc, // Continuous covariates matrix Nobs x Kc .
+ gsl_vector *y,
+ double lambdaL1,
+ double lambdaL2);
#endif // LOGISTIC_H_
diff --git a/src/prdt.cpp b/src/prdt.cpp
index 2875119..db0fa14 100644
--- a/src/prdt.cpp
+++ b/src/prdt.cpp
@@ -1,6 +1,6 @@
/*
Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011 Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
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
@@ -13,10 +13,8 @@
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 .
- */
-
-
+ along with this program. If not, see .
+*/
#include
#include
@@ -33,28 +31,16 @@
#include "gsl/gsl_linalg.h"
#include "gsl/gsl_blas.h"
-
#include "io.h"
-#include "lapack.h" //for functions EigenDecomp
+#include "lapack.h"
#include "gzstream.h"
-
-#ifdef FORCE_FLOAT
-#include "io_float.h"
-#include "prdt_float.h"
-#include "mathfunc_float.h"
-#else
#include "io.h"
#include "prdt.h"
#include "mathfunc.h"
-#endif
using namespace std;
-
-
-
-void PRDT::CopyFromParam (PARAM &cPar)
-{
+void PRDT::CopyFromParam (PARAM &cPar) {
a_mode=cPar.a_mode;
d_pace=cPar.d_pace;
@@ -81,19 +67,14 @@ void PRDT::CopyFromParam (PARAM &cPar)
return;
}
-void PRDT::CopyToParam (PARAM &cPar)
-{
+void PRDT::CopyToParam (PARAM &cPar) {
cPar.ns_test=ns_test;
cPar.time_eigen=time_eigen;
return;
}
-
-
-
-void PRDT::WriteFiles (gsl_vector *y_prdt)
-{
+void PRDT::WriteFiles (gsl_vector *y_prdt) {
string file_str;
file_str=path_out+"/"+file_out;
file_str+=".";
@@ -101,7 +82,10 @@ void PRDT::WriteFiles (gsl_vector *y_prdt)
file_str+=".txt";
ofstream outfile (file_str.c_str(), ofstream::out);
- if (!outfile) {cout<<"error writing file: "<size2; j++) {
- outfile<size, ni_total=G->size1;
gsl_matrix *Goo=gsl_matrix_alloc (ni_test, ni_test);
@@ -190,7 +172,9 @@ void PRDT::AddBV (gsl_matrix *G, const gsl_vector *u_hat, gsl_vector *y_prdt)
clock_t time_start=clock();
EigenDecomp (Goo, U, eval, 0);
for (size_t i=0; isize; i++) {
- if (gsl_vector_get(eval,i)<1e-10) {gsl_vector_set(eval, i, 0);}
+ if (gsl_vector_get(eval,i)<1e-10) {
+ gsl_vector_set(eval, i, 0);
+ }
}
time_eigen=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
@@ -198,12 +182,15 @@ void PRDT::AddBV (gsl_matrix *G, const gsl_vector *u_hat, gsl_vector *y_prdt)
gsl_blas_dgemv (CblasTrans, 1.0, U, u_hat, 0.0, Utu);
for (size_t i=0; isize; i++) {
d=gsl_vector_get(eval, i);
- if (d!=0) {d=gsl_vector_get(Utu, i)/d; gsl_vector_set(Utu, i, d);}
+ if (d!=0) {
+ d=gsl_vector_get(Utu, i)/d;
+ gsl_vector_set(Utu, i, d);
+ }
}
gsl_blas_dgemv (CblasNoTrans, 1.0, U, Utu, 0.0, eval);
gsl_blas_dgemv (CblasNoTrans, 1.0, Gfo, eval, 1.0, y_prdt);
- //free matrices
+ // Free matrices.
gsl_matrix_free(Goo);
gsl_matrix_free(Gfo);
gsl_matrix_free(U);
@@ -215,13 +202,12 @@ void PRDT::AddBV (gsl_matrix *G, const gsl_vector *u_hat, gsl_vector *y_prdt)
return;
}
-
-
-void PRDT::AnalyzeBimbam (gsl_vector *y_prdt)
-{
+void PRDT::AnalyzeBimbam (gsl_vector *y_prdt) {
igzstream infile (file_geno.c_str(), igzstream::in);
-// ifstream infile (file_geno.c_str(), ifstream::in);
- if (!infile) {cout<<"error reading genotype file:"<size==n_miss) {cout<<"snp "<size==n_miss) {
+ cout << "snp " << rs << " has missing genotype for all " <<
+ "individuals and will be ignored." << endl;
+ continue;}
+
x_mean/=(double)(x->size-n_miss);
x_train_mean/=(double)(n_train_nomiss);
@@ -303,17 +305,13 @@ void PRDT::AnalyzeBimbam (gsl_vector *y_prdt)
return;
}
-
-
-
-
-
-
-void PRDT::AnalyzePlink (gsl_vector *y_prdt)
-{
+void PRDT::AnalyzePlink (gsl_vector *y_prdt) {
string file_bed=file_bfile+".bed";
ifstream infile (file_bed.c_str(), ios::binary);
- if (!infile) {cout<<"error reading bed file:"< b;
@@ -324,11 +322,11 @@ void PRDT::AnalyzePlink (gsl_vector *y_prdt)
gsl_vector *x=gsl_vector_alloc (y_prdt->size);
- //calculate n_bit and c, the number of bit for each snp
+ // Calculate n_bit and c, the number of bit for each SNP.
if (indicator_idv.size()%4==0) {n_bit=indicator_idv.size()/4;}
else {n_bit=indicator_idv.size()/4+1; }
- //print the first three majic numbers
+ // Print the first 3 magic numbers.
for (size_t i=0; i<3; ++i) {
infile.read(ch,1);
b=ch[0];
@@ -337,39 +335,71 @@ void PRDT::AnalyzePlink (gsl_vector *y_prdt)
ns_test=0;
for (vector::size_type t=0; tsize==n_miss) {cout<<"snp "<size==n_miss) {
+ cout << "snp " << rs << " has missing genotype for all " <<
+ "individuals and will be ignored."<size-n_miss);
x_train_mean/=(double)(n_train_nomiss);
@@ -407,13 +441,10 @@ void PRDT::AnalyzePlink (gsl_vector *y_prdt)
return;
}
-
-
-
-//predict missing phenotypes using ridge regression
-//Y_hat contains fixed effects
-void PRDT::MvnormPrdt (const gsl_matrix *Y_hat, const gsl_matrix *H, gsl_matrix *Y_full)
-{
+// Predict missing phenotypes using ridge regression.
+// Y_hat contains fixed effects
+void PRDT::MvnormPrdt (const gsl_matrix *Y_hat, const gsl_matrix *H,
+ gsl_matrix *Y_full) {
gsl_vector *y_obs=gsl_vector_alloc (np_obs);
gsl_vector *y_miss=gsl_vector_alloc (np_miss);
gsl_matrix *H_oo=gsl_matrix_alloc (np_obs, np_obs);
@@ -422,20 +453,22 @@ void PRDT::MvnormPrdt (const gsl_matrix *Y_hat, const gsl_matrix *H, gsl_matrix
size_t c_obs1=0, c_obs2=0, c_miss1=0, c_miss2=0;
- //obtain H_oo, H_mo
+ // Obtain H_oo, H_mo.
c_obs1=0; c_miss1=0;
for (vector::size_type i1=0; i1::size_type j1=0; j1::size_type i2=0; i2::size_type i2=0;
+ i2::size_type j2=0; j2::size_type j2=0;
+ j2::size_type i=0; i::size_type i=0;
+ i::size_type j=0; j::size_type i=0; i::size_type i=0;
+ i::size_type j=0; j::size_type i=0; i::size_type j=0; j<2; ++j) {
- if (indicator_pheno[i][j]==1) {
- gsl_vector_set (y_obs, c_obs1, gsl_matrix_get (Y_full, i, j+k*2)-gsl_matrix_get (Y_hat, i, j) );
- c_obs1++;
- } else {
- gsl_vector_set (y_miss, c_miss1, gsl_matrix_get (Y_hat, i, j) );
- c_miss1++;
- }
- }
- }
-
- LUSolve (H_oo, pmt, y_obs, Hiy);
-
- gsl_blas_dgemv (CblasNoTrans, 1.0, H_mo, Hiy, 1.0, y_miss);
-
- //put back predicted y_miss to Y_full
- c_miss1=0;
- for (vector::size_type i=0; i::size_type j=0; j<2; ++j) {
- if (indicator_pheno[i][j]==0) {
- gsl_matrix_set (Y_full, i, j+k*2, gsl_vector_get (y_miss, c_miss1) );
- c_miss1++;
- }
- }
- }
- }
- }
-*/
- //free matrices
+
+ // Free matrices.
gsl_vector_free(y_obs);
gsl_vector_free(y_miss);
gsl_matrix_free(H_oo);
diff --git a/src/prdt.h b/src/prdt.h
index 8af2cee..2da9fd0 100644
--- a/src/prdt.h
+++ b/src/prdt.h
@@ -1,6 +1,6 @@
/*
- Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011 Xiang Zhou
+ Genome-wide Efficient Mixed Model Association (GEMMA)
+ Copyright (C) 2011-2017, Xiang Zhou
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
@@ -13,31 +13,25 @@
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 .
+ along with this program. If not, see .
*/
#ifndef __PRDT_H__
#define __PRDT_H__
-
#include
#include