[GenABEL-dev] [Genabel-commits] r1819 - in pkg/OmicABELnoMM: . examples src tests
L.C. Karssen
lennart at karssen.org
Thu Sep 11 14:35:47 CEST 2014
On 11-09-14 13:45, Frank, Alvaro Jesus wrote:
> The warnings will be removed once a few more functionality is
> present to avoid breaking existing one.
Great! Thanks.
Lennart.
> ________________________________________
> 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
> -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
--
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
L.C. Karssen
Utrecht
The Netherlands
lennart at karssen.org
http://blog.karssen.org
GPG key ID: A88F554A
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 213 bytes
Desc: OpenPGP digital signature
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20140911/d8bb84c6/attachment-0001.sig>
More information about the genabel-devel
mailing list