[Genabel-commits] r1059 - pkg/OmicABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 19 11:55:08 CET 2012
Author: dfabregat
Date: 2012-12-19 11:55:08 +0100 (Wed, 19 Dec 2012)
New Revision: 1059
Modified:
pkg/OmicABEL/src/CLAK_GWAS.c
pkg/OmicABEL/src/REML.c
pkg/OmicABEL/src/REML.h
pkg/OmicABEL/src/fgls_chol.c
pkg/OmicABEL/src/utils.c
Log:
[METHODOLOGY] Chol variant uses the REML based on the eigendecomposition
* REML.[ch]: removed every function and data structure related to chol
* CLAK_GWAS.c: chol uses eig for REML
* fgls_chol.c: Added a function to build SPD Phi out of eigenvectors and
fixed eigenvalues of Phi. A bit of a cleanup.
* utils.c updated accordingly to calculate ooc_b for chol
Modified: pkg/OmicABEL/src/CLAK_GWAS.c
===================================================================
--- pkg/OmicABEL/src/CLAK_GWAS.c 2012-12-18 15:31:36 UTC (rev 1058)
+++ pkg/OmicABEL/src/CLAK_GWAS.c 2012-12-19 10:55:08 UTC (rev 1059)
@@ -79,11 +79,11 @@
// var, nths, xtile, ytile, xb, yb, write_output
);
- if ( !strcmp( var, "chol" ) )
- {
- printf("Chol variant is broken. Please use -var eigen\n");
- exit(EXIT_FAILURE);
- }
+ /*if ( !strcmp( var, "chol" ) )*/
+ /*{*/
+ /*printf("Chol variant is broken. Please use -var eigen\n");*/
+ /*exit(EXIT_FAILURE);*/
+ /*}*/
printf("\n**************************************************************************\n");
printf("*** This is an alpha version under development. Use it at your own risk!!!\n");
@@ -116,7 +116,8 @@
printf( "Estimating GWAS parameters: heritability and variance... " );
fflush( stdout );
gettimeofday(&start, NULL);
- estimates_chol( &cf ); // Computation
+ /*estimates_chol( &cf ); // Computation*/
+ estimates_eig( &cf ); // Computation
gettimeofday(&end, NULL);
printf( "Done ( took %.3f secs )\n", get_diff_ms(&start, &end)/1000. );
fflush( stdout );
Modified: pkg/OmicABEL/src/REML.c
===================================================================
--- pkg/OmicABEL/src/REML.c 2012-12-18 15:31:36 UTC (rev 1058)
+++ pkg/OmicABEL/src/REML.c 2012-12-19 10:55:08 UTC (rev 1059)
@@ -42,219 +42,6 @@
static void eigen_mr3(int n, double *Phi, double *Z, double *W);
static void premultiplyXL( int n, int wXL, double *Z, double *XL, double *ZtXL );
-
-void estimates_chol( FGLS_config_t *cf )
-{
- estimates_chol_t data;
- int j, k;
-#if DEBUG
- FILE *f_ests;
-#endif
-
- // Dimensions
- data.n = cf->n;
- data.wXL = cf->wXL;
- // Load data
- data.Phi = cf->Phi;
- data.X = cf->XL;
-
- set_multi_threaded_BLAS(cf->num_threads);
- data.Y = fgls_malloc( data.n * sizeof(double) );
- data.beta = fgls_malloc( cf->wXL * sizeof(double) );
- for ( j = 0; j < cf->t; j++ )
- {
- // Should Y_fp
- sync_read( data.Y, cf->Y, cf->n, cf->n * j );
- // Sanity
- /*checkNoNans(cf->n, data.Y, "Y");*/
- average( data.Y, cf->n, 1, cf->threshold, "TRAIT", &cf->Y_fvi->fvi_data[(cf->n+j)*NAMELENGTH], NAMELENGTH, 1 );
- /*printf("var(Y): %.15e\n", variance(data.Y, cf->n));*/
-
- data.sigma = variance(data.Y, data.n);
-
- // [a,b] = [0,0.98], tol = 0.01
- minimize( 0.0, 0.99, 1e-8, &polygenic_REML_logLik_chol_wrapper, (void*)&data );
-#if 0 //DEBUG
- printf("h: %.15e - sigma: %.15e - res_sigma: %.15e\n", data.h, data.sigma, data.res_sigma);
-#endif
- cf->h2[j] = data.h;
- cf->sigma2[j] = data.sigma;
- cf->res_sigma2[j] = data.res_sigma;
- for ( k = 0; k < cf->wXL; k++ )
- cf->beta_ests[cf->wXL*j+k] = data.beta[k];
- }
-#if DEBUG
- // Store data - For testing purposes
-/* f_ests = fgls_fopen( cf->ests_path, "wb" );
- sync_write( cf->ests, f_ests, 3 * cf->t, 0 );
- fclose( f_ests );*/
-#endif
- free( data.Y );
- free( data.beta );
-}
-
-double polygenic_REML_logLik_chol (
- int n, int widthXL,
- double *Phi, double *X, double *Y,
- double sigma2, double h2,
- double *loglik, double *res_sigma, double *beta
-)
-{
- double //*beta, //
- *initX, *initY, // Copy of the initial X and Y
- alpha, gamma, // to scale Phi
- *K, // readibility
- *V, det, // temporary
- ZERO = 0.0, // BLAS / LAPACK
- ONE = 1.0,
- MINUS_ONE = -1.0;
-
- int wXL = widthXL,
- info,
- nn = n*n,
- iONE = 1,
- i; //, j, k;
-
- /*beta = (double *) fgls_malloc ( wXL * sizeof(double) );*/
- initX = (double *) fgls_malloc ( n * wXL * sizeof(double) );
- initY = (double *) fgls_malloc ( n * sizeof(double) );
- V = (double *) fgls_malloc ( wXL * wXL * sizeof(double) );
-
- memcpy( initX, X, n * wXL * sizeof(double) );
- memcpy( initY, Y, n * sizeof(double) );
-
- /* 1) K := h^2 Phi - (1-h^2) I */
- K = Phi;
- alpha = h2; // * sigma2;
- gamma = (1 - h2); // * sigma2;
- dscal_(&nn, &alpha, K, &iONE);
- for ( i = 0; i < n; i++ )
- K[i*n + i] = K[i*n + i] + gamma;
-
- /* 2) L * L' = K */
- int *ipiv = (int *) fgls_malloc (n * sizeof(int));
- int lwork = n * 192; // nb
- double *work = (double *) fgls_malloc ( lwork * sizeof(double));
- dsysv_(LOWER, &n, &widthXL, K, &n, ipiv, X, &n, work, &lwork, &info);
- /*dpotrf_(LOWER, &n, K, &n, &info);*/
- if (info != 0)
- {
- fprintf(stderr, __FILE__ ": dsysv_ inv(M)*X - The matrix is singular(%d)\n", info);
- exit(-1);
- }
-
- /* 3) X := inv(L) X */
- /*dtrsm_(LEFT, LOWER, NO_TRANS, NON_UNIT, &n, &wXL, &ONE, K, &n, X, &n);*/
-
- /* 4) y := inv(L) y */
- /*dtrsv_(LOWER, NO_TRANS, NON_UNIT, &n, K, &n, Y, &iONE);*/
- /*dtrsv_(LOWER, NO_TRANS, NON_UNIT, &n, K, &n, &Y[j * n], &iONE);*/
-
- /* 5) beta := X' * y */
- dgemv_(TRANS, &n, &widthXL, &ONE, X, &n, Y, &iONE, &ZERO, beta, &iONE);
-
- /* 6) V := X' * X */
- /*dsyrk_(LOWER, TRANS, &wXL, &n, &ONE, X, &n, &ZERO, V, &wXL);*/
- dgemm_(TRANS, NO_TRANS, &wXL, &wXL, &n, &ONE, initX, &n, X, &n, &ZERO, V, &wXL);
-
- /* Chol(V) */
- /*dpotrf_(LOWER, &wXL, V, &wXL, &info);*/
- dsysv_(LOWER, &wXL, &iONE, V, &wXL, ipiv, beta, &wXL, work, &lwork, &info);
- if (info != 0)
- {
- fprintf(stderr, __FILE__ ": [ERROR] inv(V) y failed - info: %d\n", info);
- exit(-1);
- }
- /* beta := inv( V ) beta */
- /*dtrsv_(LOWER, NO_TRANS, NON_UNIT, &wXL, V, &wXL, beta, &iONE);*/
- /*dtrsv_(LOWER, TRANS, NON_UNIT, &wXL, V, &wXL, beta, &iONE);*/
-
- // y - X beta
- dgemv_(NO_TRANS,
- &n, &wXL,
- &MINUS_ONE, initX, &n, beta, &iONE,
- &ONE, initY, &iONE);
- // initY stores now "y - X beta"
-
- // residual sigma and loglik
- *res_sigma = variance( initY, n );
-
- memcpy( Y, initY, n * sizeof(double) );
- dsytrs_( LOWER, &n, &iONE, K, &n, ipiv, initY, &n, &info );
- if (info != 0)
- {
- fprintf(stderr, __FILE__ ": [ERROR] inv(K) y - info: %d\n", info);
- exit(-1);
- }
- /*dtrsv_(LOWER, NO_TRANS, NON_UNIT, */
- /*&n, K, &n, initY, &iONE);*/
- *loglik = ddot_(&n, Y, &iONE, initY, &iONE) / (*res_sigma);
-
- /*
- * In principle, the det is just product of diagonal entries (D)
- * D can be 2x2 block diagonal
- * Have to check this case and compute accordingly
- */
- det = 0;
- for ( i = 0; i < n; i++ )
- {
- /*det += log(K[i*n + i]);*/
- if ( ipiv[i] > 0 )
- det += log(K[i*n + i]);
- else if ( ipiv[i] < 0 && i > 0 && ipiv[i-1] == ipiv[i] )
- det += log( K[(i-1)*n + (i-1)] * K[i*n + i] -
- K[(i-1)*n + i] * K[(i-1)*n + i] );
- }
- /*det = det * 2;*/
- *loglik = det + (*loglik) + n * log(*res_sigma);
-
- // Clean up
- /*free(beta);*/
- free(initX);
- free(initY);
- free(V);
-
- free( ipiv );
- free( work );
-
- return *loglik;
-}
-
-double polygenic_REML_logLik_chol_wrapper ( double h, void *data )
-{
- estimates_chol_t *chol_t = (estimates_chol_t*)data;
- size_t n = chol_t->n,
- wXL = chol_t->wXL;
- double *Phi, // = chol_t->Phi,
- *Y, // = chol_t->Y,
- *X; // = chol_t->X;
- double sigma = chol_t->sigma;
- double loglik;
-
- Phi = fgls_malloc( n * n * sizeof(double) );
- X = fgls_malloc( n * wXL * sizeof(double) );
- Y = fgls_malloc( n * sizeof(double) );
- memcpy( Phi, chol_t->Phi, n * n * sizeof(double) );
- memcpy( X, chol_t->X, n * wXL * sizeof(double) );
- memcpy( Y, chol_t->Y, n * sizeof(double) );
-
- polygenic_REML_logLik_chol (
- n, wXL,
- Phi, X, Y,
- sigma, h,
- &loglik, &chol_t->res_sigma, chol_t->beta
- );
-
- chol_t->h = h;
-
- free( Phi );
- free( X );
- free( Y );
-
- return loglik;
-}
-
-
/*
* Eigendecomposition-based REML
*/
@@ -287,7 +74,9 @@
j++;
}
// Phi is not used anymore
- free( cf->Phi ); cf->Phi = NULL;
+ // It is now that chol also uses this method to estimate h2 and res_sigma2
+ /*free( cf->Phi ); cf->Phi = NULL;*/
+
// Load XL
data.X = cf->XL;
cf->ZtXL = data.ZtX = fgls_malloc( data.n * data.wXL * sizeof(double) );
@@ -452,11 +241,12 @@
eig_t->h = h;
-/* printf("Loglik: %.15e\n", loglik);
+#if 0
+ printf("Loglik: %.15e\n", loglik);
int i;
for ( i = 0; i < wXL; i++ )
printf("Beta[%d]: %.15e\n", i, eig_t->beta[i]);
-*/
+#endif
return loglik;
}
@@ -485,17 +275,6 @@
exit(-1);
}
-#if 0
- int i;
- int negs = 0;
- for (i = 0; i < n; i++)
- {
- if (W[i] < 0)
- negs++;
- }
- /*printf("# of negative eigenvalues: %d\n", negs);*/
-#endif
-
free( work );
free( iwork );
free( isuppz );
Modified: pkg/OmicABEL/src/REML.h
===================================================================
--- pkg/OmicABEL/src/REML.h 2012-12-18 15:31:36 UTC (rev 1058)
+++ pkg/OmicABEL/src/REML.h 2012-12-19 10:55:08 UTC (rev 1059)
@@ -28,38 +28,6 @@
#include "GWAS.h"
/*
- * REML based on the Cholesky approach
- *
- * Data structure, Cholesky REML, and wrappers
- */
-typedef struct
-{
- // Problem dimensions
- size_t n;
- size_t wXL;
-
- // Data to fit
- double *Phi;
- double *X;
- double *Y;
- double sigma;
-
- // computed estimates
- double h;
- double res_sigma;
- double *beta;
-} estimates_chol_t;
-
-void estimates_chol( FGLS_config_t *cf );
-double polygenic_REML_logLik_chol_wrapper ( double h, void *data );
-double polygenic_REML_logLik_chol(
- int n, int widthXL,
- double *Phi, double *X, double *Y,
- double sigma2, double h2,
- double *loglik, double *res_sigma, double *beta
-);
-
-/*
* REML based on the Eigendecomposition approach
*
* Data structure, Eigen REML, and wrapper
Modified: pkg/OmicABEL/src/fgls_chol.c
===================================================================
--- pkg/OmicABEL/src/fgls_chol.c 2012-12-18 15:31:36 UTC (rev 1058)
+++ pkg/OmicABEL/src/fgls_chol.c 2012-12-19 10:55:08 UTC (rev 1059)
@@ -51,6 +51,11 @@
#include "vt_user.h"
#endif
+/*
+ * Builds Phi as an SPD matrix, after the eigenvalues were fixed
+ * during the REML estimation
+ */
+void build_SPD_Phi( int n, double *eigVecs, double *eigVals, double *Phi );
/*
* Cholesky-based solution of the
@@ -72,7 +77,6 @@
double *M;
double *ests;
double *h2;
- double *sigma2;
double *res_sigma;
double alpha;
double beta;
@@ -106,17 +110,17 @@
if ( cf.y_b != 1 )
{
fprintf(stderr, "\n[Warning] y_b not used (set to 1)\n");
- /*cf.y_b = 1;*/
+ cf.y_b = 1;
}
/* Memory allocation */
// In-core
- Phi = cf.Phi;
+ build_SPD_Phi( cf.n, cf.Z, cf.W, cf.Phi );
+ Phi = cf.Phi;
M = ( double * ) fgls_malloc ( (size_t)cf.n * cf.n * sizeof(double) );
ests = cf.ests;
h2 = ests;
- sigma2 = &ests[cf.t];
res_sigma = &ests[2*cf.t];
XL_orig = cf.XL;
@@ -180,8 +184,8 @@
#endif
/* M := sigma * ( h^2 Phi - (1 - h^2) I ) */
memcpy( M, Phi, (size_t)n * n * sizeof(double) );
- alpha = h2[j]; // * res_sigma2[j];
- beta = (1 - h2[j]); // * res_sigma2[j];
+ alpha = res_sigma[j] * h2[j];
+ beta = res_sigma[j] * (1 - h2[j]);
dscal_(&nn, &alpha, M, &iONE);
for ( i = 0; i < n; i++ )
M[i*n + i] = M[i*n + i] + beta;
@@ -282,7 +286,6 @@
#endif
for (i = 0; i < x_inc; i++)
{
- /*printf("p: %d - one record: %d\n", p, size_one_b_record);*/
id = omp_get_thread_num();
oneB = &tmpBs[ id * p ];
oneV = &tmpVs[ id * p * p ];
@@ -316,10 +319,7 @@
dpotrf_(LOWER, &p, oneV, &p, &info);
if (info != 0)
{
- /*char err[STR_BUFFER_SIZE];*/
- /*snprintf(err, STR_BUFFER_SIZE, "dpotrf(V) failed (info: %d) - i: %d", info, i);*/
- /*error_msg(err, 1);*/
- for ( k = 0; k < (p+ p*(p+1)/2); k++ )
+ for ( k = 0; k < size_one_b_record; k++ )
Bij[k] = 0.0/0.0; //nan("char-sequence");
continue;
}
@@ -334,16 +334,12 @@
snprintf(err, STR_BUFFER_SIZE, "dpotri failed (info: %d)", info);
error_msg(err, 1);
}
- int p2 = p*p;
- dscal_(&p2, &res_sigma[j], oneV, &iONE);
// Copy output
for ( k = 0; k < p; k++ )
Bij[k] = (float) oneB[k];
for ( k = 0; k < p; k++ )
- {
Bij[p+k] = (float)sqrt(oneV[k*p+k]);
- }
int idx = 0;
for ( k = 0; k < p-1; k++ ) // Cols of V
for ( l = k+1; l < p; l++ ) // Rows of V
@@ -352,11 +348,6 @@
idx++;
}
#if 0
- printf("\nh2: %.16e\n", cf.h2[j]);
- printf("s2: %.16e\n" , cf.sigma2[j]);
- printf("res_sigma2: %.16e\n", cf.res_sigma2[j]);
- for ( k = 0; k < (p+ p*(p+1)/2); k++)
- printf("%.16e\n", Bij[k]);
printf("Chi square: %.6f\n", ( (oneB[p-1] / Bij[p+p-1]) * (oneB[p-1] / Bij[p+p-1]) ) );
#endif
}
@@ -369,9 +360,7 @@
#endif
/* Wait until the previous blocks of B's and V's are written */
if ( iter > 0)
- {
double_buffering_wait( &db_B, IO_BUFF );
- }
#if VAMPIR
VT_USER_END("WAIT_BV");
#endif
@@ -419,3 +408,33 @@
return 0;
}
+
+/*
+ * Builds Phi as an SPD matrix, after the eigenvalues were fixed
+ * during the REML estimation
+ *
+ * Z = eigVecs
+ * W = eigVals
+ * Phi = Z W Z^T
+ */
+void build_SPD_Phi( int n, double *eigVecs, double *eigVals, double *Phi )
+{
+ double ONE = 1.0,
+ ZERO = 0.0,
+ *tmp;
+ int i, j;
+
+ tmp = fgls_malloc( n * n * sizeof(double) );
+
+ // tmp = Z * W
+ for ( j = 0; j < n; j++ )
+ for ( i = 0; i < n; i++ )
+ tmp[ j*n + i ] = eigVecs[ j*n + i ] * eigVals[j];
+ // Phi = tmp * Z^T
+ dgemm_( NO_TRANS, TRANS,
+ &n, &n, &n,
+ &ONE, tmp, &n, eigVecs, &n,
+ &ZERO, Phi, &n );
+
+ free( tmp );
+}
Modified: pkg/OmicABEL/src/utils.c
===================================================================
--- pkg/OmicABEL/src/utils.c 2012-12-18 15:31:36 UTC (rev 1058)
+++ pkg/OmicABEL/src/utils.c 2012-12-19 10:55:08 UTC (rev 1059)
@@ -238,7 +238,36 @@
avail_mem = avail_mem - 100 * cf->n; // Extra little vars
y_b = 1, x_b = 1024 * cf->num_threads;
- ooc_b = 0;
+ /*ooc_b = 0;*/
+ // Estimate the block size for the ooc gemms
+ ooc_b = 1024 * nths; // 1024 columns per thread
+ avail_mem = avail_mem / 2; // double buffering for ooc
+ converged = 0;
+ while ( !converged )
+ {
+ /*printf("ooc_b: %d\n", ooc_b);*/
+ if ( ooc_b * cf->n < avail_mem ) // 2 buffers, n columns (Z' Y or Z' XR)
+ {
+ cf->ooc_b = ooc_b; // * cf->num_threads;
+ converged = 1;
+ }
+ else if ( ooc_b == 1 )
+ {
+ fprintf( stderr, "[ERROR] Not enought memory (ooc_b). Please try a system with larger memory.\n" );
+ /*fprintf( stderr, " For testing purposes, please reduce the value of (sample size).\n" );*/
+ exit(EXIT_FAILURE);
+ }
+ else
+ {
+ ooc_b = ooc_b / 2;
+ if ( ooc_b < 512 )
+ {
+ fprintf( stderr, "[WARNING] ooc_b below 512, this might affect performance.\n" );
+ fprintf( stderr, " We suggest to run this test in a system with a larger amount of memory.\n" );
+ /*fprintf( stderr, " For testing purposes, if you observe delays, try reducing n (sample size).\n" );*/
+ }
+ }
+ }
converged = 0;
size_t mem_usage;
@@ -260,7 +289,7 @@
{
cf->x_b = MIN( cf->m, x_b );
cf->y_b = y_b;
- cf->ooc_b = 0;
+ cf->ooc_b = ooc_b;
converged = 1;
}
else if ( x_b == 1 )
More information about the Genabel-commits
mailing list