[Genabel-commits] r1832 - in pkg/OmicABELnoMM: . examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 8 10:55:04 CEST 2014
Author: afrank
Date: 2014-10-08 10:55:04 +0200 (Wed, 08 Oct 2014)
New Revision: 1832
Modified:
pkg/OmicABELnoMM/ChangeLog
pkg/OmicABELnoMM/examples/CreateData.R
pkg/OmicABELnoMM/examples/XL.fvd
pkg/OmicABELnoMM/examples/XL.fvi
pkg/OmicABELnoMM/examples/XR.fvd
pkg/OmicABELnoMM/examples/XR.fvi
pkg/OmicABELnoMM/examples/Y.fvd
pkg/OmicABELnoMM/examples/Y.fvi
pkg/OmicABELnoMM/examples/exclude_individuals.txt
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/main.cpp
pkg/OmicABELnoMM/tests/test.cpp
Log:
Added naming checks/warning for element ids. Changed calculation of Standard errors to numerically robust versions. Reduced memcpy overhead of XR and XR XL factors. Created efficient calculation with corrections of covariance matrices. Covariants are displayed as results with -i only when the snp/XR is significant. Added a better performing operation for the calculation of residuals. Stable user usable version.
Modified: pkg/OmicABELnoMM/ChangeLog
===================================================================
--- pkg/OmicABELnoMM/ChangeLog 2014-09-23 14:10:02 UTC (rev 1831)
+++ pkg/OmicABELnoMM/ChangeLog 2014-10-08 08:55:04 UTC (rev 1832)
@@ -9,15 +9,21 @@
-Add exclusion lists for single sets of elements of phenotypes
-Add exclusion lists for single sets of elements of genotypes
-Optimizations:
-
--Reduce memcpy overhead of XR and XR XL factors
-
-
Changes
-------------
-------------
+8-10-2014
+--------------
+
+Added naming checks/warning for element ids.
+Changed calculation of Standard errors to numerically robust versions.
+Reduced memcpy overhead of XR and XR XL factors.
+Created efficient calculation with corrections of covariance matrices.
+Covariants are displayed as results with -i only when the snp/XR is significant.
+Added a better performing operation for the calculation of residuals.
+Stable user usable version.
+
22-9-2014
--------------
Fixed countless bugs as a cause of non normal problem dimensions.
Modified: pkg/OmicABELnoMM/examples/CreateData.R
===================================================================
--- pkg/OmicABELnoMM/examples/CreateData.R 2014-09-23 14:10:02 UTC (rev 1831)
+++ pkg/OmicABELnoMM/examples/CreateData.R 2014-10-08 08:55:04 UTC (rev 1832)
@@ -6,8 +6,8 @@
l = 3 # number of covariates+1 for intercept
int = 3
r = 2
-m = r*100 # number of snps
-t = 100 # number of traits
+m = r*1000 # number of snps
+t = 10000 # number of traits
var=0.05
set.seed(1001)
@@ -27,7 +27,7 @@
Y <- matrix(rnorm(t*n),ncol=t)
-for(i in 1:(n*t)){ if(sample(1:100,1) > 85) {Y[i]=0/0} }
+for(i in 1:(n*t)){ if(sample(1:100,1) > 90) {Y[i]=0/0} }
colnames(Y) <- paste("ph",1:t,sep="")
rownames(Y) <- paste("ind",1:n,sep="")
Y_db <- matrix2databel(Y,filename="Y",type="FLOAT")
@@ -61,10 +61,16 @@
}
}
-for(i in 1:(n*m)){ if(sample(1:100,1) > 85) XR[i]=0/0}
+for(i in 1:(n*m)){ if(sample(1:100,1) > 90) XR[i]=0/0}
colnames(XR) <- paste("miss",1:m,sep="")
-for(i in 1:(m/r)){for(j in 1:r) {colnames(XR)[(i-1)*r+(j)] = paste0("snp",paste(i,j,sep="_")) }}
+for(i in 1:(m/r))
+{
+ for(j in 1:r)
+ {
+ colnames(XR)[(i-1)*r+(j)] = paste0("snp",paste(i,j,sep="_") )
+ }
+}
rownames(XR) <- paste("ind",1:n,sep="")
Modified: pkg/OmicABELnoMM/examples/XL.fvd
===================================================================
(Binary files differ)
Modified: pkg/OmicABELnoMM/examples/XL.fvi
===================================================================
(Binary files differ)
Modified: pkg/OmicABELnoMM/examples/XR.fvd
===================================================================
(Binary files differ)
Modified: pkg/OmicABELnoMM/examples/XR.fvi
===================================================================
(Binary files differ)
Modified: pkg/OmicABELnoMM/examples/Y.fvd
===================================================================
(Binary files differ)
Modified: pkg/OmicABELnoMM/examples/Y.fvi
===================================================================
(Binary files differ)
Modified: pkg/OmicABELnoMM/examples/exclude_individuals.txt
===================================================================
--- pkg/OmicABELnoMM/examples/exclude_individuals.txt 2014-09-23 14:10:02 UTC (rev 1831)
+++ pkg/OmicABELnoMM/examples/exclude_individuals.txt 2014-10-08 08:55:04 UTC (rev 1832)
@@ -1 +1 @@
-0 99
+ 1 100
Modified: pkg/OmicABELnoMM/src/AIOwrapper.cpp
===================================================================
--- pkg/OmicABELnoMM/src/AIOwrapper.cpp 2014-09-23 14:10:02 UTC (rev 1831)
+++ pkg/OmicABELnoMM/src/AIOwrapper.cpp 2014-10-08 08:55:04 UTC (rev 1832)
@@ -80,6 +80,7 @@
params.n = min((int)(ALfvi->fvi_header.numObservations),params.limit_n);
Fhandler->fileN = params.n;
+ realN = params.n;
params.m = min((int)(ARfvi->fvi_header.numVariables/params.r),params.limit_m/params.r);
@@ -104,14 +105,37 @@
params.t = min((int)(Yfvi->fvi_header.numVariables),params.limit_t);
params.l = ALfvi->fvi_header.numVariables;
+
+ //!check id names order
+ string YidNames;
+ string ALidName;
+ string ARidName;
+ string INTidName;
+
int yname_idx=0;//starting idx for names on ALfvi->data
for(int i = 0; i < params.n; i++)
{
- //Nnames.push_back(string(&(Yfvi->fvi_data[yname_idx])));
+ YidNames = string(&(Yfvi->fvi_data[yname_idx]));
+ ALidName = string(&(ALfvi->fvi_data[yname_idx]));
+ ARidName = string(&(ARfvi->fvi_data[yname_idx]));
+ if(YidNames.compare(ALidName))
+ {
+ cout << "Warning, ID names between -p file and -c file do not coincide! The files must have the same meaningful order" << endl;
+ i = params.n;//exit
+ }
+ if(YidNames.compare(ARidName))
+ {
+ cout << "Warning, ID names between -p file and -g file do not coincide! The files must have the same meaningful order" << endl;
+ i = params.n;//exit
+ }
+ if(ARidName.compare(ALidName))
+ {
+ cout << "Warning, ID names between -g file and -c file do not coincide! The files must have the same meaningful order" << endl;
+ i = params.n;//exit
+ }
yname_idx += Yfvi->fvi_header.namelength;
- //cout << i << ": " << string(&(Yfvi->fvi_data[yname_idx])) << "\t";
- //cout << i << ": " << Ynames[i] << "\t";
+
}
@@ -135,7 +159,7 @@
int opt_tb = 1000;
- int opt_mb = 100;//4WONT WORK WTF?
+ int opt_mb = 500;
params.mb = min(params.m, opt_mb);
params.tb = min(params.t, opt_tb);
@@ -165,6 +189,7 @@
+
(Fhandler->excl_List)->push_back( make_pair(0,params.n) );
if(params.fname_excludelist.size()!=0)
@@ -287,6 +312,21 @@
exit(1);
}
+ int name_idx= 0;
+ for(int i = 0; i < params.n; i++)
+ {
+ string ARidName = string(&(ARfvi->fvi_data[name_idx]));
+ string INTidName = string(&(Ifvi->fvi_data[name_idx]));
+ if(ARidName.compare(INTidName))
+ {
+ cout << "Warning, ID names between -g file and interactions file do not coincide! The files must have the same meaningful order!" << endl;
+ i = params.n;//exit
+ }
+
+ name_idx += Yfvi->fvi_header.namelength;
+ }
+
+
if(params.end_IDX_interactions != -1 || params.ini_IDX_interactions != -1)
{
Fhandler->numInter = params.end_IDX_interactions - params.ini_IDX_interactions+1;
@@ -733,6 +773,10 @@
int ar_file_pos = 0;
+ int seq_count;
+ int max_secuential_write_count= 10;
+
+
if(!Fhandler->fakefiles)
{
fp_Y = fopen((Fhandler->fnameY+".fvd").c_str(), "rb");
@@ -1200,9 +1244,9 @@
}
//B write
- while(!Fhandler->write_full_buffers.empty() && Local_not_done)
+ while(!Fhandler->write_full_buffers.empty() && Local_not_done && seq_count < max_secuential_write_count)
{
-
+ seq_count++;
//cout << "S" << " ";
pthread_mutex_lock(&(Fhandler->out_buff_upd));
@@ -1258,8 +1302,9 @@
for (list< result >::iterator it2=current.res_p.begin(); it2 != current.res_p.end(); ++it2)
{
result res_p = *it2;
- if(res_p.P <= Fhandler->min_p_disp && current.R2 >= Fhandler->min_R2_disp)
+ if((res_p.P <= Fhandler->min_p_disp && current.R2 >= Fhandler->min_R2_disp )|| Fhandler->storePInd)
{
+ Fhandler->io_overhead = "!";//let the user know a result was found
string Aname = " ";
if(res_p.AL_name_idx < Fhandler->l)
@@ -1323,7 +1368,9 @@
}
+ seq_count = 0;
+
if(!Fhandler->not_done && Fhandler->write_full_buffers.empty())
{Local_not_done = false;}
@@ -1337,7 +1384,8 @@
timeToWait.tv_nsec = time.wMilliseconds*1000 + morenanos ;
#else
clock_gettime(CLOCK_REALTIME, &timeToWait);
- timeToWait.tv_nsec += 10000000;
+ //timeToWait.tv_nsec += 10000000;
+ timeToWait.tv_nsec += 10000;
#endif
pthread_mutex_lock(&(Fhandler->m_more));
@@ -1511,7 +1559,7 @@
pthread_cond_signal( &(Fhandler->condition_more ));
pthread_mutex_unlock(&(Fhandler->m_more));
- io_overhead = "!";
+ io_overhead = "$";
pthread_mutex_lock(&(Fhandler->m_read));
pthread_cond_wait( &(Fhandler->condition_read), &(Fhandler->m_read ));
@@ -1564,7 +1612,7 @@
Fhandler->n= n;
Fhandler->Y_Amount=totalY;
Fhandler->y_to_readSize = Fhandler->Y_Amount;
- Fhandler->buff_count = min(3,(totalY+ y_blockSize - 1)/y_blockSize) ;
+ Fhandler->buff_count = max(3,1+(totalY+ y_blockSize - 1)/y_blockSize) ;
//cout << "buffcount " << Fhandler->buff_count;
@@ -1619,6 +1667,9 @@
pthread_mutex_unlock(&(Fhandler->out_buff_upd));
+ sigResults->clear();
+
+
pthread_mutex_lock(&(Fhandler->m_more));
pthread_cond_signal( &(Fhandler->condition_more ));
pthread_mutex_unlock(&(Fhandler->m_more));
@@ -1634,6 +1685,15 @@
sigResults = 0;
//cout << Fhandler->write_full_buffers.size() << endl;
+ if(Fhandler->io_overhead.compare(""))
+ {
+ io_overhead = Fhandler->io_overhead;
+ }
+ else
+ {
+ Fhandler->io_overhead = "";
+ }
+
pthread_mutex_unlock(&(Fhandler->out_buff_upd));
@@ -1740,7 +1800,7 @@
- int buff_count = min(10,2*(Fhandler->Ar_Amount+desired_blockSize-1)/desired_blockSize+1);
+ int buff_count = max(10,4*(Fhandler->Ar_Amount+desired_blockSize-1)/desired_blockSize+1);
#ifdef DEBUG
cout << "buff_count" << buff_count << endl << flush;
@@ -1769,7 +1829,7 @@
}
-void AIOwrapper::removeALmissings(list< pair<int,int> >* excl_List,struct Settings params, int &Almissings)
+void AIOwrapper::removeALmissings(list< pair<int,int> >* excl_List, struct Settings params, int &Almissings)
{
float* tempAL = new float[params.n*params.l];
@@ -1812,6 +1872,10 @@
}
+
+
+
+
bool AIOwrapper::splitpair(int value, list< pair<int,int> >* excl_List, int n)
{
@@ -2018,8 +2082,15 @@
cout << "\nPlease provide ordered pairs!\n";
exit( 1 );
}
- for(int i = first; i <= second;i++ )
+ if(first < 1 || second < 1)
{
+ cout << "\nPlease provide valid ositive indixes!\n";
+ exit( 1 );
+ }
+
+ for(int i = first-1; i <= second-1;i++ )
+ {
+
bool removed = splitpair(i,excl, n );
if(removed)
{
@@ -2082,7 +2153,7 @@
f = fopen( path, mode );
if ( f == NULL )
{
- cout << "\nerror on fgls_fopen\n";
+ cout << "\nError reading file: " << path << endl;
exit( 1 );
}
return f;
Modified: pkg/OmicABELnoMM/src/AIOwrapper.h
===================================================================
--- pkg/OmicABELnoMM/src/AIOwrapper.h 2014-09-23 14:10:02 UTC (rev 1831)
+++ pkg/OmicABELnoMM/src/AIOwrapper.h 2014-10-08 08:55:04 UTC (rev 1832)
@@ -31,7 +31,9 @@
string fnameAR;
string fnameY;
+ string io_overhead;
+
string fnameOutFiles;
string fname_dosages;
@@ -193,6 +195,7 @@
void write_OutSignificant(list < resultH >* sigResults, int min_p_disp, bool disp_cov);
string io_overhead;
+ int realN;
protected:
Modified: pkg/OmicABELnoMM/src/Algorithm.cpp
===================================================================
--- pkg/OmicABELnoMM/src/Algorithm.cpp 2014-09-23 14:10:02 UTC (rev 1831)
+++ pkg/OmicABELnoMM/src/Algorithm.cpp 2014-10-08 08:55:04 UTC (rev 1832)
@@ -29,7 +29,7 @@
}
-void Algorithm::extract_subMatrix(type_precision* source, type_precision* dest,
+void Algorithm::extract_subMatrix(float* source, float* dest,
int dim1_source, int dim2_source,
int dim1_ini, int dim1_end, int dim2_ini,
int dim2_end)
@@ -41,9 +41,9 @@
j = dim1_ini;
source_ini = i * dim1_source+j;
size = dim1_end - dim1_ini;
- memcpy( (type_precision*)&dest[idx],
- (type_precision*)&source[source_ini],
- size * sizeof(type_precision) );
+ memcpy( (float*)&dest[idx],
+ (float*)&source[source_ini],
+ size * sizeof(float) );
// for (j = dim1_ini; j<dim1_end; j++)
// {
// dest[idx] = source[i*dim1_source+j];
@@ -54,7 +54,7 @@
}
-void Algorithm::prepare_Bfinal(type_precision* bfinal, type_precision* bsource, int a_amount, int y_amount, int p)
+void Algorithm::prepare_Bfinal(float* bfinal, float* bsource, int a_amount, int y_amount, int p)
{
// // memcpy are faster version of the fors
// int i, k, w, top_idx, bot_idx;
@@ -64,9 +64,9 @@
// for (k = 0; k < dim2_b; k++)
// {
// size = k * dim1_b + (dim1_b - dim1_b_bot) - (k * dim1_b);
-// memcpy( (type_precision*)&bfinal[k * dim1_b],
-// (type_precision*)&top[top_idx],
-// size * sizeof(type_precision) );
+// memcpy( (float*)&bfinal[k * dim1_b],
+// (float*)&top[top_idx],
+// size * sizeof(float) );
//// for (i = k*dim1_b; i < k*dim1_b+(dim1_b-dim1_b_bot); i++)
//// {
//// bfinal[i] = top[top_idx];
@@ -77,9 +77,9 @@
// w = i;
//
// size = w + dim1_b_bot - w;
-// memcpy( (type_precision*)&bfinal[w],
-// (type_precision*)&bot[bot_idx],
-// size * sizeof(type_precision) );
+// memcpy( (float*)&bfinal[w],
+// (float*)&bot[bot_idx],
+// size * sizeof(float) );
//// for (i = w; i < w+dim1_b_bot; i++)
//// {
//// bfinal[i] = bot[bot_idx];
@@ -90,8 +90,8 @@
}
-void Algorithm::prepare_QY(type_precision* qy, type_precision* top,
- type_precision* bot, int dim1_QY,
+void Algorithm::prepare_QY(float* qy, float* top,
+ float* bot, int dim1_QY,
int dim2_QY, int dim1_qy_bot, int bot_blocks )
{
int i, k, w, top_idx, bot_idx;
@@ -116,10 +116,10 @@
}
-type_precision* Algorithm::extract_R(type_precision* A, int dim1_A, int dim2_A)
+float* Algorithm::extract_R(float* A, int dim1_A, int dim2_A)
{
- type_precision* R = (type_precision*)calloc(dim2_A * dim2_A,
- sizeof(type_precision));
+ float* R = (float*)calloc(dim2_A * dim2_A,
+ sizeof(float));
int i, j;
int R_idx = 0;
@@ -137,11 +137,11 @@
}
-type_precision* Algorithm::prepare_R(type_precision* RL, int dim1_A,
+float* Algorithm::prepare_R(float* RL, int dim1_A,
int dim2_AL, int dim2_AR)
{
int R_dim = (dim2_AR + dim2_AL);
- type_precision* R = new type_precision[R_dim * R_dim];
+ float* R = new float[R_dim * R_dim];
int i, j;
@@ -161,8 +161,8 @@
}
-void Algorithm::update_R(type_precision* R, type_precision* topRr,
- type_precision* botRr, int dim1, int dim2, int r)
+void Algorithm::update_R(float* R, float* topRr,
+ float* botRr, int dim1, int dim2, int r)
{
int i, j, w;
int max = dim1 * dim2;
@@ -186,8 +186,8 @@
}
-void Algorithm::build_S(type_precision* S, type_precision* Stl,
- type_precision* Str, type_precision* Sbr, int l, int r)
+void Algorithm::build_S(float* S, float* Stl,
+ float* Str, float* Sbr, int l, int r)
{
int Sidx;
int p = l + r;
@@ -224,12 +224,12 @@
}
-void Algorithm::check_result(type_precision* AL, type_precision* AR,
+void Algorithm::check_result(float* AL, float* AR,
int rowsA, int colsA, int rhs, int colsAR,
- type_precision* y, type_precision* res,struct Settings params,int iX,int iiX, int jY, int jjY)
+ float* y, float* res,struct Settings params,int iX,int iiX, int jY, int jjY)
{
- type_precision* A = (type_precision*)malloc(rowsA * colsA *
- sizeof(type_precision));
+ float* A = (float*)malloc(rowsA * colsA *
+ sizeof(float));
int i, ar_idx = 0;
for (i = 0; i < rowsA * (colsA - colsAR); i++)
@@ -246,9 +246,9 @@
//matlab_print_matrix("A", rowsA, colsA, A);
- type_precision* ynew = replicate_vec(y, rowsA*rhs);
- type_precision* new_sol = (type_precision*)malloc(colsA * rhs *
- sizeof(type_precision));
+ float* ynew = replicate_vec(y, rowsA*rhs);
+ float* new_sol = (float*)malloc(colsA * rhs *
+ sizeof(float));
lapack_int info = LAPACKE_sgels(STORAGE_TYPE, 'N', rowsA, colsA, rhs, A,
rowsA, ynew, rowsA);
@@ -266,7 +266,7 @@
index_new += colsA;
}
- type_precision* precomp_betas = new float[colsA];
+ float* precomp_betas = new float[colsA];
if(!params.use_fake_files)
{
@@ -278,36 +278,36 @@
}
else
{
- //type_precision precomp_betas[colsA];
+ //float precomp_betas[colsA];
int file_pos = jY*params.tb*params.m*colsA + jjY*params.m*colsA+iX*params.mb*colsA+iiX*colsA;
- fseek ( fp_Bprecomputed , file_pos*sizeof(type_precision) , SEEK_SET );
+ fseek ( fp_Bprecomputed , file_pos*sizeof(float) , SEEK_SET );
- size_t result = fread (precomp_betas,sizeof(type_precision),colsA,fp_Bprecomputed); result++;
+ size_t result = fread (precomp_betas,sizeof(float),colsA,fp_Bprecomputed); result++;
//matlab_print_matrix("res", colsA, 1, res);
//matlab_print_matrix("new_sol", colsA, 1, new_sol);
//matlab_print_matrix("precomp_betas", colsA, 1, precomp_betas);
cblas_saxpy(colsA, -1.0, res, 1, precomp_betas, 1);
- type_precision u_norm = cblas_snrm2(colsA, precomp_betas, 1);
+ float u_norm = cblas_snrm2(colsA, precomp_betas, 1);
//cout << "pa:"<< u_norm ;
if (fabs(u_norm) > 0.001 || isnan(u_norm))
{
- fseek ( fp_Bprecomputed , file_pos*sizeof(type_precision) , SEEK_SET );
- result = fread (precomp_betas,sizeof(type_precision),colsA,fp_Bprecomputed); result++;
+ fseek ( fp_Bprecomputed , file_pos*sizeof(float) , SEEK_SET );
+ result = fread (precomp_betas,sizeof(float),colsA,fp_Bprecomputed); result++;
cblas_saxpy(colsA, 1.0, res, 1, precomp_betas, 1);
- type_precision u_norm2 = cblas_snrm2(colsA, precomp_betas, 1);
+ float u_norm2 = cblas_snrm2(colsA, precomp_betas, 1);
//cout << u_norm2 << " ";
//cout << file_pos <<"c "<< u_norm2 << "\n";
if(fabs(u_norm2) > 0.001 || isnan(u_norm))
{
- fseek ( fp_Bprecomputed , file_pos*sizeof(type_precision) , SEEK_SET );
- result = fread (precomp_betas,sizeof(type_precision),colsA,fp_Bprecomputed); result++;
+ fseek ( fp_Bprecomputed , file_pos*sizeof(float) , SEEK_SET );
+ result = fread (precomp_betas,sizeof(float),colsA,fp_Bprecomputed); result++;
fflush(stdout);
// matlab_print_matrix("AL", rowsA, colsA-colsAR, AL);
@@ -323,8 +323,8 @@
}
- fseek ( fp_Bprecomputed , file_pos*sizeof(type_precision) , SEEK_SET );
- result = fread (precomp_betas,sizeof(type_precision),colsA,fp_Bprecomputed); result++;
+ fseek ( fp_Bprecomputed , file_pos*sizeof(float) , SEEK_SET );
+ result = fread (precomp_betas,sizeof(float),colsA,fp_Bprecomputed); result++;
cblas_saxpy(colsA, -1.0, new_sol, 1, precomp_betas, 1);
u_norm = cblas_snrm2(colsA, precomp_betas, 1);
//cout << " pl:" << u_norm ;
@@ -346,7 +346,7 @@
// printf("\n Btop=(Rtl\\(Ql'*Y))-(Rtl\\Rtr)*(Rbr\\(Qr'*Y)); \n [Btop ; Rbr\\Qr'*Y] - bcomputed \n");
cblas_saxpy(rhs * colsA, -1.0, res, 1, new_sol, 1);
- type_precision u_norm = cblas_snrm2(rhs * colsA, new_sol, 1);
+ float u_norm = cblas_snrm2(rhs * colsA, new_sol, 1);
//
if (fabs(u_norm) >= 0.1 || isnan(u_norm))
{
@@ -435,13 +435,13 @@
- //type_precision *Ytemp;
+ //float *Ytemp;
lapack_int info, n, lda, l, r, p;
cputime_type start_tick, start_tick2, start_tick3, end_tick;
- if(params.minR2disp > params.minR2store || params.storeBin)
+ if(params.minR2disp < params.minR2store || params.storeBin)
params.minR2store = params.minR2disp;
if(params.minPdisp > params.minPstore || params.storeBin)
@@ -472,11 +472,13 @@
int y_amount = params.t;
- int y_block_size = params.tb; // kk
+ int y_block_size = params.tb;
+ //int orig_y_block_size = y_block_size;
//cout << "yt:"<< y_amount << " oybz:"<<y_block_size << flush;
int a_amount = params.m;
- int a_block_size = params.mb;
+ int a_block_size = params.mb;
+ //int orig_a_block_size = a_block_size;
//cout << r << endl;
@@ -517,54 +519,65 @@
//add memalign
- //type_precision Stl[l*l];
- //type_precision Str[l*r*a_block_size];
- type_precision* Stl = new type_precision[l*l*1];
- type_precision* Str = new type_precision[l*r*a_block_size*1];
+ //float Stl[l*l];
+ //float Str[l*r*a_block_size];
+ float* Stl = new float[l*l*1];
+ float* Str = new float[l*r*a_block_size*1];
- type_precision* Sbr = new type_precision[r * r * a_block_size];
- type_precision* Ay = new type_precision[p * a_block_size];
- //type_precision* B_resorted = new type_precision[p * a_block_size*y_block_size];
- //type_precision* S = new type_precision[p * p];
+ float* Stl_corr = new float[l*l*y_block_size];
+ float* Sbr_corr = new float[r*r*a_block_size*y_block_size];
+ float* Str_corr = new float[l*r*a_block_size*y_block_size];
- type_precision* S2global = new float[p*p*max_threads];
+
+ float* Sbr = new float[r * r * a_block_size];
+ float* Ay = new float[p * a_block_size * y_block_size];
+ //float* B_resorted = new float[p * a_block_size*y_block_size];
+
+ //float* S = new float[p * p];
+
+ float* S2global = new float[p*p*max_threads];
- type_precision* Ay_top = new type_precision[l * y_amount];
- type_precision* Ay_bot = new type_precision[y_block_size * a_block_size * r];
+ float* Ay_top = new float[l * y_amount];
+ float* Ay_bot = new float[y_block_size * a_block_size * r];
- type_precision* y_residual = new type_precision[n * y_block_size ];
- type_precision* y_res_norms = new type_precision[a_block_size];
+ float* y_residual = new float[n * y_block_size ];
+ float* y_res_norms = new float[a_block_size];
//list<long int>* al_nan_idxs = new list<long int>[l];
list<long int>* y_nan_idxs = new list<long int>[y_block_size];
list<long int>* ar_nan_idxs = new list<long int>[a_block_size*r];
int* Ymiss=new int[y_block_size];
+
+ std::map <int , list < int > > y_missings_jj;
- type_precision* A = new type_precision[n * p * 1];
- type_precision* AR = new type_precision[n * r * a_block_size * 1];
+ float* A = new float[n * p * 1];
+ float* AR;// = new float[n * r * a_block_size * 1];
//Sum of squares for A parts
- type_precision* ssAR = new type_precision[r*a_block_size];
- type_precision* ssAL = new type_precision[l];
- type_precision* ssY = new type_precision[y_block_size];
+ float* ssA = new float[p*a_block_size];
- SYY = new type_precision[a_block_size];
+ float* ssY = new float[y_block_size];
+ SYY = new float[a_block_size*y_block_size];
+
+ //float* cholSAL = new float[l*l];
+
+
list < resultH >* sigResults;
-// type_precision* AL = new type_precision[n * l * 1];
- type_precision* AL = A;
+// float* AL = new float[n * l * 1];
+ float* AL = A;
- type_precision* B = Ay;
+ float* B = Ay;
- type_precision* backupAR; // = new type_precision[n*r*a_block_size];
- type_precision* backupAL; // = new type_precision[n*l];
+ float* backupAR; // = new float[n*r*a_block_size];
+ float* backupAL; // = new float[n*l];
AIOfile.load_AL(&backupAL);
@@ -580,16 +593,32 @@
// }
- sumSquares(backupAL,l,n,ssAL,0);
+ //sumSquares(backupAL,l,n,ssAL,0);
+
//LAPACKE_dgesdd()
- copy_vec(backupAL, AL, n*l);
+ copy_vec(backupAL, AL, n*l);
+
+ get_ticks(start_tick2);
+ //! Generate Stl
+ cblas_ssyrk(CblasColMajor, CblasUpper, CblasTrans,
+ l, n, 1.0, AL, lda, 0.0, Stl, l);
+
+ get_ticks(end_tick);
+ out.acc_stl += ticks2sec(end_tick,start_tick2);
+
+
+
+
+
+
- type_precision* Y;
+ float* Y;
+
// printf("\n\n%%Computations\n%%");
@@ -618,17 +647,21 @@
for (int jj = 0; jj < y_block_size; jj++)
{
-// list<long int> nans = al_nan_idxs[0];
-// list<long int> nans2 = y_nan_idxs[jj];
-// nans.merge(nans2);
-// nans.sort();//!might not be needed
-// nans.unique();
+// for (list< long int >::iterator it=y_nan_idxs[jj].begin(); it != y_nan_idxs[jj].end(); ++it)
+// {
+// list< int > jjs = y_missings_jj[*it];
+// jjs.push_back(jj);
+// y_missings_jj[*it] = jjs;
+// }
+
+
Ymiss[jj] = y_nan_idxs[jj].size();
}
//matlab_print_matrix("Y", n, y_block_size, Y);
- //out.acc_other += ticks2sec(end_tick,start_tick2);
+ get_ticks(end_tick);
+ out.acc_other += ticks2sec(end_tick,start_tick2);
get_ticks(start_tick2);
@@ -639,7 +672,14 @@
&Ay_top[j * l * y_block_size], l);
get_ticks(end_tick);
- out.acc_gemm += ticks2sec(end_tick,start_tick2);
+ out.acc_gemm += ticks2sec(end_tick,start_tick2);
+
+
+ get_ticks(start_tick2);
+ generate_correction_missings_Stl(y_block_size,1,1,n,l,l,AL,AL,Stl_corr,y_nan_idxs);
+
+ get_ticks(end_tick);
+ out.acc_scorrect += ticks2sec(end_tick,start_tick2);
for (int i = 0; i < a_iters; i++)
@@ -648,7 +688,8 @@
if (!params.ForceCheck && y_iters <= 2 &&
( (a_iters >= 10 && (i%(a_iters/(10))) == 0) || (a_iters < (10)) ))
{
- cout << "*" << flush;
+ cout << AIOfile.io_overhead << flush;
+ AIOfile.io_overhead = "*";
}
get_ticks(start_tick2);
@@ -663,28 +704,90 @@
replace_nans(ar_nan_idxs, a_block_size*r, backupAR, n , 1);
+ //sumSquares(backupAR,a_block_size*r,n,ssAR,0);
+
replace_nans_avgs(a_block_size, backupAR, n, r, ar_nan_idxs);
- sumSquares(backupAR,a_block_size*r,n,ssAR,ar_nan_idxs);
+ get_ticks(end_tick);
+ out.acc_other += ticks2sec(end_tick,start_tick2);
+ AR = backupAR;
+ get_ticks(start_tick2);
+
+ //! Generate Str
+ cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans,
+ l, r * a_block_size, n, 1.0, AL,
+ n, AR, n, 0.0, Str, l);//!
+
+ get_ticks(end_tick);
+ out.acc_str += ticks2sec(end_tick,start_tick2);
- //cout << "mb:" << a_block_size << " ";
- //matlab_print_matrix("ARb",n,a_block_size*r,backupAR);
+ get_ticks(start_tick2);
- //replace_with_zeros(al_nan_idxs, backupAR, n, r, a_block_size);
+ blas_set_num_threads(1);
+ omp_set_num_threads(max_threads);
+ #pragma omp parallel default(shared)
+ {
+
+ #pragma omp for nowait
+ for (int ii= 0; ii < a_block_size; ii++)
+ {
+ //! Generate Sbr
+ cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans,
+ r, r, n, 1.0, &AR[ii*r*n], n,
+ &AR[ii * r * n], n, 0.0,
+ &Sbr[ii * r * r], r);
+ float* S2 = &S2global[p*p*omp_get_thread_num()];
- copy_vec(backupAR, AR, n * r * a_block_size);
+ //! build S for X'X-1 for fake "variance of x"
+ build_S(S2, Stl, &Str[ii*r*l], &Sbr[ii*r*r], l, r);
+
+
+ info = LAPACKE_spotrf(STORAGE_TYPE,'U',p,S2,p);info++;
+
+ info = LAPACKE_spotri(STORAGE_TYPE,'U',p,S2,p);info++;
+
+ for (int h= 0; h < p; h++)
+ ssA[ii*p+h] = S2[h*p+h];//diagonal elements ARE the variance
+
+
+ }
+
+
+
+
+ }
+
+ blas_set_num_threads(max_threads);
+
+
+
get_ticks(end_tick);
- //out.acc_other += ticks2sec(end_tick,start_tick2);
+ out.acc_sbr += ticks2sec(end_tick,start_tick2 );
+
+ get_ticks(start_tick2);
+
+ generate_correction_missings_Str(y_block_size,a_block_size,n,l,r,AL,AR,Str_corr,Sbr_corr,y_nan_idxs);
+
+ get_ticks(end_tick);
+ out.acc_scorrect += ticks2sec(end_tick,start_tick2);
+
+
+
+
+
+ //cout << "mb:" << a_block_size << " ";
+
+ //matlab_print_matrix("ARb",n,a_block_size*r,backupAR);
+
-
get_ticks(start_tick2);
//! Ay_bot = AR'*Y
@@ -701,52 +804,8 @@
{
//int thread_id = 0 * omp_get_thread_num();//so far singel thread version only
-
-
//cout << "y_block_size:" << y_block_size << " a_block_size:" << a_block_size <<" r:" << r << "\n";
-
- get_ticks(start_tick2);
-
- copy_vec(backupAL, AL, n * l);//try to remove!
-
- replace_with_zeros(&y_nan_idxs[jj], AL, n, l, 1);
-
-
- get_ticks(end_tick);//2%
- out.acc_other += ticks2sec(end_tick,start_tick2);
-
- get_ticks(start_tick2);
-
- //! Generate Stl
- cblas_ssyrk(CblasColMajor, CblasUpper, CblasTrans,
- l, n, 1.0, AL, lda, 0.0, Stl, l);
-
- get_ticks(end_tick);
- out.acc_stl += ticks2sec(end_tick,start_tick2);
-
-
- get_ticks(start_tick2);
- copy_vec(backupAR,AR, n*r*a_block_size);//!10%//try to remove!
- replace_with_zeros(&y_nan_idxs[jj], AR, n, r, a_block_size);
- get_ticks(end_tick);
- out.acc_other += ticks2sec(end_tick,start_tick2);
-
- get_ticks(start_tick2);
-
- //! Generate Str
- cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans,
- l, r * a_block_size, n, 1.0, AL,
- n, AR, n, 0.0, Str, l);//!45
-
- get_ticks(end_tick);
- out.acc_str += ticks2sec(end_tick,start_tick2);
-
-
-
- //type_precision Sbr[r * r * a_block_size*1];
- //type_precision Ay[p * a_block_size*1];
-
blas_set_num_threads(1);
omp_set_num_threads(max_threads);
@@ -757,60 +816,56 @@
for (int ii= 0; ii < a_block_size; ii++)
{
// //cout << omp_get_thread_num() << endl << flush;
-//
-// get_ticks(start_tick2);
-//
-// copy_vec(&backupAR[ii*r*n],&AR[ii*r*n], n*r);//!10%//try to remove!
-// replace_with_zeros(&y_nan_idxs[jj], &AR[ii*r*n], n, r, 1);
-//
-// get_ticks(end_tick);
-// out.acc_other += ticks2sec(end_tick,start_tick2);
-
-
- get_ticks(start_tick2);
-
- //! Generate Sbr
- cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans,
- r, r, n, 1.0, &AR[ii*r*n], n,
- &AR[ii * r * n], n, 0.0,
- &Sbr[ii * r * r], r);
+ cputime_type start_tick2;
-
- get_ticks(end_tick);
- out.acc_sbr += ticks2sec(end_tick,start_tick2 );
-
get_ticks(start_tick2);
- copy_vec(&Ay_top[j*l*y_block_size+jj*l], &Ay[ii*p], l);
+ copy_vec(&Ay_top[j*l*y_block_size+jj*l], &Ay[jj*p*a_block_size+ ii*p], l);
copy_vec(&Ay_bot[ii*r + jj*r*a_block_size],
- &Ay[l+ii*p], r);
+ &Ay[jj*p*a_block_size+l+ii*p], r);
- type_precision* S2 = &S2global[p*p*omp_get_thread_num()];
+ float* S2 = &S2global[p*p*omp_get_thread_num()];
//! Rebuild S
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1832
More information about the Genabel-commits
mailing list