[Genabel-commits] r1045 - in pkg/OmicABEL: . inc src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Dec 11 13:26:41 CET 2012
Author: dfabregat
Date: 2012-12-11 13:26:41 +0100 (Tue, 11 Dec 2012)
New Revision: 1045
Added:
pkg/OmicABEL/Makefile
pkg/OmicABEL/inc/
pkg/OmicABEL/inc/make.inc.gnu-goto
pkg/OmicABEL/inc/make.inc.gnu-mkl
pkg/OmicABEL/inc/make.inc.icc-mkl
pkg/OmicABEL/inc/make.inc.mac
pkg/OmicABEL/inc/make.inc.static
pkg/OmicABEL/make.inc
pkg/OmicABEL/src/CLAK_GWAS.c
pkg/OmicABEL/src/GWAS.c
pkg/OmicABEL/src/GWAS.h
pkg/OmicABEL/src/REML.c
pkg/OmicABEL/src/REML.h
pkg/OmicABEL/src/blas.h
pkg/OmicABEL/src/databel.c
pkg/OmicABEL/src/databel.h
pkg/OmicABEL/src/double_buffering.c
pkg/OmicABEL/src/double_buffering.h
pkg/OmicABEL/src/fgls_chol.c
pkg/OmicABEL/src/fgls_chol.h
pkg/OmicABEL/src/fgls_eigen.c
pkg/OmicABEL/src/fgls_eigen.h
pkg/OmicABEL/src/lapack.h
pkg/OmicABEL/src/ooc_BLAS.c
pkg/OmicABEL/src/ooc_BLAS.h
pkg/OmicABEL/src/optimization.c
pkg/OmicABEL/src/optimization.h
pkg/OmicABEL/src/options.h
pkg/OmicABEL/src/statistics.c
pkg/OmicABEL/src/statistics.h
pkg/OmicABEL/src/timing.c
pkg/OmicABEL/src/timing.h
pkg/OmicABEL/src/utils.c
pkg/OmicABEL/src/utils.h
pkg/OmicABEL/src/wrappers.c
pkg/OmicABEL/src/wrappers.h
Log:
Adding source files and makefiles.
Notice makefiles are not polished. Need to adjust them to your file system.
Added: pkg/OmicABEL/Makefile
===================================================================
--- pkg/OmicABEL/Makefile (rev 0)
+++ pkg/OmicABEL/Makefile 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,50 @@
+include ./make.inc
+
+SRCDIR = ./src
+DRIVER = ./HP-GWAS
+
+CFLAGS+=-g -Wall -I $(SRCDIR)/ # -D__WORDSIZE=64
+LDLIBS += -lm -lrt
+
+SRCS = $(SRCDIR)/CLAK_GWAS.c $(SRCDIR)/fgls_chol.c $(SRCDIR)/fgls_eigen.c $(SRCDIR)/wrappers.c $(SRCDIR)/timing.c $(SRCDIR)/statistics.c $(SRCDIR)/REML.c $(SRCDIR)/optimization.c $(SRCDIR)/ooc_BLAS.c $(SRCDIR)/double_buffering.c $(SRCDIR)/utils.c $(SRCDIR)/GWAS.c $(SRCDIR)/databel.c
+OBJS = $(SRCS:.c=.o)
+
+.PHONY: all clean
+
+all: $(DRIVER)
+
+$(DRIVER): $(OBJS)
+ $(CC) $(CFLAGS) $^ $(LDFLAGS) $(LDLIBS) -o $@
+
+clean:
+ $(RM) $(OBJS)
+ $(RM) $(DRIVER) $(EXE_GPU);
+ $(RM) $(SRCDIR)/*mod*
+ $(RM) $(SRCDIR)/*opari_GPU*
+
+
+src/CLAK_GWAS.o: src/CLAK_GWAS.c src/wrappers.h src/utils.h src/GWAS.h \
+ src/databel.h src/timing.h src/REML.h src/fgls_chol.h src/fgls_eigen.h \
+ src/double_buffering.h
+src/GWAS.o: src/GWAS.c src/utils.h src/GWAS.h src/databel.h src/wrappers.h
+src/REML.o: src/REML.c src/options.h src/blas.h src/lapack.h src/ooc_BLAS.h \
+ src/wrappers.h src/utils.h src/GWAS.h src/databel.h src/statistics.h \
+ src/optimization.h src/REML.h
+src/databel.o: src/databel.c src/databel.h src/wrappers.h
+src/double_buffering.o: src/double_buffering.c src/GWAS.h src/databel.h \
+ src/wrappers.h src/double_buffering.h
+src/fgls_chol.o: src/fgls_chol.c src/blas.h src/lapack.h src/options.h \
+ src/GWAS.h src/databel.h src/wrappers.h src/timing.h \
+ src/double_buffering.h src/utils.h src/fgls_chol.h
+src/fgls_eigen.o: src/fgls_eigen.c src/blas.h src/lapack.h src/options.h \
+ src/GWAS.h src/databel.h src/wrappers.h src/timing.h \
+ src/double_buffering.h src/ooc_BLAS.h src/utils.h src/fgls_eigen.h
+src/ooc_BLAS.o: src/ooc_BLAS.c src/blas.h src/lapack.h src/options.h \
+ src/GWAS.h src/databel.h src/wrappers.h src/utils.h \
+ src/double_buffering.h src/ooc_BLAS.h
+src/optimization.o: src/optimization.c src/optimization.h
+src/statistics.o: src/statistics.c src/statistics.h
+src/timing.o: src/timing.c src/timing.h
+src/utils.o: src/utils.c \
+ src/GWAS.h src/databel.h src/wrappers.h src/utils.h
+src/wrappers.o: src/wrappers.c src/GWAS.h src/databel.h src/wrappers.h
Added: pkg/OmicABEL/inc/make.inc.gnu-goto
===================================================================
--- pkg/OmicABEL/inc/make.inc.gnu-goto (rev 0)
+++ pkg/OmicABEL/inc/make.inc.gnu-goto 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,5 @@
+CC = gcc
+CFLAGS = -D_FILE_OFFSET_BITS=64 -pthread -fopenmp -DGOTO
+
+LDFLAGS =
+LDLIBS = $(HOME)/lapack/lapack_LINUX.a $(HOME)/GotoBLAS2/libgoto2.a -lgfortran -gfortranbegin -lpthread
Added: pkg/OmicABEL/inc/make.inc.gnu-mkl
===================================================================
--- pkg/OmicABEL/inc/make.inc.gnu-mkl (rev 0)
+++ pkg/OmicABEL/inc/make.inc.gnu-mkl 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,5 @@
+CC = gcc
+CFLAGS = -D_FILE_OFFSET_BITS=64 -pthread -fopenmp -DMKL
+
+LDFLAGS = -L$(MKLROOT)/lib/intel64/ -L$(MKLROOT)/lib/em64t/ -L$(INTELROOT)/compiler/lib/intel64
+LDLIBS = -lrwthmkl -liomp5
Added: pkg/OmicABEL/inc/make.inc.icc-mkl
===================================================================
--- pkg/OmicABEL/inc/make.inc.icc-mkl (rev 0)
+++ pkg/OmicABEL/inc/make.inc.icc-mkl 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,5 @@
+CC = icc
+CFLAGS = -D_FILE_OFFSET_BITS=64 -pthread -openmp -DMKL
+
+LDFLAGS = -L$(MKLROOT)/lib/intel64/ -L$(INTELROOT)/compiler/lib/intel64
+LDLIBS = -lrwthmkl
Added: pkg/OmicABEL/inc/make.inc.mac
===================================================================
--- pkg/OmicABEL/inc/make.inc.mac (rev 0)
+++ pkg/OmicABEL/inc/make.inc.mac 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,5 @@
+CC = icc
+CFLAGS = -D_FILE_OFFSET_BITS=64 -pthread -openmp -DMKL
+
+LDFLAGS = -L$(MKLROOT)/lib/ -L/opt/intel/composerxe-2011.1.122/compiler/lib/
+LDLIBS = -lmkl_intel -lmkl_intel_thread -lmkl_core -liomp5 -lmkl_lapack95_ilp64 # -lguide
Added: pkg/OmicABEL/inc/make.inc.static
===================================================================
--- pkg/OmicABEL/inc/make.inc.static (rev 0)
+++ pkg/OmicABEL/inc/make.inc.static 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,5 @@
+CC = gcc
+CFLAGS = -static -static-libgcc -D_FILE_OFFSET_BITS=64 -pthread -fopenmp -DMKL
+
+#LDFLAGS = -L$(MKLROOT)/lib/intel64/ -L$(MKLROOT)/lib/em64t/ -L$(INTELROOT)/compiler/lib/intel64
+LDLIBS = $(MKLROOT)/lib/intel64/librwthmkl.a $(INTELROOT)/compiler/lib/intel64/libiomp5.a
Added: pkg/OmicABEL/make.inc
===================================================================
--- pkg/OmicABEL/make.inc (rev 0)
+++ pkg/OmicABEL/make.inc 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,4 @@
+CC = gcc
+CFLAGS = -static -static-libgcc -D_FILE_OFFSET_BITS=64 -fopenmp -DMKL
+
+LDLIBS = $(MKLROOT)/lib/intel64/librwthmkl.a $(INTELROOT)/compiler/lib/intel64/libiomp5.a
Added: pkg/OmicABEL/src/CLAK_GWAS.c
===================================================================
--- pkg/OmicABEL/src/CLAK_GWAS.c (rev 0)
+++ pkg/OmicABEL/src/CLAK_GWAS.c 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,589 @@
+/*
+ * Copyright (c) 2010-2012, Diego Fabregat-Traver and Paolo Bientinesi.
+ * All rights reserved.
+ *
+ * This file is part of OmicABEL.
+ *
+ * OmicABEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * OmicABEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with OmicABEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ *
+ * Coded by:
+ * Diego Fabregat-Traver (fabregat at aices.rwth-aachen.de)
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <sys/stat.h>
+
+#include "wrappers.h"
+#include "utils.h"
+#include "timing.h"
+#include "REML.h"
+#include "fgls_chol.h"
+#include "fgls_eigen.h"
+
+void usage( void );
+int parse_input( int argc, char *argv[], char *var,
+ char *cov_base, char *phi_base, char *snp_base, char *pheno_base, char *out_base,
+ int *nths, int *thres, int *xtile, int *ytile, int *xb, int *yb, int *write_output );
+void check_input_integrity( FGLS_config_t *cf, char *var,
+ char *cov_base, char *phi_base, char *snp_base, char *pheno_base, char *out_base,
+ int nths, int thres, int xtile, int ytile, int xb, int yb, int write_output );
+void print_info( FGLS_config_t *cf );
+void append_estimates( FGLS_config_t *cf );
+void write_output_info_file( FGLS_config_t *cf );
+void cleanup( FGLS_config_t *cf );
+
+int main( int argc, char *argv[] )
+{
+ char cov_base[STR_BUFFER_SIZE] = "",
+ phi_base[STR_BUFFER_SIZE] = "",
+ snp_base[STR_BUFFER_SIZE] = "",
+ pheno_base[STR_BUFFER_SIZE] = "",
+ out_base[STR_BUFFER_SIZE] = "",
+ var[STR_BUFFER_SIZE] = "chol";
+ int nths = 1, write_output = 1;
+ int xb = -1, yb = -1, xtile = 160, ytile = 160;
+ int thres = 95;
+ int status;
+
+ FGLS_config_t cf;
+ struct timeval start, end;
+
+ status = parse_input( argc, argv, var,
+ cov_base, phi_base, snp_base, pheno_base, out_base,
+ &nths, &thres, &xtile, &ytile, &xb, &yb, &write_output );
+ if ( status == 0 )
+ exit(EXIT_FAILURE);
+
+ check_input_integrity( &cf, var,
+ cov_base, phi_base, snp_base, pheno_base, out_base,
+ nths, thres, xtile, ytile, xb, yb, write_output );
+
+ initialize_config(
+ &cf,
+ cov_base, phi_base, snp_base, pheno_base, out_base //,
+ // 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);
+ }
+
+ printf("\n**************************************************************************\n");
+ printf("*** This is an alpha version under development. Use it at your own risk!!!\n");
+ printf("**************************************************************************\n\n");
+ print_info( &cf );
+
+ // Compute h and sigma
+ //
+ // Perform GWAS
+ if ( !strcmp( var, "eigen" ) )
+ {
+ printf( "Estimating GWAS parameters: heritability and variance... " );
+ fflush( stdout );
+ gettimeofday(&start, NULL);
+ estimates_eig( &cf ); // Computation
+ gettimeofday(&end, NULL);
+ printf( "Done ( took %.3f secs )\n", get_diff_ms(&start, &end)/1000. );
+ fflush( stdout );
+
+ printf( "Performing the study... " );
+ fflush( stdout );
+ gettimeofday(&start, NULL);
+ fgls_eigen( &cf ); // Computation
+ gettimeofday(&end, NULL);
+ printf( "Done ( took %.3f secs )\n", get_diff_ms(&start, &end)/1000. );
+ fflush( stdout );
+ }
+ else if ( !strcmp( var, "chol" ) )
+ {
+ printf( "Estimating GWAS parameters: heritability and variance... " );
+ fflush( stdout );
+ gettimeofday(&start, NULL);
+ estimates_chol( &cf ); // Computation
+ gettimeofday(&end, NULL);
+ printf( "Done ( took %.3f secs )\n", get_diff_ms(&start, &end)/1000. );
+ fflush( stdout );
+
+ printf( "Performing the study... " );
+ fflush( stdout );
+ gettimeofday(&start, NULL);
+ fgls_chol( cf ); // Computation
+ gettimeofday(&end, NULL);
+ printf( "Done ( took %.3f secs )\n", get_diff_ms(&start, &end)/1000. );
+ fflush( stdout );
+ }
+
+ append_estimates( &cf );
+ write_output_info_file( &cf );
+ if ( !strcmp( var, "eigen" ) )
+ cleanup( &cf );
+
+ finalize_config( &cf );
+
+ return 0;
+}
+
+void usage( void )
+{
+ fprintf(stderr, "\nUsage: HP-GWAS [options]\n\n");
+ fprintf(stderr, "Following options may be given:\n\n");
+ fprintf(stderr, " -var [chol | eigen] Default is chol.\n");
+ fprintf(stderr, " -cov base path to covariates file \n");
+ fprintf(stderr, " -phi base path to kinship? file \n");
+ fprintf(stderr, " -snp base path to SNPs file \n");
+ fprintf(stderr, " -pheno base path to phenotypes file \n");
+ fprintf(stderr, " -out base path to output file \n");
+ fprintf(stderr, " -nths # Default is 1 thread.\n");
+ fprintf(stderr, " -thres # Default is 95%%.\n");
+ /*fprintf(stderr, " -no-output Default is to write the output files.\n");*/
+ fprintf(stderr, " -h Show this help and exit\n\n");
+}
+
+int parse_input( int argc, char *argv[], char *var,
+ char *cov_base, char *phi_base, char *snp_base, char *pheno_base, char *out_base,
+ int *nths, int *thres, int *xtile, int *ytile, int *xb, int *yb, int *write_output )
+{
+ int arg = 1;
+ int ok = 1;
+ char *err = NULL;
+
+ while (arg < argc && ok)
+ {
+ if ( !strcmp(argv[arg], "-h") )
+ {
+ usage();
+ ok = 0;
+ }
+ else if ( ! strcmp(argv[arg], "-var") )
+ {
+ if ( ++arg < argc )
+ {
+ if ( !strcmp( argv[arg], "chol" ) )
+ {
+ strncpy( var, argv[arg], STR_BUFFER_SIZE );
+ /*fprintf( stderr, "Chol variant is currently broken\n" );*/
+ /*exit( EXIT_FAILURE );*/
+ }
+ else if ( !strcmp( argv[arg], "eigen" ) )
+ strncpy( var, argv[arg], STR_BUFFER_SIZE );
+ else
+ {
+ fprintf( stderr, "Unrecognized variant: %s\n\n", argv[arg] );
+ ok = 0;
+ }
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-cov") )
+ {
+ if ( ++arg < argc )
+ {
+ /*dir = argv[arg];*/
+ strncpy( cov_base, argv[arg], STR_BUFFER_SIZE );
+/* if ( ! ( stat( dir, &st ) == 0 && S_ISDIR( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a directory\n\n", argv[arg] );
+ ok = 0;
+ }*/
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-phi") )
+ {
+ if ( ++arg < argc )
+ {
+ strncpy( phi_base, argv[arg], STR_BUFFER_SIZE );
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-snp") )
+ {
+ if ( ++arg < argc )
+ {
+ strncpy( snp_base, argv[arg], STR_BUFFER_SIZE );
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-pheno") )
+ {
+ if ( ++arg < argc )
+ {
+ strncpy( pheno_base, argv[arg], STR_BUFFER_SIZE );
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-out") )
+ {
+ if ( ++arg < argc )
+ {
+ strncpy( out_base, argv[arg], STR_BUFFER_SIZE );
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-nths") )
+ {
+ if ( ++arg < argc )
+ {
+ *nths = atoi(argv[arg]);
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-thres") )
+ {
+ if ( ++arg < argc )
+ {
+ *thres = atoi(argv[arg]);
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-xtile") )
+ {
+ if ( ++arg < argc )
+ {
+ *xtile = atoi(argv[arg]);
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-ytile") )
+ {
+ if ( ++arg < argc )
+ {
+ *ytile = atoi(argv[arg]);
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-xb") )
+ {
+ if ( ++arg < argc )
+ {
+ *xb = atoi(argv[arg]);
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-yb") )
+ {
+ if ( ++arg < argc )
+ {
+ *yb = atoi(argv[arg]);
+ arg++;
+ }
+ else
+ {
+ ok = 0;
+ err = argv[arg-1];
+ }
+ }
+ else if ( ! strcmp(argv[arg], "-no-output") )
+ {
+ *write_output = 0;
+ arg++;
+ }
+ else
+ {
+ fprintf(stderr, "Unrecognized input: %s\n\n", argv[arg]);
+ ok = 0;
+ arg++;
+ }
+
+ if ( err )
+ {
+ fprintf(stderr, "Error: Missing value for argument '%s'\n\n", err);
+ err = NULL;
+ }
+ }
+ return ok;
+}
+
+void check_input_integrity( FGLS_config_t *cf, char *var,
+ char *cov_base, char *phi_base, char *snp_base, char *pheno_base, char *out_base,
+ int nths, int thres, int xtile, int ytile, int xb, int yb, int write_output )
+{
+ struct stat st;
+
+ // Variant: chol or eigen
+ strncpy( cf->var, var, STR_BUFFER_SIZE );
+
+ // Paths provided?
+ if ( strcmp( phi_base, "" ) == 0 || strcmp( cov_base, "" ) == 0 ||
+ strcmp( snp_base, "" ) == 0 || strcmp( pheno_base, "" ) == 0 )
+ {
+ fprintf( stderr, "You must provide all of -cov -phi -snp -pheno\n\n" );
+ exit( EXIT_FAILURE );
+ }
+ if ( strcmp( out_base, "" ) == 0 && write_output )
+ {
+ fprintf( stderr, "Output base path required\n\n" );
+ exit( EXIT_FAILURE );
+ }
+
+ // Phi data
+ snprintf( cf->Phi_data_path, STR_BUFFER_SIZE, "%s.fvd", phi_base );
+ if ( ! ( stat( cf->Phi_data_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->Phi_data_path );
+ exit( EXIT_FAILURE );
+ }
+ // Phi info
+ snprintf( cf->Phi_info_path, STR_BUFFER_SIZE, "%s.fvi", phi_base );
+ if ( ! ( stat( cf->Phi_info_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->Phi_info_path );
+ exit( EXIT_FAILURE );
+ }
+ // Covariates (XL) data
+ snprintf( cf->XL_data_path, STR_BUFFER_SIZE, "%s.fvd", cov_base );
+ if ( ! ( stat( cf->XL_data_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->XL_data_path );
+ exit( EXIT_FAILURE );
+ }
+ // Covariates (XL) info
+ snprintf( cf->XL_info_path, STR_BUFFER_SIZE, "%s.fvi", cov_base );
+ if ( ! ( stat( cf->XL_info_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->XL_info_path );
+ exit( EXIT_FAILURE );
+ }
+ // SNPs (XR) data
+ snprintf( cf->XR_data_path, STR_BUFFER_SIZE, "%s.fvd", snp_base );
+ if ( ! ( stat( cf->XR_data_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->XR_data_path );
+ exit( EXIT_FAILURE );
+ }
+ // SNPs (XR) info
+ snprintf( cf->XR_info_path, STR_BUFFER_SIZE, "%s.fvi", snp_base );
+ if ( ! ( stat( cf->XR_info_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->XR_info_path );
+ exit( EXIT_FAILURE );
+ }
+ // Phenotypes (Y) data
+ snprintf( cf->Y_data_path, STR_BUFFER_SIZE, "%s.fvd", pheno_base );
+ if ( ! ( stat( cf->Y_data_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->Y_data_path );
+ exit( EXIT_FAILURE );
+ }
+ // SNPs (XR) info
+ snprintf( cf->Y_info_path, STR_BUFFER_SIZE, "%s.fvi", pheno_base );
+ if ( ! ( stat( cf->Y_info_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->Y_info_path );
+ exit( EXIT_FAILURE );
+ }
+ // Output (B) data
+ if ( write_output )
+ {
+ snprintf( cf->B_data_path, STR_BUFFER_SIZE, "%s.out", out_base );
+ snprintf( cf->B_info_path, STR_BUFFER_SIZE, "%s.iout", out_base );
+ // What to do if exists?
+/* if ( ! ( stat( cf->Y_data_path, &st ) == 0 && S_ISREG( st.st_mode ) ) )
+ {
+ fprintf( stderr, "%s does not exist or it is not a regular file\n\n", cf->Y_data_path );
+ exit( EXIT_FAILURE );
+ }
+*/ }
+ else
+ {
+ sprintf( cf->B_data_path, "/dev/null" );
+ sprintf( cf->B_info_path, "/dev/null" );
+ }
+
+ // Threads
+ cf->num_threads = nths;
+ if ( cf->num_threads < 1 )
+ {
+ fprintf( stderr, "Number of threads (%d) must be >= 1\n\n", cf->num_threads );
+ exit( EXIT_FAILURE );
+ }
+
+ // Threshold
+ cf->threshold = thres;
+ if ( cf->threshold < 1 || cf->threshold > 100 )
+ {
+ fprintf( stderr, "Threshold (%d) must be an integer value between 1 and 100\n\n", cf->threshold );
+ exit( EXIT_FAILURE );
+ }
+
+ // Tiling and Blocking (hidden options for performance testing)
+ cf->x_b = xb;
+ cf->y_b = yb;
+ cf->x_tile = xtile;
+ cf->y_tile = ytile;
+}
+
+void print_info( FGLS_config_t *cf )
+{
+ printf( "\nRunning a Genome-Wide Association Study of the following size:\n" );
+ printf( " sample size: %zu\n", cf->n);
+ printf( " # of covariates: %zu\n", cf->p-2);
+ printf( " # of SNPs: %zu\n", cf->m);
+ printf( " # of phenotypes: %zu\n", cf->t);
+ fflush( stdout );
+
+ if ( cf->var[0] == 'c' )
+ printf( "\nWill use OmicABEL-Chol with the following parameters\n" );
+ else
+ printf( "\nWill use OmicABEL-Eigen with the following parameters:\n" );
+ printf( " x_b: %zu\n", cf->x_b );
+ printf( " y_b: %zu\n", cf->y_b );
+ printf( " o_b: %zu\n", cf->ooc_b );
+ printf( " x_tile: %zu\n", cf->x_tile );
+ printf( " y_tile: %zu\n", cf->y_tile );
+ printf( " # of threads: %d\n", cf->num_threads );
+ printf( " Available memory %f GBs (out of %f GBs)\n\n", ((double)cf->availMem) / (1L<<30), ((double)cf->totalMem) / (1L<<30) );
+ fflush( stdout );
+}
+
+void write_output_info_file( FGLS_config_t *cf )
+{
+ FILE *f;
+ int type = 6; //DOUBLE;
+ int nbytes = sizeof(double); // bytes per record
+ /*int size_one_output_record = cf->p + cf->p*(cf->p+1)/2;*/
+ int namelength = 32;
+ char buff[STR_BUFFER_SIZE];
+ int i, j;
+
+ f = fgls_fopen( cf->B_info_path, "wb" );
+
+ fwrite( &type, sizeof(int), 1, f );
+ fwrite( &nbytes, sizeof(int), 1, f );
+ fwrite( &cf->p, sizeof(int), 1, f );
+ /*fwrite( &size_one_output_record, sizeof(int), 1, f );*/
+ fwrite( &cf->m, sizeof(int), 1, f );
+ fwrite( &cf->t, sizeof(int), 1, f );
+ fwrite( &cf->x_b, sizeof(int), 1, f );
+ fwrite( &cf->y_b, sizeof(int), 1, f );
+ fwrite( &namelength, sizeof(int), 1, f );
+ // Output labels
+ // beta
+ for ( i = 0; i < cf->wXL; i++ )
+ {
+ size_t cov_pos = (cf->XL_fvi->fvi_header.numObservations + i)*namelength;
+ snprintf( buff, STR_BUFFER_SIZE, "beta_%s", &cf->XL_fvi->fvi_data[cov_pos] );
+ fwrite( buff, namelength, 1, f );
+ }
+ snprintf( buff, STR_BUFFER_SIZE, "beta_SNP" );
+ fwrite( buff, namelength, 1, f );
+ // se
+ for ( i = 0; i < cf->wXL; i++ )
+ {
+ size_t cov_pos = (cf->XL_fvi->fvi_header.numObservations + i)*namelength;
+ snprintf( buff, STR_BUFFER_SIZE, "se_%s", &cf->XL_fvi->fvi_data[cov_pos] );
+ fwrite( buff, namelength, 1, f );
+ }
+ snprintf( buff, STR_BUFFER_SIZE, "se_SNP" );
+ fwrite( buff, namelength, 1, f );
+
+ // cov
+ for ( j = 0; j < cf->wXL; j++ )
+ {
+ size_t cov_pos_j = (cf->XL_fvi->fvi_header.numObservations + j)*namelength;
+ for ( i = j+1; i < cf->wXL; i++ )
+ {
+ size_t cov_pos_i = (cf->XL_fvi->fvi_header.numObservations + i)*namelength;
+ snprintf( buff, STR_BUFFER_SIZE, "cov_%s_%s",
+ &cf->XL_fvi->fvi_data[cov_pos_i],
+ &cf->XL_fvi->fvi_data[cov_pos_j] );
+ fwrite( buff, namelength, 1, f );
+ }
+ snprintf( buff, STR_BUFFER_SIZE, "cov_SNP_%s",
+ &cf->XL_fvi->fvi_data[cov_pos_j] );
+ fwrite( buff, namelength, 1, f );
+ }
+
+ size_t vars_pos = cf->XR_fvi->fvi_header.numObservations*namelength;
+ fwrite( &cf->XR_fvi->fvi_data[vars_pos], namelength, cf->m, f );
+ vars_pos = cf->Y_fvi->fvi_header.numObservations*namelength;
+ fwrite( &cf->Y_fvi->fvi_data[vars_pos], namelength, cf->t, f );
+
+ fclose( f );
+}
+
+void append_estimates( FGLS_config_t *cf )
+{
+ int size_of_one_record = cf->p + cf->p*(cf->p+1)/2;
+ sync_write( cf->ests, cf->B, (3+cf->wXL)*cf->t, size_of_one_record * cf->m * cf->t );
+}
+
+void cleanup( FGLS_config_t *cf )
+{
+ if ( remove( cf->ZtXR_path ) == -1 )
+ perror("[ERROR] Couldn't delete temporary file for XR");
+ if ( remove( cf->ZtY_path ) == -1 )
+ perror("[ERROR] Couldn't delete temporary file for Y");
+}
Added: pkg/OmicABEL/src/GWAS.c
===================================================================
--- pkg/OmicABEL/src/GWAS.c (rev 0)
+++ pkg/OmicABEL/src/GWAS.c 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,198 @@
+/*
+ * Copyright (c) 2010-2012, Diego Fabregat-Traver and Paolo Bientinesi.
+ * All rights reserved.
+ *
+ * This file is part of OmicABEL.
+ *
+ * OmicABEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * OmicABEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with OmicABEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ *
+ * Coded by:
+ * Diego Fabregat-Traver (fabregat at aices.rwth-aachen.de)
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <sys/time.h>
+
+#include <assert.h>
+
+#include "utils.h"
+#include "wrappers.h"
+#include "GWAS.h"
+
+static void free_and_nullify( double **p );
+static void close_and_nullify( FILE **f );
+static void load_databel_info( FGLS_config_t *cf );
+
+//
+// GWAS config + data
+//
+void initialize_config(
+ FGLS_config_t *cf,
+ char *cov_base, char *phi_base, char *snp_base, char *pheno_base, char *out_base//,
+// char *var, int num_threads, int xtile, int ytile, int xb, int yb, int write_output
+)
+{
+ load_databel_info( cf );
+ // Problem dimensions
+ cf->n = cf->XL_fvi->fvi_header.numObservations;
+ cf->p = cf->XL_fvi->fvi_header.numVariables + 1; // Intercept included
+ cf->m = cf->XR_fvi->fvi_header.numVariables;
+ cf->t = cf->Y_fvi->fvi_header.numVariables;
+ // Assuming wXR = 1
+ cf->wXL = cf->p - 1;
+ cf->wXR = 1;
+ // Algorithm parameters
+ /*cf->x_b = x_b;*/
+ /*cf->y_b = y_b;*/
+ /*cf->num_threads = num_threads;*/
+ get_main_memory_size( &cf->totalMem, &cf->availMem );
+ estimate_block_sizes( cf, cf->var[0], cf->x_b == -1 || cf->y_b == -1 ); // if any is -1, estimate them
+/* if ( !(cf->x_b == -1 || cf->y_b == -1) )
+ {
+ cf->x_b = xb;
+ cf->y_b = yb;
+ }
+ cf->x_tile = xtile;
+ cf->y_tile = ytile;
+*/
+ /*cf->x_b = 10;*/
+ /*cf->y_b = 10;*/
+ // In/Out Files
+ /* sprintf( cf->Phi_data_path, "%s.fvd", phi_base );
+ sprintf( cf->Phi_info_path, "%s.fvi", phi_base );
+ sprintf( cf->XL_data_path, "%s.fvd", cov_base );
+ sprintf( cf->XL_info_path, "%s.fvi", cov_base );
+ sprintf( cf->XR_data_path, "%s.fvd", snp_base );
+ sprintf( cf->XR_info_path, "%s.fvi", snp_base );
+ sprintf( cf->Y_data_path, "%s.fvd", pheno_base );
+ sprintf( cf->Y_info_path, "%s.fvi", pheno_base );
+ if ( write_output )
+ {
+ sprintf( cf->B_data_path, "%s.out", out_base );
+ sprintf( cf->B_info_path, "%s.iout", out_base );
+ }
+ else
+ {
+ sprintf( cf->B_data_path, "/dev/null" );
+ sprintf( cf->B_info_path, "/dev/null" );
+ }
+*/
+ // Temporary files
+ snprintf( cf->ZtXL_path, STR_BUFFER_SIZE, "%s.tmp", cov_base );
+ snprintf( cf->ZtXR_path, STR_BUFFER_SIZE, "%s.tmp", snp_base );
+ snprintf( cf->ZtY_path, STR_BUFFER_SIZE, "%s.tmp", pheno_base );
+
+ // Allocate memory and load in-core data
+ // NOTE: only data which is input to GWAS
+ //
+ // Intermediate data as ZtXL will be handled on the fly
+ // by the corresponding routines:
+ // * ZtXL, only available after REML_eigen (if so)
+ // * Z and W will be allocated and filled up during
+ // the first eigen-decomposition of Phi, (if so)
+
+ // ests needed in both cases chol and eigen
+ cf->ests = (double *) fgls_malloc ( (3 + cf->wXL) * cf->t * sizeof(double) );
+ cf->h2 = cf->ests;
+ cf->sigma2 = &cf->ests[ cf->t ];
+ cf->res_sigma2 = &cf->ests[ 2 * cf->t ];
+ cf->beta_ests = &cf->ests[ 3 * cf->t ];
+
+ cf->Phi = (double *) fgls_malloc ( cf->n * cf->n * sizeof(double) );
+ load_file ( cf->Phi_data_path, "rb", cf->Phi, cf->n * cf->n );
+ // Sanity check (Phi)
+ checkNoNans(cf->n * cf->n, cf->Phi, "[ERROR] NaNs not allowed in Phi\n");
+ cf->XL = (double *) fgls_malloc ( cf->n * cf->wXL * sizeof(double) );
+ load_file ( cf->XL_data_path, "rb", cf->XL, cf->n * cf->wXL );
+ // Sanity check (XL)
+ average( cf->XL, cf->n, cf->wXL, cf->threshold, "Covariate",
+ &cf->XL_fvi->fvi_data[cf->n*NAMELENGTH], NAMELENGTH, 1 );
+ cf->ZtXL = NULL;
+ cf->Z = NULL;
+ cf->W = NULL;
+
+ // Open files for out-of-core data
+ // Again, only for data input/output to GWAS. Intermediate, on the fly
+ cf->XR = fgls_fopen( cf->XR_data_path, "rb" );
+ cf->ZtXR = NULL;
+ cf->Y = fgls_fopen( cf->Y_data_path, "rb" );
+ cf->ZtY = NULL;
+ cf->B = fgls_fopen( cf->B_data_path, "wb" );
+
+ return;
+}
+
+void finalize_config( FGLS_config_t *cf )
+{
+ // Free databel headers (fvi)
+ free_databel_fvi( &cf->Phi_fvi );
+ free_databel_fvi( &cf->XL_fvi );
+ free_databel_fvi( &cf->XR_fvi );
+ free_databel_fvi( &cf->Y_fvi );
+ // Free in-core data buffers
+ free_and_nullify( &cf->ests );
+ cf->h2 = NULL;
+ cf->sigma2 = NULL;
+ cf->res_sigma2 = NULL;
+ cf->beta_ests = NULL;
+ free_and_nullify( &cf->Phi );
+ free_and_nullify( &cf->XL );
+ free_and_nullify( &cf->ZtXL );
+ free_and_nullify( &cf->Z );
+ free_and_nullify( &cf->W );
+ // Close out-of-core data files
+ close_and_nullify( &cf->XR );
+ close_and_nullify( &cf->ZtXR );
+ close_and_nullify( &cf->Y );
+ close_and_nullify( &cf->ZtY );
+ close_and_nullify( &cf->B );
+// close_and_nullify( &cf->V );
+}
+
+static void free_and_nullify( double **p )
+{
+ if ( *p )
+ {
+ free( *p );
+ *p = NULL;
+ }
+}
+
+static void close_and_nullify( FILE **f )
+{
+ if ( *f )
+ {
+ fclose( *f );
+ *f = NULL;
+ }
+}
+
+
+void load_databel_info( FGLS_config_t *cf )
+{
+ cf->Phi_fvi = load_databel_fvi( cf->Phi_info_path );
+ cf->XL_fvi = load_databel_fvi( cf->XL_info_path );
+ cf->XR_fvi = load_databel_fvi( cf->XR_info_path );
+ cf->Y_fvi = load_databel_fvi( cf->Y_info_path );
+ // Check integrity, i.e., dimension n matches
+ assert( cf->XL_fvi->fvi_header.numObservations == cf->XR_fvi->fvi_header.numObservations );
+ assert( cf->XL_fvi->fvi_header.numObservations == cf->Y_fvi->fvi_header.numObservations );
+ assert( cf->XL_fvi->fvi_header.numObservations == cf->Phi_fvi->fvi_header.numObservations );
+ assert( cf->XL_fvi->fvi_header.numObservations == cf->Phi_fvi->fvi_header.numVariables );
+ // Assert type is double
+}
Added: pkg/OmicABEL/src/GWAS.h
===================================================================
--- pkg/OmicABEL/src/GWAS.h (rev 0)
+++ pkg/OmicABEL/src/GWAS.h 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,117 @@
+/*
+ * Copyright (c) 2010-2012, Diego Fabregat-Traver and Paolo Bientinesi.
+ * All rights reserved.
+ *
+ * This file is part of OmicABEL.
+ *
+ * OmicABEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * OmicABEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with OmicABEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ *
+ * Coded by:
+ * Diego Fabregat-Traver (fabregat at aices.rwth-aachen.de)
+ */
+
+#ifndef GWAS_H
+#define GWAS_H
+
+#define DEBUG 0
+#define TIMING 0
+#define VAMPIR 0
+#define CHOL_MIX_PARALLELISM 1
+
+#define MIN(x,y) ( (x) < (y) ? (x) : (y) )
+#define STR_BUFFER_SIZE 256
+
+#include <stdio.h>
+#include "databel.h"
+
+
+// GWAS configuration and data carried throughout the computaiton
+typedef struct
+{
+ // Paths to each data files
+ // Zt* are temporary, no databel info + data
+ char Phi_info_path[STR_BUFFER_SIZE];
+ char Phi_data_path[STR_BUFFER_SIZE];
+ char XL_info_path[STR_BUFFER_SIZE];
+ char XL_data_path[STR_BUFFER_SIZE];
+ char XR_info_path[STR_BUFFER_SIZE];
+ char XR_data_path[STR_BUFFER_SIZE];
+ char Y_info_path[STR_BUFFER_SIZE];
+ char Y_data_path[STR_BUFFER_SIZE];
+ char B_info_path[STR_BUFFER_SIZE];
+ char B_data_path[STR_BUFFER_SIZE];
+ // Temporary data
+ char ZtXL_path[STR_BUFFER_SIZE];
+ char ZtXR_path[STR_BUFFER_SIZE];
+ char ZtY_path[STR_BUFFER_SIZE];
+
+ // Headers
+ databel_fvi *Phi_fvi;
+ databel_fvi *XL_fvi;
+ databel_fvi *XR_fvi;
+ databel_fvi *Y_fvi;
+ // Data ready for runtime computations
+ // in-core data -> double *
+ // out-of-core data -> FILE *
+ double *ests;
+ double *h2;
+ double *sigma2;
+ double *res_sigma2; // est_sigma2
+ double *beta_ests;
+ double *Phi;
+ double *XL;
+ double *ZtXL;
+ FILE *XR;
+ FILE *ZtXR;
+ FILE *Y;
+ FILE *ZtY;
+ FILE *B;
+ // Temporary data
+ double *Z; // Eigenvectors of Phi
+ double *W; // Eigenvalues of Phi
+
+ // Problem dimensions
+ size_t n;
+ size_t p;
+ size_t m;
+ size_t t;
+ size_t wXL; // width of XL
+ size_t wXR; // width of XR
+ size_t x_b; // Block-size along X
+ size_t y_b; // Block-size along Y
+ size_t ooc_b; // Block-size for ooc gemms
+
+ // TO-DO
+ // For now, the size of the openmp block
+ // Should be the size of the ooc tile, and xb yb the size of the omp block
+ size_t x_tile;
+ size_t y_tile;
+
+ size_t totalMem;
+ size_t availMem;
+ int num_threads;
+ char var[STR_BUFFER_SIZE];
+ int threshold;
+} FGLS_config_t;
+
+void initialize_config(
+ FGLS_config_t *cf,
+ char *cov_base, char *phi_base, char *snp_base, char *pheno_base, char *out_base//,
+// char *var, int num_threads, int xtile, int ytile, int xb, int yb, int write_output
+);
+
+void finalize_config( FGLS_config_t *cf );
+
+#endif // GWAS_H
Added: pkg/OmicABEL/src/REML.c
===================================================================
--- pkg/OmicABEL/src/REML.c (rev 0)
+++ pkg/OmicABEL/src/REML.c 2012-12-11 12:26:41 UTC (rev 1045)
@@ -0,0 +1,607 @@
+/*
+ * Copyright (c) 2010-2012, Diego Fabregat-Traver and Paolo Bientinesi.
+ * All rights reserved.
+ *
+ * This file is part of OmicABEL.
+ *
+ * OmicABEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * OmicABEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1045
More information about the Genabel-commits
mailing list