Commit df3ad671 authored by Pilar Cossio's avatar Pilar Cossio
Browse files

Cross-Correlation Optimization and more options

parent a0cc70c4
......@@ -70,6 +70,11 @@ namespace po = boost::program_options;
using namespace std;
/* For dvl nodes in hydra with problem in boost
namespace std {
typedef decltype(nullptr) nullptr_t;
}*/
// A helper function of Boost
template<class T>
ostream& operator<<(ostream& os, const vector<T>& v)
......@@ -265,8 +270,12 @@ int bioem::configure(int ac, char* av[])
cout << "Remark: look at file COORDREAD to confirm that the Model coordinates are correct\n";
if (DebugOutput >= 2 && mpi_rank == 0) printf("Reading Input Data %f\n", timer.GetCurrentElapsedTime());
}
if(param.param_device.writeCC && mpi_size>1){
cout << "Exiting::: WRITE CROSS-CORRELATION ONLY VAILD FOR 1 MPI PROCESS\n";
exit(1);
}
}
#ifdef WITH_MPI
if (mpi_size > 1)
{
......@@ -395,13 +404,15 @@ int bioem::run()
// ****************** Declarying class of Probability Pointer *************************
if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
// Inizialzing Probabilites to zero and constant to -Infinity
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
pProbMap.Total = 0.0;
pProbMap.Constoadd = -9999999;
pProbMap.Constoadd = -FLT_MAX; //Problem if using double presicion
if (param.param_device.writeAngles)
{
for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
......@@ -409,19 +420,26 @@ int bioem::run()
bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
pProbAngle.forAngles = 0.0;
pProbAngle.ConstAngle = -99999999;
pProbAngle.ConstAngle = -FLT_MAX;
}
}
if (param.param_device.writeCC)
if (param.param_device.writeCC)
{ int cc=0;
for (int cent_x = 0; cent_x < param.param_device.NumberPixels -param.param_device.CCdisplace ; cent_x = cent_x + param.param_device.CCdisplace)
for (int cent_x = 0; cent_x < param.param_device.NumberPixels; cent_x = cent_x + param.param_device.CCdisplace)
{
for (int cent_y = 0; cent_y < param.param_device.NumberPixels - param.param_device.CCdisplace ; cent_y = cent_y + param.param_device.CCdisplace)
for (int cent_y = 0; cent_y < param.param_device.NumberPixels; cent_y = cent_y + param.param_device.CCdisplace)
{
bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
pProbCC.forCC = 0.0;
cc++;
//Debuggin:: cout << iRefMap << " " << cc << " " << cent_x << " " << cent_y << "\n";
if(!param.param_device.CCwithBayes) {
pProbCC.forCC=-FLT_MAX;
}else {
pProbCC.forCC = 0.0;
pProbCC.ConstCC=-FLT_MAX;
}
cc++;
}
}
}
......@@ -626,30 +644,46 @@ int bioem::run()
if (mpi_rank == 0)
{
ofstream angProbfile;
// Output for Angle Probability File
ofstream angProbfile;
if(param.param_device.writeAngles)
{
angProbfile.open ("ANG_PROB");
angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
angProbfile <<" RefMap: MapNumber ; alpha - beta - gamma - log Probability\n" ;
angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
}
// Output for Cross Correlation File
ofstream ccProbfile;
if(param.param_device.writeCC)
{
ccProbfile.open ("CROSS_CORRELATION");
ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
ccProbfile <<" RefMap: MapNumber ; Pixel x - Pixel y - Cross-Correlation \n";
ccProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
}
ofstream outputProbFile;
outputProbFile.open ("Output_Probabilities");
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n";
outputProbFile << " RefMap: MapNumber ; LogProb: natural logarithm of posterior Probability ; Constant: Numerical Const. for adding Probabilities \n";
// Output for Standard Probability
ofstream outputProbFile;
outputProbFile.open ("Output_Probabilities");
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";
if(param.writeCTF) outputProbFile << " RefMap: MapNumber ; CTFMaxParm: defocus - b-factor (ref. Penzeck 2010)\n";
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";
if(param.writeCTF) outputProbFile << " RefMap: MapNumber ; CTFMaxParm: defocus - b-factor (ref. Penzeck 2010)\n";
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";
// Loop over reference maps
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
// **** Total Probability ***
bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
//Controll for Value of Total Probability
// 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";
outputProbFile << "RefMap: " << iRefMap << " Maximizing Param: ";
......@@ -667,17 +701,22 @@ int bioem::run()
outputProbFile << pProbMap.max.max_prob_norm << " [ ] " ;
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 << "RefMap: " << iRefMap << "Most probabily Numerical Constant Out of Float-point range: " << pProbMap.Constoadd << "\n"; }
// Writing out CTF parameters if requiered
if(param.writeCTF){
myfloat_t denomi;
denomi = param.CtfParam[pProbMap.max.max_prob_conv].pos[1] * param.CtfParam[pProbMap.max.max_prob_conv].pos[1] + param.CtfParam[pProbMap.max.max_prob_conv].pos[2] * param.CtfParam[pProbMap.max.max_prob_conv].pos[2];
outputProbFile << "RefMap: " << iRefMap << " CTFMaxParam: ";
outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel << " [micro-m]; ";
outputProbFile << "2*(" << sqrt(4*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi) << ")² [1./A²] \n";
}
// *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap);
myfloat_t denomi;
denomi = param.CtfParam[pProbMap.max.max_prob_conv].pos[1] * param.CtfParam[pProbMap.max.max_prob_conv].pos[1] +
param.CtfParam[pProbMap.max.max_prob_conv].pos[2] * param.CtfParam[pProbMap.max.max_prob_conv].pos[2];
outputProbFile << "RefMap: " << iRefMap << " CTFMaxParam: ";
outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel << " [micro-m]; ";
outputProbFile << "2*(" << sqrt(4*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi) << ")² [1./A²] \n";
}
// Writing Individual Angle probabilities
if(param.param_device.writeAngles)
{
for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
......@@ -687,8 +726,10 @@ int bioem::run()
angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << log(pProbAngle.forAngles) + pProbAngle.ConstAngle + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
}
}
if(param.param_device.writeCC)
{
// Writing Cross-Correlations if requiered
if(param.param_device.writeCC){
int cc=0;
int halfPix;
int rx=0;
......@@ -698,20 +739,22 @@ int bioem::run()
halfPix = param.param_device.NumberPixels / 2 ;
// Ordering the centers of the Cross Correlation
for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels -param.param_device.CCdisplace ; rx = rx + param.param_device.CCdisplace)
for (int rx = 0; rx < param.param_device.NumberPixels ; rx++)
{
for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels - param.param_device.CCdisplace ; ry = ry + param.param_device.CCdisplace)
for (int ry = 0; ry < param.param_device.NumberPixels ; ry++)
{
localcc[ rx * param.param_device.NumberPixels + ry ] = 0.0;
}
}
for (int cent_x = 0; cent_x < param.param_device.NumberPixels -param.param_device.CCdisplace ; cent_x = cent_x + param.param_device.CCdisplace)
for (int cent_x = 0; cent_x < param.param_device.NumberPixels ; cent_x = cent_x + param.param_device.CCdisplace)
{
for (int cent_y = 0; cent_y < param.param_device.NumberPixels - param.param_device.CCdisplace ; cent_y = cent_y + param.param_device.CCdisplace)
for (int cent_y = 0; cent_y < param.param_device.NumberPixels ; cent_y = cent_y + param.param_device.CCdisplace)
{
//localcc[ rx * param.param_device.NumberPixels + ry ] = 0.0;
bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
// Applying Periodic boundary conditions to the CC
if(cent_x < halfPix && cent_y < halfPix){
// ccProbfile << " " << iRefMap << " " << (myfloat_t) halfPix - cent_x << " " << halfPix - cent_y << " " << pProbCC.forCC <<"\n";
rx = halfPix - cent_x;
......@@ -729,20 +772,36 @@ int bioem::run()
rx = 3 * halfPix - cent_x;
ry = 3 * halfPix - cent_y;}
// cout << " TT " << cent_x << " " << rx << " " << cent_y << " " << ry << " " << pProbCC.forCC << "\n";
localcc[ rx * param.param_device.NumberPixels + ry ] = pProbCC.forCC;
if(!param.param_device.CCwithBayes){
localcc[ rx * param.param_device.NumberPixels + ry ] = pProbCC.forCC;
}else{
localcc[ rx * param.param_device.NumberPixels + ry ] = log(pProbCC.forCC)+pProbCC.ConstCC;
}
cc++;
}
// ccProbfile << "\n";
}
for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels -param.param_device.CCdisplace ; rx = rx + param.param_device.CCdisplace)
{
for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels - param.param_device.CCdisplace ; ry = ry + param.param_device.CCdisplace)
{
ccProbfile << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
}
ccProbfile << "\n";
}
if(!param.ignoreCCoff){
for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx = rx + param.param_device.CCdisplace)
{
for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry = ry + param.param_device.CCdisplace)
{
ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
// cout << " cc " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] <<"\n" ;
}
ccProbfile << "\n";
}
}else{
for (int rx = param.param_device.CCdisplace; rx < param.param_device.NumberPixels ; rx++)
{
for (int ry = param.param_device.CCdisplace; ry < param.param_device.NumberPixels ; ry++)
{
if(localcc[ rx * param.param_device.NumberPixels + ry ]!=0.0) ccProbfile << "RefMap: "<< iRefMap << " " << rx << " " << ry << " " << localcc[ rx * param.param_device.NumberPixels + ry ] << "\n" ;
}
ccProbfile << "\n";
}
}
}
}
......@@ -807,9 +866,9 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t
#ifdef PILAR_DEBUG
if (param.param_device.writeCC)
{ int cc=0;
for (int cent_x = 0; cent_x < param.param_device.NumberPixels -param.param_device.CCdisplace ; cent_x = cent_x + param.param_device.CCdisplace)
for (int cent_x = 0; cent_x < param.param_device.NumberPixels ; cent_x = cent_x + param.param_device.CCdisplace)
{
for (int cent_y = 0; cent_y < param.param_device.NumberPixels - param.param_device.CCdisplace ; cent_y = cent_y + param.param_device.CCdisplace)
for (int cent_y = 0; cent_y < param.param_device.NumberPixels ; cent_y = cent_y + param.param_device.CCdisplace)
{
cout << "CHECKCC " << " " << cent_x << " " << cent_y <<" " << lCC[cent_x * param.param_device.NumberPixels + cent_y] / (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels ) << "\n";
cc++;
......
......@@ -77,7 +77,7 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
if(param.debugterm){
// printf("Separate cc: %f c: %f oo: %f o: %f co: %f logP: %f\n",sumsquare,sum,sumsquareref,sumref,crossproMapConv, logpro );
// printf("Separate cc: %f c: %f oo: %f o: %f co: %f logP: %f\n",sumsquare,sum,sumsquareref,sumref,crossproMapConv, logpro );
}
return(logpro);
}
......@@ -109,6 +109,10 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
}
pProbMap.Total += exp(logpro - pProbMap.Constoadd);
if(param.debugterm){
printf("Separate Ptot: %f Const: %f logProb %f \n",pProbMap.Total,pProbMap.Constoadd,logpro);
}
if (param.writeAngles)
{
bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
......@@ -151,32 +155,45 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient,
if (param.writeCC)
{
int cc=0;
for (int cent_x = 0; cent_x < param.NumberPixels -param.CCdisplace ; cent_x = cent_x + param.CCdisplace)
for (int cent_x = 0; cent_x < param.NumberPixels ; cent_x = cent_x + param.CCdisplace)
{
for (int cent_y = 0; cent_y < param.NumberPixels - param.CCdisplace ; cent_y = cent_y + param.CCdisplace)
for (int cent_y = 0; cent_y < param.NumberPixels ; cent_y = cent_y + param.CCdisplace)
{
bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
myfloat_t ttmp,ttmp2;
ttmp2 = (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels);
if(not param.flipped){
//Here we are inverting the sign of the cross-correlation for the images that are not flipped
ttmp=-ttmp2; }
else{ ttmp=ttmp2; }
if(!param.CCwithBayes){
// Storing Only the Maximum both for flipped and not flipped
if(pProbCC.forCC < ttmp) pProbCC.forCC = ttmp;
}else {
// Storing the Cross-correlation with Bayesian formalism
if(pProbCC.ConstCC < ttmp)
{
pProbCC.forCC = pProbCC.forCC * exp(-ttmp + pProbCC.ConstCC);
pProbCC.ConstCC = ttmp;
}
pProbCC.forCC += exp(ttmp - pProbCC.ConstCC);
bioem_Probability_cc& pProbCC = pProb.getProbCC(iRefMap, cc);
myfloat_t ttmp;
ttmp= (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels);
if(not param.flipped){
// Storing the lowest cross correlations
if(pProbCC.forCC > ttmp) pProbCC.forCC = ttmp;
}
else{
//Storing the highest cross correlations
if(pProbCC.forCC < ttmp) pProbCC.forCC = ttmp;
}
// pProbCC.forCC = ttmp;
cc++;
}
// printf("Separate Ptot: %f Const: %f logProb %f \n",pProbCC.forCC,pProbCC.ConstCC,ttmp);
cc++;
}
}
}
}
template <int GPUAlgo, class RefT>
......
......@@ -105,7 +105,8 @@ public:
class bioem_Probability_cc
{
public:
myfloat_t forCC;
myfloat_t forCC;
myfloat_t ConstCC;
};
class bioem_Probability
......@@ -139,8 +140,8 @@ public:
{
ptr_map = (bioem_Probability_map*) ptr;
ptr_angle = (bioem_Probability_angle*) (&ptr_map[nMaps]);
ptr_cc = (bioem_Probability_cc*) (&ptr_angle[nMaps * nAngles]);
// ptr_cc = (bioem_Probability_cc*) (&ptr_angle[nMaps * 2000 ]);
// ptr_cc = (bioem_Probability_cc*) (&ptr_angle[nMaps * nAngles]);
ptr_cc = (bioem_Probability_cc*) (&ptr_map[nMaps]);
}
};
......
......@@ -27,6 +27,7 @@ public:
bool flipped;
bool debugterm;
int CCdisplace;
bool CCwithBayes;
};
......@@ -43,6 +44,7 @@ public:
void PrepareFFTs();
bool doaaradius;
bool writeCTF;
bool ignoreCCoff;
myfloat_t elecwavel;
bioem_param_device param_device;
......
This diff is collapsed.
......@@ -77,8 +77,10 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
param_device.flipped=false;
param_device.debugterm=false;
param_device.writeCC=false;
param_device.CCwithBayes=false;
writeCTF=false;
elecwavel=0.0;
ignoreCCoff=false;
//Assigning name of angles
//Storing angle file name if existing.
......@@ -241,11 +243,21 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
if (param_device.CCdisplace < 0 ) { cout << "*** Error: Negative CROSSCOR_DISPLACE "; exit(1);}
cout << "Writing Cross Correlation Displacement " << param_device.CCdisplace << " \n";
}
else if (strcmp(token, "CROSSCOR_WITHBAYESIAN") == 0)
{
param_device.CCwithBayes=true;
cout << "Using Bayesian Analysis to write Cross Correlation \n";
}
else if (strcmp(token, "FLIPPED") == 0) //Key word if images are flipped for cross-correlation
{
param_device.flipped = true;
cout << "Micrograph Flipped Intensities \n";
}
else if (strcmp(token, "IGNORE_CROSSCORR_OFFSET") == 0) //Key word if images are flipped for cross-correlation
{
ignoreCCoff = true;
cout << "Ignoring Cross-Correlation offset \n";
}
else if (strcmp(token, "DO_AA_RADIUS") == 0) //If projecting CA with amino-acid radius
{
doaaradius = true;
......@@ -296,6 +308,10 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
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);};
if( !param_device.writeCC && param_device.CCwithBayes ){ cout << "**** INPUT MISSING: WRITE_CROSSCOR keyword \n"; exit(1);}
if( param_device.writeCC) {if(!param_device.CCwithBayes ){ cout << "Remark:: Not Using Bayesian method to store Cross-Correlation.\n Only Printing out Maximum\n";}
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;
......@@ -645,7 +661,7 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
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=(myfloat_t) (param_device.NumberPixels * param_device.NumberPixels) / (myfloat_t) ( param_device.CCdisplace * param_device.CCdisplace );
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