Skip to content

GitLab

B
BioEM
Commit 987c9a7d authored by Pilar Cossio's avatar Pilar Cossio
Browse files
Options

Cleaning/debuging input

parent f83b2d54
Hide whitespace changes
Inline Side-by-side
Showing with 111 additions and 59 deletions
bioem.cpp
View file @ 987c9a7d
......@@ -115,6 +115,7 @@ int bioem::configure(int ac, char* av[])
RefMap.readMRC = false;
RefMap.readMultMRC = false;
param.notuniformangles=false;
yesoutfilename=false;
// *************************************************************************************
cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
......@@ -135,6 +136,7 @@ int bioem::configure(int ac, char* av[])
("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
("DumpMaps", "(Optional) Dump maps after they were red from maps file")
("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
("OutputFile", po::value< std::string>(), "(Optional) For changing the outputfile name")
("help", "(Optional) Produce help message")
;
......@@ -150,6 +152,7 @@ int bioem::configure(int ac, char* av[])
p.add("PrintBestCalMap",-1);
p.add("DumpMaps", -1);
p.add("LoadMapDump", -1);
p.add("OutputFile",-1);
po::variables_map vm;
po::store(po::command_line_parser(ac, av).
......@@ -195,6 +198,12 @@ int bioem::configure(int ac, char* av[])
Inputanglefile = vm["ReadOrientation"].as< std::string >();
param.notuniformangles=true;
}
if (vm.count("OutputFile"))
{
OutfileName = vm["OutputFile"].as< std::string >();
cout << "Writing OUTPUT to: " << vm["OutputFile"].as< std::string >() << "\n";
yesoutfilename=true;
}
if (vm.count("PrintBestCalMap"))
{
cout << "Reading Euler Angles from file: "
......@@ -655,7 +664,8 @@ int bioem::run()
{
angProbfile.open ("ANG_PROB");
angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
angProbfile <<" RefMap: MapNumber ; alpha - beta - gamma - log Probability\n" ;
if(!param.doquater){ angProbfile <<" RefMap: MapNumber ; alpha[rad] - beta[rad] - gamma[rad] - log Probability\n" ;}
else { angProbfile <<" RefMap: MapNumber ; q1 - q2 -q3 - log Probability\n" ;};
angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
}
// Output for Cross Correlation File
......@@ -670,10 +680,14 @@ int bioem::run()
// Output for Standard Probability
ofstream outputProbFile;
outputProbFile.open ("Output_Probabilities");
if(!yesoutfilename)OutfileName="Output_Probabilities";
outputProbFile.open (OutfileName.c_str());
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n";
outputProbFile << " RefMap: MapNumber ; LogProb: natural logarithm of posterior Probability ; Constant: Numerical Const. for adding Probabilities \n";
outputProbFile << " RefMap: MapNumber ; Maximizing Param: MaxLogProb - alpha - beta - gamma - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
outputProbFile << "Notation= RefMap: MapNumber ; LogProb: natural logarithm of posterior Probability ; Constant: Numerical Const. for adding Probabilities \n";
if(!param.doquater){
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - alpha[rad] - beta[rad] - gamma[rad] - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";}
else {
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";}
if(param.writeCTF) outputProbFile << " RefMap: MapNumber ; CTFMaxParm: defocus - b-factor (ref. Penzeck 2010)\n";
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";
......@@ -685,7 +699,7 @@ int bioem::run()
bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
//Controll for Value of Total Probability
// cout << pProbMap.Total << " " << pProbMap.Constoadd << " " << FLT_MAX <<" " << log(FLT_MAX) << "\n";
// cout << pProbMap.Total << " " << pProbMap.Constoadd << " " << FLT_MAX <<" " << log(FLT_MAX) << "\n";
if(pProbMap.Total>1.e-38){
outputProbFile << "RefMap: " << iRefMap << " LogProb: " << log(pProbMap.Total) + pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant: " << pProbMap.Constoadd << "\n";
......@@ -694,9 +708,9 @@ int bioem::run()
// *** Param that maximize probability****
outputProbFile << (pProbMap.Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[0] << " [rad] ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[1] << " [rad] ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[2] << " [rad] ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[0] << " [ ] ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[1] << " [ ] ";
outputProbFile << param.angles[pProbMap.max.max_prob_orient].pos[2] << " [ ] ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[0] << " [ ] ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << " [1./A²] ";
outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " [1./A²] ";
......@@ -706,7 +720,7 @@ int bioem::run()
outputProbFile << pProbMap.max.max_prob_mu << " [ ] ";
outputProbFile << "\n";
}else{
outputProbFile << "PROBLEM!! with Map " << iRefMap << "Numerical Integrated Probability = 0.0; Try re-normalizing experimental map\n";
outputProbFile << "PROBLEM!! with Map " << iRefMap << "Numerical Integrated Probability = 0.0; \n - Try re-normalizing experimental map \n-looking at your model density. Maybe RADIUS IS LESS THAN PIXEL SIZE \n If so try KEYWORD NO_PROJECT_RADIUS in parameter inputfile";
outputProbFile << "RefMap: " << iRefMap << "Most probabily Numerical Constant Out of Float-point range: " << pProbMap.Constoadd << "\n"; }
// Writing out CTF parameters if requiered
......
include/bioem.h
View file @ 987c9a7d
......@@ -36,6 +36,10 @@ public:
bioem_Probability pProb;
string OutfileName;
bool yesoutfilename;
protected:
virtual int deviceInit();
virtual int deviceStartRun();
......
include/param.h
View file @ 987c9a7d
......@@ -29,6 +29,7 @@ public:
int CCdisplace;
bool CCwithBayes;
};
class bioem_param
......@@ -68,6 +69,7 @@ public:
int GridPointsQuatern;
bool doquater;
myfloat_t voluang;
bool notuniformangles;
int NotUn_angles;
......@@ -80,13 +82,16 @@ public:
// Grid sampling for the convolution kernel
//ENVELOPE
myfloat_t startGridEnvelop;
myfloat_t endGridEnvelop;
int numberGridPointsEnvelop;
myfloat_t gridEnvelop;
//CTF=Amp*cos(phase*x)-sqrt(1-Amp**2)*sin(phase*x)
myfloat_t startGridCTF_phase;
myfloat_t endGridCTF_phase;
int numberGridPointsCTF_phase;
myfloat_t gridCTF_phase;
myfloat_t startGridCTF_amp;
myfloat_t endGridCTF_amp;
int numberGridPointsCTF_amp;
myfloat_t gridCTF_amp;
// Others
......
model.cpp
View file @ 987c9a7d
......@@ -14,6 +14,7 @@
#include <cstring>
#include "model.h"
#include "param.h"
using namespace std;
......@@ -27,7 +28,7 @@ bioem_model::~bioem_model()
if (points) free(points);
}
int bioem_model::readModel(const char* filemodel)
int bioem_model::readModel(bioem_param& param, const char* filemodel)
{
// **************************************************************************************
// ***************Reading reference Models either PDB or x,y,z,r,d format****************
......@@ -211,6 +212,8 @@ int bioem_model::readModel(const char* filemodel)
//Moving to Model to its center of density mass:
myfloat3_t r_cm;
if(not(param.nocentermass)){ //by default it is normally done
for(int n = 0; n < 3; n++)r_cm.pos[n] = 0.0;
for(int n = 0; n < nPointsModel; n++)
......@@ -230,7 +233,7 @@ int bioem_model::readModel(const char* filemodel)
points[n].point.pos[2] -= r_cm.pos[2];
}
}
return(0);
}
......
param.cpp
View file @ 987c9a7d
......@@ -149,8 +149,8 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
}
else if (strcmp(token, "QUATERNIONS") == 0)
{
token = strtok(NULL, " ");
GridPointsQuatern = int(atoi(token));
if(!notuniformangles){token = strtok(NULL, " ");
GridPointsQuatern = int(atoi(token));}
doquater= true;
cout << "Gridpoints Quaternions " << GridPointsQuatern << "\n";
......@@ -172,12 +172,13 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
cout << "Start Envelope " << startGridEnvelop << "\n";
yesSTEnv = true ;
}
else if (strcmp(token, "GRIDSPACE_ENVELOPE") == 0)
else if (strcmp(token, "END_ENVELOPE") == 0)
{
token = strtok(NULL, " ");
gridEnvelop = atof(token);
if (gridEnvelop < 0 ) { cout << "*** Error: Negative GRIDSPACE_ENVELOPE "; exit(1);}
cout << "Grid spacing Envelope " << gridEnvelop << "\n";
endGridEnvelop = atof(token);
// gridEnvelop = atof(token);
if ( endGridEnvelop < 0 ) { cout << "*** Error: Negative END_ENVELOPE "; exit(1);}
cout << "End Envelope range " << endGridEnvelop << "\n";
yesGSPEnv = true ;
}
else if (strcmp(token,"GRIDPOINTS_PSF_PHASE")==0)
......@@ -194,11 +195,11 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
cout << "Start PSF " << startGridCTF_phase << "\n";
yesSTpha = true ;
}
else if (strcmp(token,"GRIDSPACE_PSF_PHASE")==0)
else if (strcmp(token,"END_PSF_PHASE")==0)
{
token = strtok(NULL," ");
gridCTF_phase=atof(token);
cout << "Grid Space PSF " << gridCTF_phase << "\n";
endGridCTF_phase=atof(token);
cout << "End Grid phase PSF " << endGridCTF_phase << "\n";
yesGSPpha = true ;
}
else if (strcmp(token,"GRIDPOINTS_PSF_AMP")==0)
......@@ -215,11 +216,11 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
cout << "Start PSF " << startGridCTF_amp << "\n";
yesSTamp = true ;
}
else if (strcmp(token,"GRIDSPACE_PSF_AMP")==0)
else if (strcmp(token,"END_PSF_AMP")==0)
{
token = strtok(NULL," ");
gridCTF_amp=atof(token);
cout << "Grid Space PSF " << gridCTF_amp << "\n";
endGridCTF_amp=atof(token);
cout << "End amplitud grid " << endGridCTF_amp << "\n";
yesGSPamp = true ;
}
else if (strcmp(token, "MAX_D_CENTER") == 0)
......@@ -333,9 +334,9 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
if( not ( yesSTEnv ) ) { cout << "**** INPUT MISSING: Please provide START_ENVELOPE \n" ; exit (1);};
if( not ( yesSTamp ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_AMP \n" ; exit (1);};
if( not ( yesSTpha ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_PHASE \n" ; exit (1);};
if( not ( yesGSPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_AMP \n" ; exit (1);};
if( not ( yesGSPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_ENVELOPE \n" ; exit (1);};
if( not ( yesGSPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_PHASE \n" ; exit (1);};
if( not ( yesGSPamp ) ) { cout << "**** INPUT MISSING: Please provide END_PSF_AMP \n" ; exit (1);};
if( not ( yesGSPEnv ) ) { cout << "**** INPUT MISSING: Please provide END_ENVELOPE \n" ; exit (1);};
if( not ( yesGSPpha ) ) { cout << "**** INPUT MISSING: Please provide END_PSF_PHASE \n" ; exit (1);};
if( not ( yesMDC ) ) { cout << "**** INPUT MISSING: Please provide MAX_D_CENTER \n" ; exit (1);};
if( not ( yesGCen ) ) { cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ; exit (1);};
if( param_device.writeCC && param_device.CCdisplace < 1 ){ cout << "**** INPUT MISSING: Please provide CROSSCOR_DISPLACE \n" ; exit (1);};
......@@ -344,27 +345,8 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
if(param_device.flipped){ cout << "Remark:: Micrographs are Flipped = Particles are white\n";} else { cout << "Remark:: Micrographs are NOT Flipped = Particles are dark\n";}
}
//if only one grid point for PSF kernel:
if( (myfloat_t) numberGridPointsCTF_amp == 1 ) gridCTF_amp = startGridCTF_amp;
if( (myfloat_t) numberGridPointsCTF_phase == 1 ) gridCTF_phase = startGridCTF_phase;
if( (myfloat_t) numberGridPointsEnvelop == 1 ) gridEnvelop = startGridEnvelop;
//More checks with input parameters
// Envelope should not have a standard deviation greater than Npix/2
if(sqrt(1./( (myfloat_t) numberGridPointsDisplaceCenter * gridEnvelop + startGridEnvelop))> float(param_device.NumberPixels)/2.0) {
cout << "MAX Standard deviation of envelope is larger than Allowed KERNEL Length\n";
exit(1);
}
// Envelop param should be positive
if( startGridCTF_amp < 0 || startGridCTF_amp > 1){
cout << "PSF Amplitud should be between 0 and 1\n";
exit(1);
}
if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) (numberGridPointsCTF_amp - 1) * gridCTF_amp + startGridCTF_amp > 1){
cout << "Value should be between 0 and 1\n" ;
exit(1);
}
param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
FFTMapSize = param_device.NumberPixels * param_device.NumberFFTPixels1D;
......@@ -592,7 +574,7 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
// **************** Routine that pre-calculates Euler angle grids **********************
// ************************************************************************************
myfloat_t voluang;
// myfloat_t voluang;
// Euler angle GRID
......@@ -834,19 +816,7 @@ if(!notuniformangles){
cout << "Analysis with Quaternions. Total number of quaternions " << nTotGridAngles << "\n";
}
// ********** Calculating normalized volumen element *********
param_device.volu = voluang * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
* gridCTF_phase * gridCTF_amp * gridEnvelop / ((2.f * (myfloat_t) param_device.maxDisplaceCenter+1.)) / (2.f * (myfloat_t) (param_device.maxDisplaceCenter + 1.)) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp )
/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop );
// cout << "VOLU " << param_device.volu << " " << gridCTF_amp << "\n";
// *** Number of total pixels***
param_device.Ntotpi = (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
param_device.NtotDist = (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1);
nTotCC = (int) ((myfloat_t) param_device.NumberPixels / (myfloat_t) param_device.CCdisplace + 1) * (int) ((myfloat_t) param_device.NumberPixels / (myfloat_t) param_device.CCdisplace + 1);
return(0);
}
......@@ -874,6 +844,48 @@ int bioem_param::CalculateRefCTF()
myfloat_t normctf;
gridCTF_amp = (endGridCTF_amp - startGridCTF_amp) / (myfloat_t) numberGridPointsCTF_amp;
gridCTF_phase = (endGridCTF_phase - startGridCTF_phase) / (myfloat_t) numberGridPointsCTF_phase;
gridEnvelop = (endGridEnvelop - startGridEnvelop) / (myfloat_t) numberGridPointsEnvelop;
//if only one grid point for PSF kernel:
if( (myfloat_t) numberGridPointsCTF_amp == 1 ) {
gridCTF_amp = startGridCTF_amp;}
else if ( (endGridCTF_amp - startGridCTF_amp) < 0. ){
cout << "Error: Interval of amplitude in PSF is Negative"; exit(1);
}
if( (myfloat_t) numberGridPointsCTF_phase == 1 ) {
gridCTF_phase = startGridCTF_phase;
}else if ( (endGridCTF_phase - startGridCTF_phase) < 0.){
cout << "Error: Interval of PHASE in PSF is Negative"; exit(1);
}
if( (myfloat_t) numberGridPointsEnvelop == 1 ) {
gridEnvelop = startGridEnvelop;
} else if ( (endGridEnvelop - startGridEnvelop) < 0.) {
cout << "Error: Interval of Envelope in PSF is Negative"; exit(1);
}
//More checks with input parameters
// Envelope should not have a standard deviation greater than Npix/2
if(sqrt(1./( (myfloat_t) numberGridPointsEnvelop * gridEnvelop + startGridEnvelop))> float(param_device.NumberPixels)/2.0) {
cout << "MAX Standard deviation of envelope is larger than Allowed KERNEL Length\n";
exit(1);
}
// Envelop param should be positive
if( startGridCTF_amp < 0 || startGridCTF_amp > 1){
cout << "Error: PSF Amplitud should be between 0 and 1\n";
exit(1);
}
if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) (numberGridPointsCTF_amp - 1) * gridCTF_amp + startGridCTF_amp > 1){
cout << "Error: alues of amplitud in PSF should be between 0 and 1\n" ;
exit(1);
}
for (int iamp = 0; iamp < numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
{
amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
......@@ -943,6 +955,20 @@ int bioem_param::CalculateRefCTF()
exit(1);
}
// ********** Calculating normalized volumen element *********
param_device.volu = voluang * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
* gridCTF_phase * gridCTF_amp * gridEnvelop / ((2.f * (myfloat_t) param_device.maxDisplaceCenter+1.)) / (2.f * (myfloat_t) (param_device.maxDisplaceCenter + 1.)) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp )
/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop );
// cout << "VOLU " << param_device.volu << " " << gridCTF_amp << "\n";
// *** Number of total pixels***
param_device.Ntotpi = (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
param_device.NtotDist = (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1);
nTotCC = (int) ((myfloat_t) param_device.NumberPixels / (myfloat_t) param_device.CCdisplace + 1) * (int) ((myfloat_t) param_device.NumberPixels / (myfloat_t) param_device.CCdisplace + 1);
return(0);
}
......
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