[GenABEL-dev] [Genabel-commits] r1819 - in pkg/OmicABELnoMM: . examples src tests
Frank, Alvaro Jesus
alvaro.frank at rwth-aachen.de
Thu Sep 11 13:45:13 CEST 2014
The warnings will be removed once a few more functionality is present to avoid breaking existing one.
________________________________________
From: genabel-devel-bounces at lists.r-forge.r-project.org [genabel-devel-bounces at lists.r-forge.r-project.org] on behalf of L.C. Karssen [lennart at karssen.org]
Sent: Wednesday, September 10, 2014 10:39 PM
To: genabel-devel at lists.r-forge.r-project.org
Subject: Re: [GenABEL-dev] [Genabel-commits] r1819 - in pkg/OmicABELnoMM: . examples src tests
Hi Alvaro,
On 09-09-14 15:54, noreply at r-forge.r-project.org wrote:
> Author: afrank
> Date: 2014-09-09 15:54:05 +0200 (Tue, 09 Sep 2014)
> New Revision: 1819
>
> Added:
> pkg/OmicABELnoMM/examples/dosages_2.txt
> Modified:
> pkg/OmicABELnoMM/ChangeLog
> pkg/OmicABELnoMM/configure.ac
> pkg/OmicABELnoMM/src/AIOwrapper.cpp
> pkg/OmicABELnoMM/src/AIOwrapper.h
> pkg/OmicABELnoMM/src/Algorithm.cpp
> pkg/OmicABELnoMM/src/Algorithm.h
> pkg/OmicABELnoMM/src/Definitions.h
> pkg/OmicABELnoMM/src/Utility.cpp
> pkg/OmicABELnoMM/src/main.cpp
> pkg/OmicABELnoMM/tests/test.cpp
> Log:
> Fixed bug related to reusing the same instance of the solver.
> AIOwrapper is now recreated on every call. Added Additive,Recessive,
> Dominant models. Added option for Custom Models. Custom Additive Model
> uses custom factors. Custom Linear Model uses custom models with beta
> coefficients for each column of the independent variable.
Great to see you implemented new genetic models to the code. That's a
great addition and makes OmicABELnoMM more feature-comparable to ProbABEL.
I also noticed a steady increase in the number of cpplint warnings in
Jenkins (see
http://jenkins.genabel.org/jenkins/ob/OmicABELnoMM/39/violations/). A
lot of them seem to have to do with code layout issues like lines that
are longer than 80 characters and missing or too many spaces. It would
be great if you could fix these as it makes the code easier to read (and
thus to maintain). Of course it would be great if you tackle (some of)
the other cpplint issues as well!
Thanks a lot for all the good work,
Lennart.
>
> Modified: pkg/OmicABELnoMM/ChangeLog
> ===================================================================
> --- pkg/OmicABELnoMM/ChangeLog 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/ChangeLog 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -5,20 +5,24 @@
> -Add exclusion lists for single sets of elements of phenotypes
> -Add exclusion lists for single sets of elements of genotypes
> -Compare ID lists of all dvi files to assure correct ordering
> --Allow for runtime dosage models
>
> Optimizations:
>
> -Reduce memcpy overhead of XR and XR XL factors
> --Reduce computation time of XR and XR XL factors (do GEMMS)
>
>
>
> -
> Changes
> -------------
> -------------
>
> +9-9-2014
> +--------------
> +Fixed bug related to reusing the same instance of the solver. AIOwrapper is now recreated on every call.
> +Added Additive,Recessive, Dominant models.
> +Added option for Custom Models. Custom Additive Model uses custom factors.
> +Custom Linear Model uses custom models with beta coefficients for each column of the independent variable.
> +
> 8-9-2014
> --------------
> Removed individuals with covariates missing
>
> Modified: pkg/OmicABELnoMM/configure.ac
> ===================================================================
> --- pkg/OmicABELnoMM/configure.ac 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/configure.ac 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -18,8 +18,8 @@
> # Set some default compile flags
> if test -z "$CXXFLAGS"; then
> # User did not set CXXFLAGS, so we can put in our own defaults
> - CXXFLAGS=" -O3 -march=corei7 -mfpmath=sse -mtune=corei7 -flto -funroll-loops"
> - #CXXFLAGS="-g -ggdb"
> + #CXXFLAGS=" -O3 -march=corei7 -mfpmath=sse -mtune=corei7 -flto -funroll-loops"
> + CXXFLAGS="-g -ggdb"
> fi
> if test -z "$CPPFLAGS"; then
> # User did not set CPPFLAGS, so we can put in our own defaults
> @@ -37,7 +37,7 @@
> AC_OPENMP
> AC_SUBST(AM_CXXFLAGS, "$OPENMP_CFLAGS")
>
> -AM_CXXFLAGS="-static -O3 -march=corei7 -mfpmath=sse -mtune=corei7 -flto -funroll-loops -I../libs/include -I./libs/include $AM_CXXFLAGS"
> +AM_CXXFLAGS="-static -g -ggdb -I../libs/include -I./libs/include $AM_CXXFLAGS"
> #AM_CXXFLAGS="-static -I../libs/include -I./libs/include $AM_CXXFLAGS"
> # Checks for libraries.
> # pthread library
>
> Added: pkg/OmicABELnoMM/examples/dosages_2.txt
> ===================================================================
> --- pkg/OmicABELnoMM/examples/dosages_2.txt (rev 0)
> +++ pkg/OmicABELnoMM/examples/dosages_2.txt 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -0,0 +1 @@
> +2 1
> \ No newline at end of file
>
> Modified: pkg/OmicABELnoMM/src/AIOwrapper.cpp
> ===================================================================
> --- pkg/OmicABELnoMM/src/AIOwrapper.cpp 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/src/AIOwrapper.cpp 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -31,8 +31,17 @@
> Fhandler->fakefiles = params.use_fake_files;
>
>
> + Fhandler->use_dosages = params.dosages;
> + if(params.dosages && Fhandler->model ==-1)
> + {
> + cout << "Requested dosages model wihtout a valid model!" << endl;
> + exit(1);
> + }
> + Fhandler->not_done = true;
> + Fhandler->model = params.model;
> + Fhandler->fname_dosages = params.fname_dosages;
> +
>
> - Fhandler->not_done = true;
>
> if(!Fhandler->fakefiles)
> {
> @@ -47,8 +56,9 @@
> Fhandler->storePInd = params.storePInd;
>
> Fhandler->min_p_disp = params.minPdisp;
> - Fhandler->min_R2_disp = params.minR2disp;
> + Fhandler->min_R2_disp = params.minR2disp;
>
> +
> Yfvi = load_databel_fvi( (Fhandler->fnameY+".fvi").c_str() );
> ALfvi = load_databel_fvi( (Fhandler->fnameAL+".fvi").c_str() );
> ARfvi = load_databel_fvi( (Fhandler->fnameAR+".fvi").c_str() );
> @@ -56,7 +66,8 @@
>
>
> params.n = ALfvi->fvi_header.numObservations;
> - Fhandler->fileN = params.n;
> + Fhandler->fileN = params.n;
> + Fhandler->fileR = params.r;
> params.m = ARfvi->fvi_header.numVariables/params.r;
> params.t = Yfvi->fvi_header.numVariables;
> params.l = ALfvi->fvi_header.numVariables;
> @@ -81,12 +92,24 @@
>
>
> int Aname_idx=params.n*ARfvi->fvi_header.namelength;//skip the names of the rows
> - for(int i = 0; i < params.m*params.r; i++)
> + if(Fhandler->use_dosages)
> {
> - Fhandler->ARnames.push_back(string(&(ARfvi->fvi_data[Aname_idx])));
> - Aname_idx += ARfvi->fvi_header.namelength;
> + for(int i = 0; i < params.m; i++)
> + {
> + Fhandler->ARnames.push_back(string(&(ARfvi->fvi_data[Aname_idx])));
> + Aname_idx += ARfvi->fvi_header.namelength*Fhandler->fileR;
> + }
> }
> + else
> + {
> + for(int i = 0; i < params.m*params.r; i++)
> + {
> + Fhandler->ARnames.push_back(string(&(ARfvi->fvi_data[Aname_idx])));
> + Aname_idx += ARfvi->fvi_header.namelength;
> + }
> + }
>
> +
> Aname_idx=params.n*ALfvi->fvi_header.namelength;
> for(int i = 0; i < params.l; i++)
> {
> @@ -100,18 +123,17 @@
>
>
> int opt_tb = 1000;
> - int opt_mb = 1000;
> + int opt_mb = 100;
>
> - params.mb = min(params.m, opt_tb);
> - params.tb = min(params.t, opt_mb);
> + params.mb = min(params.m, opt_mb);
> + params.tb = min(params.t, opt_tb);
>
> -
>
>
> }
> else
> {
> -
> + //other params come from outside
> }
>
> //params.fname_excludelist = "exclfile.txt";
> @@ -137,7 +159,60 @@
>
> }
>
> - params.n -= (excl_count + Almissings);
> + params.n -= (excl_count + Almissings);
> +
> + if(params.dosages)
> + {
> +
> + Fhandler->ArDosage = new float[Fhandler->fileR*params.n];
> + Fhandler->dosages = new float[Fhandler->fileR];
> +
> +
> + switch (Fhandler->model)
> + {
> + case -1://nomodel
> +
> + break;
> + case 0://add
> + if(Fhandler->fileR != 3)
> + {
> + cout << "The amount of columns per Independent Variable (--ngpred) is not 3 for a valid Additive Model!" << endl;
> + exit(1);
> + }
> + Fhandler->dosages[0] = 2;Fhandler->dosages[1] = 1;Fhandler->dosages[2] = 0;
> + params.r = 1;
> + Fhandler->add_dosages = true;
> + break;
> + case 1://dom
> + if(Fhandler->fileR != 3)
> + {
> + cout << "The amount of columns per Independent Variable (--ngpred) is not 3 for a valid Dominant Model!" << endl;
> + exit(1);
> + }
> + Fhandler->dosages[0] = 1;Fhandler->dosages[1] = 1;Fhandler->dosages[2] = 0;
> + params.r = 1;
> + Fhandler->add_dosages = true;
> + break;
> + case 2://rec
> + if(Fhandler->fileR != 3)
> + {
> + cout << "The amount of columns per Independent Variable (--ngpred) is not 3 for a valid Recessive Model!" << endl;
> + exit(1);
> + }
> + Fhandler->dosages[0] = 1;Fhandler->dosages[1] = 0;Fhandler->dosages[2] = 0;
> + params.r = 1;
> + Fhandler->add_dosages = true;
> + break;
> + case 3://linear
> + read_dosages(params.fname_dosages,Fhandler->fileR,Fhandler->dosages);
> + break;
> + case 4://additive
> + read_dosages(params.fname_dosages,Fhandler->fileR,Fhandler->dosages);
> + params.r = 1;
> + Fhandler->add_dosages = true;
> + break;
> + }
> + }
>
> params.p = params.l + params.r;
>
> @@ -174,7 +249,18 @@
> fp_InfoResults.write( (char*)&ALfvi->fvi_data[Aname_idx],ALfvi->fvi_header.namelength*(params.l-1)*sizeof(char));
>
> Aname_idx=params.n*ARfvi->fvi_header.namelength;//skip the names of the rows
> - fp_InfoResults.write( (char*)&ARfvi->fvi_data[Aname_idx],ARfvi->fvi_header.namelength*params.r*params.m*sizeof(char));
> + if(Fhandler->use_dosages)
> + {
> + for(int i = 0; i < params.m; i++)
> + {
> + fp_InfoResults.write( (char*)&ARfvi->fvi_data[Aname_idx],ARfvi->fvi_header.namelength*sizeof(char));
> + Aname_idx += Fhandler->fileR*ARfvi->fvi_header.namelength*sizeof(char);
> + }
> + }
> + else
> + {
> + fp_InfoResults.write( (char*)&ARfvi->fvi_data[Aname_idx],Fhandler->fileR*params.m*ARfvi->fvi_header.namelength*sizeof(char));
> + }
>
> int Yname_idx=params.n*Yfvi->fvi_header.namelength;//skip the names of the rows
> fp_InfoResults.write( (char*)&Yfvi->fvi_data[Yname_idx],Yfvi->fvi_header.namelength*params.t*sizeof(char));
> @@ -190,8 +276,8 @@
> // int opt_tb = max(4*2000,opt_block);
> // int opt_mb = max(2000,opt_block);
> //
> -// params.mb = min(params.m,opt_tb);
> -// params.tb = min(params.t,opt_mb);
> + params.mb = min(params.m,params.mb);
> + params.tb = min(params.t,params.tb);
>
> prepare_AL(params.l,params.n);
> prepare_AR( params.mb, params.n, params.m, params.r);
> @@ -231,6 +317,11 @@
> pthread_cond_destroy(&(Fhandler->condition_read));
>
> delete Fhandler->excl_List;
> + if(Fhandler->use_dosages)
> + {
> + delete [](Fhandler->ArDosage);
> + delete [](Fhandler->dosages);
> + }
>
>
>
> @@ -361,7 +452,8 @@
> Fhandler->empty_buffers.pop();
>
>
> - tobeFilled->size = tmp_y_blockSize;
> + tobeFilled->size = tmp_y_blockSize;
> + //cout << "tbz:" << tmp_y_blockSize << " " << flush;
>
> if(Fhandler->fakefiles)
> {
> @@ -454,21 +546,74 @@
> int chunk_size_buff;
> int buff_pos=0;
> int file_pos;
> + float* destination = Fhandler->ArDosage;
>
> - for(int i = 0; i < tmp_ar_blockSize*Fhandler->r; i++)
> + if(Fhandler->use_dosages)
> {
> - for (list< pair<int,int> >::iterator it=excl_List->begin(); it != excl_List->end(); ++it)
> +
> + if(!Fhandler->add_dosages)
> {
> - file_pos = Fhandler->fileN*i + it->first;
> - fseek ( fp_Ar , file_pos*sizeof(type_precision) , SEEK_SET );
> + destination = tobeFilled->buff;//no need to use temp variable
> + }
>
> - chunk_size_buff = it->second;
> - size_t result = fread (&tobeFilled->buff[buff_pos],sizeof(type_precision),chunk_size_buff,fp_Ar); result++;
> - buff_pos += chunk_size_buff;
> + for(int i = 0; i < tmp_ar_blockSize; i++)
> + {
> + buff_pos=0;
> + for(int ii = 0; ii < Fhandler->fileR; ii++)
> + {
> + for (list< pair<int,int> >::iterator it=excl_List->begin(); it != excl_List->end(); ++it)
> + {
> + file_pos = Fhandler->fileN*i + it->first;
> + fseek ( fp_Ar , file_pos*sizeof(type_precision) , SEEK_SET );
>
> + chunk_size_buff = it->second;
>
> + size_t result = fread (&(destination[buff_pos]),sizeof(type_precision),chunk_size_buff,fp_Ar); result++;
> + buff_pos += chunk_size_buff;
> + }
> + }
> +
> + if(Fhandler->add_dosages)
> + {
> + cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
> + Fhandler->n, 1, Fhandler->fileR, 1.0, Fhandler->ArDosage, Fhandler->n, Fhandler->dosages,Fhandler->fileR ,
> + 0.0, &(tobeFilled->buff[i*Fhandler->n]), Fhandler->n);
> + }
> + else
> + {
> + for(int ii = 0; ii < Fhandler->fileR; ii++)
> + {
> + for(int k=0; k < Fhandler->n; k++)
> + {
> + destination[Fhandler->n*ii+k] *= Fhandler->dosages[ii];
> + }
> + }
> + }
> +
> }
> +
> +
> +
> }
> + else
> + {
> + for(int i = 0; i < tmp_ar_blockSize*Fhandler->r; i++)
> + {
> + for (list< pair<int,int> >::iterator it=excl_List->begin(); it != excl_List->end(); ++it)
> + {
> + file_pos = Fhandler->fileN*i + it->first;
> + fseek ( fp_Ar , file_pos*sizeof(type_precision) , SEEK_SET );
> +
> + chunk_size_buff = it->second;
> + size_t result = fread (&tobeFilled->buff[buff_pos],sizeof(type_precision),chunk_size_buff,fp_Ar); result++;
> + buff_pos += chunk_size_buff;
> +
> +
> + }
> + }
> + }
> +
> +
>
>
> }
> @@ -702,6 +847,7 @@
> Fhandler->write_empty_buffers.pop();
> delete tmp2;
> }
> +
> }
>
>
> @@ -1016,7 +1162,8 @@
> void AIOwrapper::prepare_AR( int desired_blockSize, int n, int totalR, int columnsAR)
> {
>
> - Fhandler->Ar = new type_precision[desired_blockSize*columnsAR*n];
> + Fhandler->Ar = new type_precision[desired_blockSize*columnsAR*n];
> +
> Fhandler->Ar_blockSize = desired_blockSize;
> Fhandler->r = columnsAR;
> Fhandler->Ar_Amount = totalR;
> @@ -1352,6 +1499,29 @@
>
> }
>
> +
> +void AIOwrapper::read_dosages(string fname_dosages, int expected_count, float* vec)
> +{
> + ifstream fp_dos(fname_dosages.c_str());
> + if(fp_dos == 0)
> + {
> + cout << "Error reading dosages file."<< endl;
> + exit(1);
> + }
> + int i;
> + for (i=0; i < expected_count && !fp_dos.eof(); i++)
> + {
> + fp_dos >> vec[i];
> + //cout << vec[i];
> + }
> + if(i!=expected_count)
> + {
> + cout << "not enough factor for the dosage model! required " << expected_count << endl;
> + exit(1);
> + }
> +
> +}
> +
>
> void AIOwrapper::free_databel_fvi( struct databel_fvi **fvi )
> {
>
> Modified: pkg/OmicABELnoMM/src/AIOwrapper.h
> ===================================================================
> --- pkg/OmicABELnoMM/src/AIOwrapper.h 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/src/AIOwrapper.h 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -29,8 +29,10 @@
>
>
> string fnameOutFiles;
> + string fname_dosages;
>
>
> +
> list< pair<int,int> >* excl_List;
>
>
> @@ -48,7 +50,8 @@
> vector< string > ALnames;
>
> type_precision* Yb;
> - type_precision* Ar;
> + type_precision* Ar;
> + type_precision* ArDosage;
> type_precision* AL;
> type_precision* B;
> type_buffElement* currentReadBuff;
> @@ -66,11 +69,14 @@
> queue<type_buffElement*> ar_full_buffers;
>
> int index;
> - int fileN;
> + int fileN;
> + int fileR;
> int n;
> int r;
> int l;
> - int p;
> + int p;
> +
> + int model;
>
> int Ar_Amount;
> int Ar_blockSize;
> @@ -84,10 +90,17 @@
> int max_b_blockSize;
>
> bool not_done;
> - bool reset_wait;
> + bool reset_wait;
> + bool use_dosages;
> + bool add_dosages;
>
> int seed;
> - int Aseed;
> + int Aseed;
> +
> + float* dosages;
> + vector< vector <float> > cov_2_Terms;
> + vector< vector <float> > x_Terms;
> + vector< vector <float> > xcov_2_Terms;
>
> pthread_mutex_t m_more ;
> pthread_cond_t condition_more ;
> @@ -165,7 +178,8 @@
>
> private:
>
> - void read_excludeList(list< pair<int,int> >* excl, int &excl_count, int max_excl, string fname_excludeList);
> + void read_excludeList(list< pair<int,int> >* excl, int &excl_count, int max_excl, string fname_excludeList);
> + void read_dosages(string fname_dosages, int expected_count, float* vec);
>
>
> void prepare_AR( int desired_blockSize, int n, int totalR, int columnsR);
>
> Modified: pkg/OmicABELnoMM/src/Algorithm.cpp
> ===================================================================
> --- pkg/OmicABELnoMM/src/Algorithm.cpp 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/src/Algorithm.cpp 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -396,8 +396,10 @@
> params.disp_cov = false;
> params.storePInd = false;
> params.storeBin = false;
> + params.dosages = false;
> params.threads = 1;
> params.r = 1;
> + params.model = -1;
>
>
> params.minR2store = 0.00001;
> @@ -434,7 +436,7 @@
> if(params.minPdisp > params.minPstore || params.storeBin)
> params.minPstore = params.minPdisp;
>
> -
> + AIOwrapper AIOfile;//leave here to avoid memory errors of reusing old threads
> AIOfile.initialize(params);//THIS HAS TO BE DONE FIRST! ALWAYS
>
> //cout << params.n << "\n";
> @@ -455,7 +457,8 @@
>
>
> int y_amount = params.t;
> - int y_block_size = params.tb; // kk
> + int y_block_size = params.tb; // kk
> + //cout << "yt:"<< y_amount << " oybz:"<<y_block_size << flush;
>
> int a_amount = params.m;
> int a_block_size = params.mb;
> @@ -464,7 +467,7 @@
>
> int y_iters = (y_amount + y_block_size - 1) / y_block_size;
>
> - //cout << y_iters << " " << a_iters << endl;
> + //cout << "yiters:" << y_iters << " aiters:" << a_iters << endl;
>
>
> lda = n;
> @@ -581,11 +584,13 @@
> get_ticks(start_tick2);
>
> AIOfile.load_Yblock(&Y, y_block_size);
> + //cout << "ybz:"<< y_block_size << " " << flush;
>
> get_ticks(end_tick);
> out.acc_loady += ticks2sec(end_tick,start_tick2);
>
> get_ticks(start_tick2);
> +
> replace_nans(&y_nan_idxs[0],y_block_size, Y, n,1);
> sumSquares(Y,y_block_size,n,ssY,y_nan_idxs);
>
>
> Modified: pkg/OmicABELnoMM/src/Algorithm.h
> ===================================================================
> --- pkg/OmicABELnoMM/src/Algorithm.h 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/src/Algorithm.h 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -50,8 +50,8 @@
> protected:
> private:
>
> - AIOwrapper AIOfile;
>
> +
> list < resultH > sigResults;
>
> int max_threads;
>
> Modified: pkg/OmicABELnoMM/src/Definitions.h
> ===================================================================
> --- pkg/OmicABELnoMM/src/Definitions.h 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/src/Definitions.h 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -164,13 +164,16 @@
> float minR2disp;
> float minR2store;
> bool storePInd;
> - bool disp_cov;
> + bool disp_cov;
> + bool dosages;
> + int model;//recessive additive dominant etc
>
> string fnameAL;
> string fnameAR;
> string fnameY;
> string fnameOutFiles;
> - string fname_excludelist;
> + string fname_excludelist;
> + string fname_dosages;
>
> bool doublefileType;
>
>
> Modified: pkg/OmicABELnoMM/src/Utility.cpp
> ===================================================================
> --- pkg/OmicABELnoMM/src/Utility.cpp 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/src/Utility.cpp 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -141,10 +141,12 @@
> int idx = k*cols*rows+i*rows+j;
>
> //cout << idx;
> + if(idx >= rows*cols*vec_blocksize)
> + exit(1);
> +
>
> - if(/*idx < rows*cols*vec_blocksize &&*/ isnan( vec[idx] ))
> - {
> -
> + if(isnan( vec[idx] ))
> + {
> vec[idx] = 0;
> if(indexs_vec)
> {
>
> Modified: pkg/OmicABELnoMM/src/main.cpp
> ===================================================================
> --- pkg/OmicABELnoMM/src/main.cpp 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/src/main.cpp 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -26,7 +26,7 @@
> Optional: \n\t\
> -n --ngpred \t <#SNPcols> Number of columns in the geno file that represent a single SNP \n\t\
> -t --thr \t <#CPUs> Number of computing threads to use to speed computations \n\t\
> --x --excl \t <path/filename> file containing list of individuals to exclude from input files \n\t\
> +-x --excl \t <path/filename> file containing list of individuals to exclude from input files, (see example file) \n\t\
> -d --pdisp \t <0.0~1.0> Value to use as maximum threshold for significance.\n\t\
> \t\t Results with P-values UNDER this threshold will be displayed in the putput .txt file \n\t\
> -r --rdisp \t <-10.0~1.0> Value to use as minimum threshold for R2. \n\t\
> @@ -35,7 +35,17 @@
> -s --psto \t <0.0~1.0> Results with P-values UNDER this threshold will be displayed in the putput binary files \n\t\
> -e --rsto \t <-10.0~1.0> Results with R2-values ABOVE this threshold will be stored in the putput binary files \n\t\
> -i --fdcov \t Flag that forces to include covariates as part of the results that are stored in .txt and binary files \n\t\
> --f --fdgen \t Flag that forces to consider all included results (causes the analisis to ignores ALL threshold values).";
> +-f --fdgen \t Flag that forces to consider all included results (causes the analisis to ignores ALL threshold values). \n\t\
> +-j --additive \t Flag that runs the analisis with an Additive Model with (2*AA,1*AB,0*BB) effects \n\t\
> +-k --dominant \t Flag that runs the analisis with an Dominant Model with (1*AA,1*AB,0*BB) effects \n\t\
> +-l --recessive \t Flag that runs the analisis with an Recessive Model with (1*AA,0*AB,0*BB) effects \n\t\
> +-z --mylinear \t <path/filename> to read Factors 'f_i' for a Custom Linear Model with f1*X1,f2*X2,f3*X3...fn*X_ngpred as effects,\n\t\
> + \t each column of each independent variable will be multiplied with the specified factors. \n\t\
> + \t Formula: y~alpha*cov + beta_1*f1*X1 + beta_2*f2*X2 +...+ beta_n*fn*Xn, (see example files!) \n\t\
> +-y --myaddit \t <path/filename> to read Factors 'f_i' for a Custom Additive Model with (f1*X1,f2*X2,f3*X3...fn*X_ngpred) as effects,\n\t\
> + \t each column of each independent variable will be multiplied with the specified factors and then added together. \n\t\
> + \t Formula: y~alpha*cov + beta*(f1*X1 + f2*X2 +...+ fn*Xn), (see example files!) \n\t\
> +";
>
>
>
> @@ -89,6 +99,11 @@
> {"rsto", required_argument, 0, 'e'},//
> {"fdcov", no_argument, 0, 'i'},//
> {"fdgen", no_argument, 0, 'f'},//
> + {"additive", no_argument, 0, 'j'},//
> + {"dominant", no_argument, 0, 'k'},//
> + {"recessive", no_argument, 0, 'l'},//
> + {"mylinear", required_argument, 0, 'z'},//
> + {"myaddit", required_argument, 0, 'y'},//
> {"help", no_argument, 0, 'h'},//
> {0, 0, 0, 0}
> };
> @@ -96,7 +111,7 @@
> // getopt_long stores the option index here.
> int option_index = 0;
>
> - c = getopt_long(argc, argv, "c:p:g:o:n:t:m:x:d:s:r:e:fibh", long_options, &option_index);
> + c = getopt_long(argc, argv, "c:p:g:o:n:t:m:x:d:s:r:e:z:y:fibhjkl", long_options, &option_index);
>
>
> // Detect the end of the options.
> @@ -220,9 +235,45 @@
> cout << "-f Forcing all included results to be considered independently of max P-val or min R2. (SLOW!)"<< endl;
> break;
>
> + case 'j':
> + params.model = 1;
> + params.dosages = true;
> +
> + cout << "-j Using Additive Model with (2*AA,1*AB,0*BB) effects"<< endl;
> + break;
> +
> + case 'k':
> + params.model = 2;
> + params.dosages = true;
> +
> + cout << "-j Using Dominant Model with (1*AA,1*AB,0*BB) effects"<< endl;
> + break;
> +
> + case 'l':
> + params.model = 3;
> + params.dosages = true;
> +
> + cout << "-j Using Recessive Model with (0*AA,0*AB,1*BB) effects"<< endl;
> + break;
> +
> + case 'z':
> + params.model = 4;
> + params.dosages = true;
> +
> + cout << "-z Using Custom Linear Model with parameters read from the file "<< params.fname_dosages << endl;
> + break;
> +
> + case 'y':
> + params.model = 5;
> + params.dosages = true;
> +
> + cout << "-z Using Custom Additive Model with parameters read from the file "<< params.fname_dosages << endl;
> + break;
> +
> case 'b':
> - params.storeBin = true;
> + params.storeBin = true;
>
> +
> cout << "-b Results will be stored in binary format too"<< endl;
> break;
>
>
> Modified: pkg/OmicABELnoMM/tests/test.cpp
> ===================================================================
> --- pkg/OmicABELnoMM/tests/test.cpp 2014-09-08 14:36:26 UTC (rev 1818)
> +++ pkg/OmicABELnoMM/tests/test.cpp 2014-09-09 13:54:05 UTC (rev 1819)
> @@ -95,9 +95,9 @@
> int factor = 0;
> params.n=2000; params.l=3; params.r=1;
> params.t=800; params.tb=min(800,params.t); params.m=1600; params.mb=min(1600,params.m);
> - alg.solve(params, out2, P_NEQ_B_OPT_MD);
> + //alg.solve(params, out2, P_NEQ_B_OPT_MD);
>
> - print_output(out2, gemm_gflopsPsec);
> + //print_output(out2, gemm_gflopsPsec);
>
>
> cout << "\nDone\n";
> @@ -117,6 +117,9 @@
> params.fnameAR="examples/XR";
> params.fnameY="examples/Y";
> params.fnameOutFiles="resultsSig";
> +// params.dosages = true;
> +// params.model = 4;
> +// params.fname_dosages = "examples/dosages_2.txt";
>
>
> for(int th = 0; th < max_threads; th++)
> @@ -138,6 +141,8 @@
>
> max_threads = 2;
> int iters = 10;
> +
> + //cout << "misc tests" << endl;
>
> for (int th = 1; th < max_threads+1; th++)
> {
>
> _______________________________________________
> Genabel-commits mailing list
> Genabel-commits at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-commits
>
--
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
L.C. Karssen
Utrecht
The Netherlands
lennart at karssen.org
http://blog.karssen.org
GPG key ID: A88F554A
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
More information about the genabel-devel
mailing list