Commit e73c526d authored by Pilar Cossio's avatar Pilar Cossio

bioem.cpp

parent 7b6dd3b4
...@@ -86,7 +86,15 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -86,7 +86,15 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
if(strcmp(tmpm,"XXX")!=0) if(strcmp(tmpm,"XXX")!=0)
{ {
indifile=strline.c_str(); indifile=strline.c_str();
cout << "FILE Name: " << indifile << "\n";
size_t foundpos= strline.find(".mrc");
size_t endpos = strline.find_last_not_of(" \t");
if(foundpos > endpos){
cout << "Warining:::: .mrc extension NOT dectected in file name::" << filemap <<" \n";
cout << "Warining:::: Are you sure you want to read an MRC? \n";
}
//Reading Multiple MRC //Reading Multiple MRC
read_MRC(indifile,param); read_MRC(indifile,param);
} }
...@@ -99,6 +107,17 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -99,6 +107,17 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
} }
else else
{ {
string strfilename(filemap);
size_t foundpos= strfilename.find(".mrc");
size_t endpos = strfilename.find_last_not_of(" \t");
if(foundpos > endpos){
cout << "Warining:::: .mrc extension NOT dectected in file name::" << filemap <<" \n";
cout << "Warining:::: Are you sure you want to read an MRC? \n";
}
read_MRC(filemap,param); read_MRC(filemap,param);
cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n"; cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ; cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ;
...@@ -120,6 +139,8 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -120,6 +139,8 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
char tmpLine[512] = {0}; char tmpLine[512] = {0};
bool first=true; bool first=true;
int countpix;
while (!input.eof()) while (!input.eof())
{ {
input.getline(line, 511); input.getline(line, 511);
...@@ -135,9 +156,11 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -135,9 +156,11 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
first=false; first=false;
} }
if (strcmp(token, "PARTICLE") == 0) // to count the number of maps if (strcmp(token, "PARTICLE") == 0) // to count the number of maps
{ {
nummap++; nummap++;
countpix=0;
if (allocsize == 0) if (allocsize == 0)
{ {
allocsize = 64; allocsize = 64;
...@@ -162,7 +185,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -162,7 +185,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
{ {
int i, j; int i, j;
float z; float z;
char tmpVals[36] = {0}; char tmpVals[36] = {0};
strncpy (tmpVals, line, 8); strncpy (tmpVals, line, 8);
...@@ -171,14 +194,16 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -171,14 +194,16 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
strncpy (tmpVals, line + 8, 8); strncpy (tmpVals, line + 8, 8);
sscanf (tmpVals, "%d", &j); sscanf (tmpVals, "%d", &j);
strncpy (tmpVals, line + 16, 12); strncpy (tmpVals, line + 16, 16);
sscanf (tmpVals, "%f", &z); sscanf (tmpVals, "%f", &z);
//checking for Map limits //checking for Map limits
if(i > 0 && i - 1 < param.param_device.NumberPixels && j > 0 && j - 1 < param.param_device.NumberPixels) if(i > -1 && i < param.param_device.NumberPixels && j > -1 && j < param.param_device.NumberPixels)
{ {
maps[nummap * refMapSize + (i - 1) * numPixels + (j - 1)] = (myfloat_t) z; countpix++;
maps[nummap * refMapSize + i * numPixels + j] = (myfloat_t) z;
lasti = i; lasti = i;
lastj = j; lastj = j;
// cout << countpix << " " << param.param_device.NumberPixels*param.param_device.NumberPixels << "\n";
} }
else else
{ {
...@@ -187,7 +212,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -187,7 +212,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
} }
} }
} }
if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels ) if(lasti != param.param_device.NumberPixels-1 || lastj != param.param_device.NumberPixels-1 || countpix != param.param_device.NumberPixels*param.param_device.NumberPixels +1 )
{ {
cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n"; cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
exit(1); exit(1);
...@@ -197,20 +222,23 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) ...@@ -197,20 +222,23 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap); maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
cout << "Particle Maps read from Standard File \n"; cout << "Particle Maps read from Standard File \n";
if (param.dumpMap)
{
FILE* fp = fopen("maps.dump", "w+b");
if (fp == NULL)
{
cout << "Error opening dump file\n";
exit(1);
}
fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
fwrite(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
fclose(fp);
}
} }
//Dumping Maps
if (param.dumpMap)
{
FILE* fp = fopen("maps.dump", "w+b");
if (fp == NULL)
{
cout << "Error opening dump file\n";
exit(1);
}
fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
fwrite(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
fclose(fp);
}
if (getenv("BIOEM_DEBUG_NMAPS")) if (getenv("BIOEM_DEBUG_NMAPS"))
{ {
ntotRefMap = atoi(getenv("BIOEM_DEBUG_NMAPS")); ntotRefMap = atoi(getenv("BIOEM_DEBUG_NMAPS"));
...@@ -299,7 +327,8 @@ void bioem_Probability::init(size_t maps, size_t angles, size_t cc, bioem& bio) ...@@ -299,7 +327,8 @@ void bioem_Probability::init(size_t maps, size_t angles, size_t cc, bioem& bio)
nAngles = angles; nAngles = angles;
nCC = cc; nCC = cc;
ptr = bio.malloc_device_host(get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC)); ptr = bio.malloc_device_host(get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC));
cout << "Allocation #Maps " << maps << " #Angles " << angles << " #cross.cor " << cc << "== " << get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC)<< "\n"; cout << "Allocation #Maps " << maps << " #Angles " << angles << " #cross.cor " << cc << "\n";
//<< " == " << get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC)<< "\n";
set_pointers(); set_pointers();
} }
...@@ -483,7 +512,8 @@ int bioem_RefMap::read_MRC(const char* filename,bioem_param& param) ...@@ -483,7 +512,8 @@ int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
st2 = sqrt(st2 / float(nr * nc) - st * st); st2 = sqrt(st2 / float(nr * nc) - st * st);
for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){ for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){
maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] / st2 - st/st2; maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] / st2 - st/st2;
//if(nmap+ntotRefMap==300) cout << i << " " << j << " " << nmap+ntotRefMap << " " << Ref[nmap+ntotRefMap].points[i][j] << "\n"; //(nmap+ntotRefMap==300)
//cout <<"MAP:: " << i << " " << j << " " << maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] << "\n";
} }
} }
ntotRefMap += ns ; ntotRefMap += ns ;
......
...@@ -35,6 +35,8 @@ int bioem_model::readModel(const char* filemodel) ...@@ -35,6 +35,8 @@ int bioem_model::readModel(const char* filemodel)
ofstream exampleReadCoor; ofstream exampleReadCoor;
exampleReadCoor.open ("COORDREAD"); exampleReadCoor.open ("COORDREAD");
exampleReadCoor << "Text --- Number ---- x ---- y ---- z ---- radius ---- density\n";
int allocsize = 0; int allocsize = 0;
...@@ -53,6 +55,17 @@ int bioem_model::readModel(const char* filemodel) ...@@ -53,6 +55,17 @@ int bioem_model::readModel(const char* filemodel)
int numres = 0; int numres = 0;
NormDen = 0.0; NormDen = 0.0;
string strfilename(filemodel);
size_t foundpos= strfilename.find(".pdb");
size_t endpos = strfilename.find_last_not_of(" \t");
if(foundpos > endpos){
cout << "Warining:::: .pdb extension NOT dectected in file name \n";
cout << "Warining:::: Are you sure you want to read a PDB? \n";
}
// cout << " HERE " << filemodel ; // cout << " HERE " << filemodel ;
// for eachline in the file // for eachline in the file
while (!input.eof()) while (!input.eof())
...@@ -127,7 +140,7 @@ int bioem_model::readModel(const char* filemodel) ...@@ -127,7 +140,7 @@ int bioem_model::readModel(const char* filemodel)
points[numres].point.pos[0] = (myfloat_t) x; points[numres].point.pos[0] = (myfloat_t) x;
points[numres].point.pos[1] = (myfloat_t) y; points[numres].point.pos[1] = (myfloat_t) y;
points[numres].point.pos[2] = (myfloat_t) z; points[numres].point.pos[2] = (myfloat_t) z;
exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].density << "\n"; exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].radius << " " << points[numres].density << "\n";
numres++; numres++;
} }
} }
...@@ -182,7 +195,7 @@ int bioem_model::readModel(const char* filemodel) ...@@ -182,7 +195,7 @@ int bioem_model::readModel(const char* filemodel)
points[numres].radius = (myfloat_t) tmpval[3]; points[numres].radius = (myfloat_t) tmpval[3];
points[numres].density = (myfloat_t) tmpval[4]; points[numres].density = (myfloat_t) tmpval[4];
exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].density << "\n"; exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].radius << " " << points[numres].density << "\n";
NormDen += points[numres].density; NormDen += points[numres].density;
numres++; numres++;
} }
......
...@@ -77,7 +77,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) ...@@ -77,7 +77,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
param_device.flipped=false; param_device.flipped=false;
param_device.debugterm=false; param_device.debugterm=false;
param_device.writeCC=false; param_device.writeCC=false;
param_device.CCwithBayes=false; param_device.CCwithBayes=true;
writeCTF=false; writeCTF=false;
elecwavel=220; elecwavel=220;
ignoreCCoff=false; ignoreCCoff=false;
...@@ -236,16 +236,16 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) ...@@ -236,16 +236,16 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
param_device.writeCC = true; param_device.writeCC = true;
cout << "Writing CrossCorrelations \n"; cout << "Writing CrossCorrelations \n";
} }
else if (strcmp(token, "CROSSCOR_DISPLACE") == 0) else if (strcmp(token, "CROSSCOR_GRID_SPACE") == 0)
{ {
token = strtok(NULL, " "); token = strtok(NULL, " ");
param_device.CCdisplace=int(atoi(token)); param_device.CCdisplace=int(atoi(token));
if (param_device.CCdisplace < 0 ) { cout << "*** Error: Negative CROSSCOR_DISPLACE "; exit(1);} if (param_device.CCdisplace < 0 ) { cout << "*** Error: Negative CROSSCOR_DISPLACE "; exit(1);}
cout << "Writing Cross Correlation Displacement " << param_device.CCdisplace << " \n"; cout << "Writing Cross Correlation Displacement " << param_device.CCdisplace << " \n";
} }
else if (strcmp(token, "CROSSCOR_WITHBAYESIAN") == 0) else if (strcmp(token, "CROSSCOR_NOTBAYESIAN") == 0)
{ {
param_device.CCwithBayes=true; param_device.CCwithBayes=false;
cout << "Using Bayesian Analysis to write Cross Correlation \n"; cout << "Using Bayesian Analysis to write Cross Correlation \n";
} }
else if (strcmp(token, "FLIPPED") == 0) //Key word if images are flipped for cross-correlation else if (strcmp(token, "FLIPPED") == 0) //Key word if images are flipped for cross-correlation
...@@ -258,7 +258,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) ...@@ -258,7 +258,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
ignoreCCoff = true; ignoreCCoff = true;
cout << "Ignoring Cross-Correlation offset \n"; cout << "Ignoring Cross-Correlation offset \n";
} }
else if (strcmp(token, "DO_AA_RADIUS") == 0) //If projecting CA with amino-acid radius else if (strcmp(token, "PROJECT_RADIUS") == 0) //If projecting CA with amino-acid radius
{ {
doaaradius = true; doaaradius = true;
cout << "Projecting corresponding radius \n"; cout << "Projecting corresponding radius \n";
...@@ -300,8 +300,10 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) ...@@ -300,8 +300,10 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
// Checks for ALL INPUT // Checks for ALL INPUT
if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);}; if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);};
if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);}; if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);};
if( not ( yesGPal ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);}; if(!notuniformangles){
if( not ( yesGPbe ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ; exit (1);}; if( not ( yesGPal) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);};
if( not ( yesGPbe )) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ; exit (1);};
}
if( not ( yesGPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ; exit (1);}; if( not ( yesGPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ; exit (1);};
if( not ( yesGPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ; exit (1);}; if( not ( yesGPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ; exit (1);};
if( not ( yesGPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ; exit (1);}; if( not ( yesGPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ; exit (1);};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment