Commit 9582b044 authored by David Rohr's avatar David Rohr
Browse files

fix comment style

parent 1e3203f7
...@@ -47,13 +47,13 @@ bioem::~bioem() ...@@ -47,13 +47,13 @@ bioem::~bioem()
int bioem::configure(int ac, char* av[]) int bioem::configure(int ac, char* av[])
{ {
//************************************************************************************** // **************************************************************************************
//**** Configuration Routine using boost for extracting parameters, models and maps **** // **** Configuration Routine using boost for extracting parameters, models and maps ****
//************************************************************************************** // **************************************************************************************
//****** And Precalculating necessary grids, map crosscorrelations and kernels ******** // ****** And Precalculating necessary grids, map crosscorrelations and kernels ********
//************************************************************************************* // *************************************************************************************
/*** Inizialzing default variables ***/ // *** Inizialzing default variables ***
std::string infile, modelfile, mapfile; std::string infile, modelfile, mapfile;
Model.readPDB = false; Model.readPDB = false;
param.writeAngles = false; param.writeAngles = false;
...@@ -63,11 +63,11 @@ int bioem::configure(int ac, char* av[]) ...@@ -63,11 +63,11 @@ int bioem::configure(int ac, char* av[])
RefMap.readMultMRC = false; RefMap.readMultMRC = false;
//************************************************************************************* // *************************************************************************************
cout << "+++++++++++++ FROM COMMAND LINE +++++++++++\n\n"; cout << "+++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
//************************************************************************************* // *************************************************************************************
//********************* Command line reading input with BOOST ************************ // ********************* Command line reading input with BOOST ************************
try { try {
po::options_description desc("Command line inputs"); po::options_description desc("Command line inputs");
...@@ -164,23 +164,23 @@ int bioem::configure(int ac, char* av[]) ...@@ -164,23 +164,23 @@ int bioem::configure(int ac, char* av[])
return 1; return 1;
} }
//********************* Reading Parameter Input *************************** // ********************* Reading Parameter Input ***************************
// copying inputfile to param class // copying inputfile to param class
param.fileinput = infile.c_str(); param.fileinput = infile.c_str();
param.readParameters(); param.readParameters();
//********************* Reading Model Input ****************************** // ********************* Reading Model Input ******************************
// copying modelfile to model class // copying modelfile to model class
Model.filemodel = modelfile.c_str(); Model.filemodel = modelfile.c_str();
Model.readModel(); Model.readModel();
//********************* Reading Particle Maps Input ********************** // ********************* Reading Particle Maps Input **********************
//********* HERE: PROBLEM if maps dont fit on the memory!! *************** // ********* HERE: PROBLEM if maps dont fit on the memory!! ***************
// copying mapfile to ref map class // copying mapfile to ref map class
RefMap.filemap = mapfile.c_str(); RefMap.filemap = mapfile.c_str();
RefMap.readRefMaps(param); RefMap.readRefMaps(param);
//****************** Precalculating Necessary Stuff ********************* // ****************** Precalculating Necessary Stuff *********************
precalculate(); precalculate();
if (getenv("BIOEM_DEBUG_BREAK")) if (getenv("BIOEM_DEBUG_BREAK"))
...@@ -196,9 +196,9 @@ int bioem::configure(int ac, char* av[]) ...@@ -196,9 +196,9 @@ int bioem::configure(int ac, char* av[])
int bioem::precalculate() int bioem::precalculate()
{ {
//************************************************************************************** // **************************************************************************************
//* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels * //* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels *
//************************************************************************************** // **************************************************************************************
// Generating Grids of orientations // Generating Grids of orientations
param.CalculateGridsParam(); param.CalculateGridsParam();
...@@ -226,12 +226,12 @@ int bioem::precalculate() ...@@ -226,12 +226,12 @@ int bioem::precalculate()
int bioem::run() int bioem::run()
{ {
//************************************************************************************** // **************************************************************************************
//**** Main BioEM routine, projects, convolutes and compares with Map using OpenMP **** // **** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****
//************************************************************************************** // **************************************************************************************
//**** If we want to control the number of threads -> omp_set_num_threads(XX); ****** // **** If we want to control the number of threads -> omp_set_num_threads(XX); ******
//****************** Declarying class of Probability Pointer ************************* // ****************** Declarying class of Probability Pointer *************************
pProb = new bioem_Probability[RefMap.ntotRefMap]; pProb = new bioem_Probability[RefMap.ntotRefMap];
printf("\tInitializing\n"); printf("\tInitializing\n");
...@@ -247,12 +247,12 @@ int bioem::run() ...@@ -247,12 +247,12 @@ int bioem::run()
pProb[iRefMap].ConstAngle[iOrient] = -99999999; pProb[iRefMap].ConstAngle[iOrient] = -99999999;
} }
} }
//************************************************************************************** // **************************************************************************************
deviceStartRun(); deviceStartRun();
//******************************** MAIN CYCLE ****************************************** // ******************************** MAIN CYCLE ******************************************
//*** Declaring Private variables for each thread ***** // *** Declaring Private variables for each thread *****
mycomplex_t* proj_mapFFT; mycomplex_t* proj_mapFFT;
bioem_map conv_map; bioem_map conv_map;
mycomplex_t* conv_mapFFT; mycomplex_t* conv_mapFFT;
...@@ -268,25 +268,25 @@ int bioem::run() ...@@ -268,25 +268,25 @@ int bioem::run()
printf("\tInner Loop Count (%d %d %d) %lld\n", param.param_device.maxDisplaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, (long long int) (param.param_device.NumberPixels * param.param_device.NumberPixels * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1) * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1))); printf("\tInner Loop Count (%d %d %d) %lld\n", param.param_device.maxDisplaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, (long long int) (param.param_device.NumberPixels * param.param_device.NumberPixels * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1) * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1)));
for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++) for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
{ {
//*************************************************************************************** // ***************************************************************************************
//***** Creating Projection for given orientation and transforming to Fourier space ***** // ***** Creating Projection for given orientation and transforming to Fourier space *****
timer.ResetStart(); timer.ResetStart();
createProjection(iProjectionOut, proj_mapFFT); createProjection(iProjectionOut, proj_mapFFT);
printf("Time Projection %d: %f\n", iProjectionOut, timer.GetCurrentElapsedTime()); printf("Time Projection %d: %f\n", iProjectionOut, timer.GetCurrentElapsedTime());
//*************************************************************************************** // ***************************************************************************************
//***** **** Internal Loop over convolutions **** ***** // ***** **** Internal Loop over convolutions **** *****
for (int iConv = 0; iConv < param.nTotCTFs; iConv++) for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
{ {
printf("\t\tConvolution %d %d\n", iProjectionOut, iConv); printf("\t\tConvolution %d %d\n", iProjectionOut, iConv);
//*** Calculating convolutions of projection map and crosscorrelations *** // *** Calculating convolutions of projection map and crosscorrelations ***
timer.ResetStart(); timer.ResetStart();
createConvolutedProjectionMap(iProjectionOut, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV); createConvolutedProjectionMap(iProjectionOut, iConv, proj_mapFFT, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
printf("Time Convolution %d %d: %f\n", iProjectionOut, iConv, timer.GetCurrentElapsedTime()); printf("Time Convolution %d %d: %f\n", iProjectionOut, iConv, timer.GetCurrentElapsedTime());
//*************************************************************************************** // ***************************************************************************************
//*** Comparing each calculated convoluted map with all experimental maps *** // *** Comparing each calculated convoluted map with all experimental maps ***
timer.ResetStart(); timer.ResetStart();
compareRefMaps(iProjectionOut, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV); compareRefMaps(iProjectionOut, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
...@@ -307,9 +307,9 @@ int bioem::run() ...@@ -307,9 +307,9 @@ int bioem::run()
deviceFinishRun(); deviceFinishRun();
//************* Writing Out Probabilities *************** // ************* Writing Out Probabilities ***************
//*** Angular Probability *** // *** Angular Probability ***
// if(param.writeAngles){ // if(param.writeAngles){
ofstream angProbfile; ofstream angProbfile;
...@@ -321,12 +321,12 @@ int bioem::run() ...@@ -321,12 +321,12 @@ int bioem::run()
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++) for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{ {
//**** Total Probability *** // **** Total Probability ***
outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProb[iRefMap].Total) + pProb[iRefMap].Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd << "\n"; outputProbFile << "RefMap " << iRefMap << " Probability " << log(pProb[iRefMap].Total) + pProb[iRefMap].Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd << "\n";
outputProbFile << "RefMap " << iRefMap << " Maximizing Param: "; outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
//*** Param that maximize probability**** // *** Param that maximize probability****
outputProbFile << (pProb[iRefMap].max_prob + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " "; outputProbFile << (pProb[iRefMap].max_prob + 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[pProb[iRefMap].max_prob_orient].pos[0] << " "; outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[0] << " ";
outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[1] << " "; outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[1] << " ";
...@@ -338,7 +338,7 @@ int bioem::run() ...@@ -338,7 +338,7 @@ int bioem::run()
outputProbFile << pProb[iRefMap].max_prob_cent_y; outputProbFile << pProb[iRefMap].max_prob_cent_y;
outputProbFile << "\n"; outputProbFile << "\n";
//*** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap); // *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap);
if(param.writeAngles) if(param.writeAngles)
{ {
...@@ -378,8 +378,8 @@ int bioem::run() ...@@ -378,8 +378,8 @@ int bioem::run()
int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap) int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
{ {
//*************************************************************************************** // ***************************************************************************************
//***** BioEM routine for comparing reference maps to convoluted maps ***** // ***** BioEM routine for comparing reference maps to convoluted maps *****
if (FFTAlgo) //IF using the fft algorithm if (FFTAlgo) //IF using the fft algorithm
{ {
...@@ -418,8 +418,8 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m ...@@ -418,8 +418,8 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC) inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
{ {
//*************************************************************************************** // ***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm ***** // ***** Calculating cross correlation in FFTALGOrithm *****
const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.RefMapSize]; const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.RefMapSize];
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++) for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
...@@ -435,10 +435,10 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t ...@@ -435,10 +435,10 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t
int bioem::createProjection(int iMap, mycomplex_t* mapFFT) int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
{ {
//************************************************************************************** // **************************************************************************************
//**** BioEM Create Projection routine in Euler angle predefined grid**************** // **** BioEM Create Projection routine in Euler angle predefined grid****************
//********************* and turns projection into Fourier space ********************** // ********************* and turns projection into Fourier space **********************
//************************************************************************************** // **************************************************************************************
myfloat3_t RotatedPointsModel[Model.nPointsModel]; myfloat3_t RotatedPointsModel[Model.nPointsModel];
myfloat_t rotmat[3][3]; myfloat_t rotmat[3][3];
...@@ -452,9 +452,9 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -452,9 +452,9 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
beta = param.angles[iMap].pos[1]; beta = param.angles[iMap].pos[1];
gam = param.angles[iMap].pos[2]; gam = param.angles[iMap].pos[2];
//**** To see how things are going: cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; *** // **** To see how things are going: cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***
//********** Creat Rotation with pre-defiend grid of orientations********** // ********** Creat Rotation with pre-defiend grid of orientations**********
rotmat[0][0] = cos(gam) * cos(alpha) - cos(beta) * sin(alpha) * sin(gam); rotmat[0][0] = cos(gam) * cos(alpha) - cos(beta) * sin(alpha) * sin(gam);
rotmat[0][1] = cos(gam) * sin(alpha) + cos(beta) * cos(alpha) * sin(gam); rotmat[0][1] = cos(gam) * sin(alpha) + cos(beta) * cos(alpha) * sin(gam);
...@@ -485,7 +485,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -485,7 +485,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
int i, j; int i, j;
//************ Projection over the Z axis******************** // ************ Projection over the Z axis********************
for(int n = 0; n < Model.nPointsModel; n++) for(int n = 0; n < Model.nPointsModel; n++)
{ {
//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 ) //Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
...@@ -495,7 +495,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -495,7 +495,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
localproj[i * param.param_device.NumberPixels + j] += Model.densityPointsModel[n] / Model.NormDen; localproj[i * param.param_device.NumberPixels + j] += Model.densityPointsModel[n] / Model.NormDen;
} }
//**** Output Just to check**** // **** Output Just to check****
if(iMap == 10) if(iMap == 10)
{ {
ofstream myexamplemap; ofstream myexamplemap;
...@@ -513,8 +513,8 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -513,8 +513,8 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
myexampleRot.close(); myexampleRot.close();
} }
//***** Converting projection to Fourier Space for Convolution later with kernel**** // ***** Converting projection to Fourier Space for Convolution later with kernel****
//********** Omp Critical is necessary with FFTW******* // ********** Omp Critical is necessary with FFTW*******
myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT); myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT);
return(0); return(0);
...@@ -522,18 +522,18 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -522,18 +522,18 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC) int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
{ {
//************************************************************************************** // **************************************************************************************
//**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier ********** // **** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
//**************** calculated Projection with convoluted precalculated Kernel********** // **************** calculated Projection with convoluted precalculated Kernel**********
//*************** and Backtransforming it to real Space ****************************** // *************** and Backtransforming it to real Space ******************************
//************************************************************************************** // **************************************************************************************
myfloat_t* localconvFFT; myfloat_t* localconvFFT;
localconvFFT = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels); localconvFFT = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
mycomplex_t* tmp; mycomplex_t* tmp;
tmp = (mycomplex_t*) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D); tmp = (mycomplex_t*) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
//**** Multiplying FFTmap with corresponding kernel **** // **** Multiplying FFTmap with corresponding kernel ****
const mycomplex_t* refCTF = &param.refCTF[iConv * param.RefMapSize]; const mycomplex_t* refCTF = &param.refCTF[iConv * param.RefMapSize];
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++) for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
...@@ -546,10 +546,10 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -546,10 +546,10 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
//FFTW_C2R will destroy the input array, so we have to work on a copy here //FFTW_C2R will destroy the input array, so we have to work on a copy here
memcpy(tmp, localmultFFT, sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D); memcpy(tmp, localmultFFT, sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
//**** Bringing convoluted Map to real Space **** // **** Bringing convoluted Map to real Space ****
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, localconvFFT); myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, localconvFFT);
//****Asigning convolution fftw_complex to bioem_map **** // ****Asigning convolution fftw_complex to bioem_map ****
for(int i = 0; i < param.param_device.NumberPixels ; i++ ) for(int i = 0; i < param.param_device.NumberPixels ; i++ )
{ {
for(int j = 0; j < param.param_device.NumberPixels ; j++ ) for(int j = 0; j < param.param_device.NumberPixels ; j++ )
...@@ -558,7 +558,7 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -558,7 +558,7 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
} }
} }
//*** Calculating Cross-correlations of cal-convoluted map with its self ***** // *** Calculating Cross-correlations of cal-convoluted map with its self *****
sumC = 0; sumC = 0;
sumsquareC = 0; sumsquareC = 0;
for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++) for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
...@@ -566,14 +566,14 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -566,14 +566,14 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
sumC += localconvFFT[i]; sumC += localconvFFT[i];
sumsquareC += localconvFFT[i] * localconvFFT[i]; sumsquareC += localconvFFT[i] * localconvFFT[i];
} }
//*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier *** // *** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***
// Normalizing // Normalizing
myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels); myfloat_t norm2 = (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels);
myfloat_t norm4 = norm2 * norm2; myfloat_t norm4 = norm2 * norm2;
sumC = sumC / norm2; sumC = sumC / norm2;
sumsquareC = sumsquareC / norm4; sumsquareC = sumsquareC / norm4;
//**** Freeing fftw_complex created (dont know if omp critical is necessary) **** // **** Freeing fftw_complex created (dont know if omp critical is necessary) ****
myfftw_free(localconvFFT); myfftw_free(localconvFFT);
myfftw_free(tmp); myfftw_free(tmp);
...@@ -582,7 +582,7 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -582,7 +582,7 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
int bioem::calcross_cor(bioem_map& localmap, myfloat_t& sum, myfloat_t& sumsquare) int bioem::calcross_cor(bioem_map& localmap, myfloat_t& sum, myfloat_t& sumsquare)
{ {
//*********************** Routine to calculate Cross correlations*********************** // *********************** Routine to calculate Cross correlations***********************
sum = 0.0; sum = 0.0;
sumsquare = 0.0; sumsquare = 0.0;
......
...@@ -14,8 +14,8 @@ ...@@ -14,8 +14,8 @@
template <int GPUAlgo> template <int GPUAlgo>
__device__ static inline void update_prob(const myfloat_t logpro, const int iRefMap, const int iOrient, const int iConv, const int cent_x, const int cent_y, bioem_Probability* pProb, myfloat_t* buf3 = NULL, int* bufint = NULL) __device__ static inline void update_prob(const myfloat_t logpro, const int iRefMap, const int iOrient, const int iConv, const int cent_x, const int cent_y, bioem_Probability* pProb, myfloat_t* buf3 = NULL, int* bufint = NULL)
{ {
/******* Summing total Probabilities *************/ // ******* Summing total Probabilities *************
/******* Need a constant because of numerical divergence*****/ // ******* Need a constant because of numerical divergence*****
if(pProb[iRefMap].Constoadd < logpro) if(pProb[iRefMap].Constoadd < logpro)
{ {
pProb[iRefMap].Total = pProb[iRefMap].Total * exp(-logpro + pProb[iRefMap].Constoadd); pProb[iRefMap].Total = pProb[iRefMap].Total * exp(-logpro + pProb[iRefMap].Constoadd);
...@@ -35,7 +35,7 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef ...@@ -35,7 +35,7 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef
pProb[iRefMap].forAngles[iOrient] += exp(logpro - pProb[iRefMap].ConstAngle[iOrient]); pProb[iRefMap].forAngles[iOrient] += exp(logpro - pProb[iRefMap].ConstAngle[iOrient]);
} }
/********** Getting parameters that maximize the probability ***********/ // ********** Getting parameters that maximize the probability ***********
if(pProb[iRefMap].max_prob < logpro) if(pProb[iRefMap].max_prob < logpro)
{ {
pProb[iRefMap].max_prob = logpro; pProb[iRefMap].max_prob = logpro;
...@@ -64,7 +64,7 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param, ...@@ -64,7 +64,7 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
const myfloat_t firstele = param.Ntotpi * (sumsquareref * sumsquare - crossproMapConv * crossproMapConv) + const myfloat_t firstele = param.Ntotpi * (sumsquareref * sumsquare - crossproMapConv * crossproMapConv) +
2 * sumref * sum * crossproMapConv - sumsquareref * sum * sum - sumref * sumref * sumsquare; 2 * sumref * sum * crossproMapConv - sumsquareref * sum * sum - sumref * sumref * sumsquare;
//******* Calculating log of Prob*********/ // ******* Calculating log of Prob*********
// As in fortran code: logpro=(3-Ntotpi)*0.5*log(firstele/pConvMap[iOrient].ForLogProbfromConv[iConv])+(Ntotpi*0.5-2)*log(Ntotpi-2)-0.5*log(pConvMap[iOrient].ForLogProbfromConv[iConv])+0.5*log(PI)+(1-Ntotpi*0.5)*(log(2*PI)+1); // As in fortran code: logpro=(3-Ntotpi)*0.5*log(firstele/pConvMap[iOrient].ForLogProbfromConv[iConv])+(Ntotpi*0.5-2)*log(Ntotpi-2)-0.5*log(pConvMap[iOrient].ForLogProbfromConv[iConv])+0.5*log(PI)+(1-Ntotpi*0.5)*(log(2*PI)+1);
const myfloat_t logpro = (3 - param.Ntotpi) * 0.5 * log(firstele) + (param.Ntotpi * 0.5 - 2) * log((param.Ntotpi - 2) * ForLogProb); const myfloat_t logpro = (3 - param.Ntotpi) * 0.5 * log(firstele) + (param.Ntotpi * 0.5 - 2) * log((param.Ntotpi - 2) * ForLogProb);
return(logpro); return(logpro);
...@@ -72,9 +72,9 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param, ...@@ -72,9 +72,9 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
__device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, float value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap) __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, float value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
{ {
/********************************************************/ // ********************************************************
/*********** Calculates the BioEM probability ***********/ // *********** Calculates the BioEM probability ***********
/********************************************************/ // ********************************************************
const myfloat_t logpro = calc_logpro(param, sumC, sumsquareC, value, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]); const myfloat_t logpro = calc_logpro(param, sumC, sumsquareC, value, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]);
...@@ -134,13 +134,13 @@ template <int GPUAlgo, class RefT> ...@@ -134,13 +134,13 @@ template <int GPUAlgo, class RefT>
__device__ static inline void compareRefMap(const int iRefMap, const int iOrient, const int iConv, const bioem_map& Mapconv, bioem_Probability* pProb, const bioem_param_device& param, const RefT& RefMap, __device__ static inline void compareRefMap(const int iRefMap, const int iOrient, const int iConv, const bioem_map& Mapconv, bioem_Probability* pProb, const bioem_param_device& param, const RefT& RefMap,
const int cent_x, const int cent_y, const int myShift = 0, const int nShifts2 = 0, const int myRef = 0, const bool threadActive = true) const int cent_x, const int cent_y, const int myShift = 0, const int nShifts2 = 0, const int myRef = 0, const bool threadActive = true)
{ {
/**************************************************************************************/ // **************************************************************************************
/********************** Calculating BioEM Probability ********************************/ // ********************** Calculating BioEM Probability ********************************
/************************* Loop of center displacement here ***************************/ // ************************* Loop of center displacement here ***************************
// Taking into account the center displacement // Taking into account the center displacement
/*** Inizialzing crosscorrelations of calculated projected convolutions ***/ // *** Inizialzing crosscorrelations of calculated projected convolutions ***
#ifdef SSECODE #ifdef SSECODE
myfloat_t sum, sumsquare, crossproMapConv; myfloat_t sum, sumsquare, crossproMapConv;
__m128 sum_v = _mm_setzero_ps(), sumsquare_v = _mm_setzero_ps(), cross_v = _mm_setzero_ps(), d1, d2; __m128 sum_v = _mm_setzero_ps(), sumsquare_v = _mm_setzero_ps(), cross_v = _mm_setzero_ps(), d1, d2;
...@@ -149,7 +149,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient ...@@ -149,7 +149,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
myfloat_t sumsquare = 0.0; myfloat_t sumsquare = 0.0;
myfloat_t crossproMapConv = 0.0; myfloat_t crossproMapConv = 0.0;
#endif #endif
/****** Loop over Pixels to calculate dot product and cross-correlations of displaced Ref Conv. Map***/ // ****** Loop over Pixels to calculate dot product and cross-correlations of displaced Ref Conv. Map***
myfloat_t logpro; myfloat_t logpro;
if (GPUAlgo != 2 || threadActive) if (GPUAlgo != 2 || threadActive)
{ {
...@@ -217,7 +217,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient ...@@ -217,7 +217,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
crossproMapConv = _mm_cvtss_f32(cross_v); crossproMapConv = _mm_cvtss_f32(cross_v);
#endif #endif
/********** Calculating elements in BioEM Probability formula ********/ // ********** Calculating elements in BioEM Probability formula ********
logpro = calc_logpro(param, sum, sumsquare, crossproMapConv, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]); logpro = calc_logpro(param, sum, sumsquare, crossproMapConv, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]);
} }
else else
...@@ -380,7 +380,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient ...@@ -380,7 +380,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
else else
#endif #endif
/***** Summing & Storing total/Orientation Probabilites for each map ************/ // ***** Summing & Storing total/Orientation Probabilites for each map ************
{ {
update_prob < -1 > (logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb); update_prob < -1 > (logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb);
} }
......
...@@ -21,9 +21,9 @@ ...@@ -21,9 +21,9 @@
int main(int argc, char* argv[]) int main(int argc, char* argv[])
{ {
//************************************************************************************* // *************************************************************************************
//********************************* Main BioEM code ********************************** // ********************************* Main BioEM code **********************************
//************************************************************************************* // *************************************************************************************
#ifdef _MM_SET_DENORMALS_ZERO_MODE #ifdef _MM_SET_DENORMALS_ZERO_MODE
#pragma omp parallel #pragma omp parallel
...@@ -45,17 +45,17 @@ int main(int argc, char* argv[]) ...@@ -45,17 +45,17 @@ int main(int argc, char* argv[])
bio = new bioem; bio = new bioem;
} }
//************ Configuration and Pre-calculating necessary objects ********************* // ************ Configuration and Pre-calculating necessary objects *********************
printf("Configuring\n"); printf("Configuring\n");
bio->configure(argc, argv); bio->configure(argc, argv);
//******************************* Run BioEM routine ************************************