diff --git a/bioem.cpp b/bioem.cpp index c07e4e09ffae304fd3eb24bf9d281822ddbff965..04d53b5619e1c8647f26440ab0c7f64b0e48ddcf 100644 --- a/bioem.cpp +++ b/bioem.cpp @@ -92,7 +92,7 @@ int bioem::configure(int ac, char* av[]) // *** Inizialzing default variables *** std::string infile, modelfile, mapfile; Model.readPDB = false; - param.writeAngles = false; + param.param_device.writeAngles = false; param.dumpMap = false; param.loadMap = false; RefMap.readMRC = false; @@ -280,12 +280,15 @@ int bioem::run() pProbMap.Total = 0.0; pProbMap.Constoadd = -9999999; pProbMap.max_prob = -9999999; - for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++) + if (param.param_device.writeAngles) { - bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); + for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++) + { + bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); - pProbAngle.forAngles = 0.0; - pProbAngle.ConstAngle = -99999999; + pProbAngle.forAngles = 0.0; + pProbAngle.ConstAngle = -99999999; + } } } // ************************************************************************************** @@ -377,7 +380,7 @@ int bioem::run() // *** Angular Probability *** ofstream angProbfile; - if(param.writeAngles) + if(param.param_device.writeAngles) { angProbfile.open ("ANG_PROB"); } @@ -407,7 +410,7 @@ int bioem::run() // *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap); - if(param.writeAngles) + if(param.param_device.writeAngles) { for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++) { @@ -418,7 +421,7 @@ int bioem::run() } } - if(param.writeAngles) + if(param.param_device.writeAngles) { angProbfile.close(); } @@ -540,6 +543,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) } // **** Output Just to check**** +#ifdef PILAR_DEBUG if(iMap == 10) { ofstream myexamplemap; @@ -556,6 +560,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) myexamplemap.close(); myexampleRot.close(); } +#endif // ***** Converting projection to Fourier Space for Convolution later with kernel**** // ********** Omp Critical is necessary with FFTW******* diff --git a/bioem_algorithm.h b/bioem_algorithm.h index edab147dfdbfada44289920d8dccb9c91e907032..d8ab2a48c588e1cd6b8ad647c0d93c1798085933 100644 --- a/bioem_algorithm.h +++ b/bioem_algorithm.h @@ -21,30 +21,30 @@ #endif template -__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, bool doAngle, myfloat_t* buf3 = NULL, int* bufint = NULL) { // ******* Summing total Probabilities ************* // ******* Need a constant because of numerical divergence***** bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap); - bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); - if(pProbMap.Constoadd < logpro) { pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd); pProbMap.Constoadd = logpro; } + if (GPUAlgo != 2) pProbMap.Total += exp(logpro - pProbMap.Constoadd); - //Summing probabilities for each orientation - if(pProbAngle.ConstAngle < logpro) + if (doAngle) { - pProbAngle.forAngles = pProbAngle.forAngles * exp(-logpro + pProbAngle.ConstAngle); - pProbAngle.ConstAngle = logpro; - } + bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); - if (GPUAlgo != 2) - { - pProbMap.Total += exp(logpro - pProbMap.Constoadd); - pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle); + //Summing probabilities for each orientation + if(pProbAngle.ConstAngle < logpro) + { + pProbAngle.forAngles = pProbAngle.forAngles * exp(-logpro + pProbAngle.ConstAngle); + pProbAngle.ConstAngle = logpro; + } + + if (GPUAlgo != 2) pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle); } // ********** Getting parameters that maximize the probability *********** @@ -90,11 +90,9 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo const myfloat_t logpro = calc_logpro(param, sumC, sumsquareC, value, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]); - bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap); - bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); - //GCC is too stupid to inline properly, so the code is copied here - //update_prob<-1>(logpro, iRefMap, iOrient, iConv, disx, disy, pProb); + //update_prob<-1>(logpro, iRefMap, iOrient, iConv, disx, disy, pProb, param.writeAngles); + bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap); if(pProbMap.Constoadd < logpro) { pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd); @@ -102,13 +100,18 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo } pProbMap.Total += exp(logpro - pProbMap.Constoadd); - if(pProbAngle.ConstAngle < logpro) + if (param.writeAngles) { - pProbAngle.forAngles = pProbAngle.forAngles * exp(-logpro + pProbAngle.ConstAngle); - pProbAngle.ConstAngle = logpro; - } - pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle); + bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient); + if(pProbAngle.ConstAngle < logpro) + { + pProbAngle.forAngles = pProbAngle.forAngles * exp(-logpro + pProbAngle.ConstAngle); + pProbAngle.ConstAngle = logpro; + } + pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle); + } + if(pProbMap.max_prob < logpro) { pProbMap.max_prob = logpro; @@ -291,7 +294,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient if (myShift == 0 && iRefMap < RefMap.ntotRefMap) { const myfloat_t logpro_max = vbuf[myThreadIdxX]; - update_prob(logpro_max, iRefMap, iOrient, iConv, -1, -1, pProb, buf3, bufint); + update_prob(logpro_max, iRefMap, iOrient, iConv, -1, -1, pProb, param.writeAngles, buf3, bufint); } } @@ -399,7 +402,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient else #endif { - update_prob < -1 > (logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb); + update_prob < -1 > (logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb, param.writeAngles); } } diff --git a/bioem_cuda.cu b/bioem_cuda.cu index 6bb247dd4582c191097af9741d389be36fa4e187..09e1db5504d7a0926b88d4ab24a5ecd405cdcf79 100644 --- a/bioem_cuda.cu +++ b/bioem_cuda.cu @@ -374,7 +374,7 @@ int bioem_cuda::deviceInit() gpumap->sum_RefMap = sum; gpumap->sumsquare_RefMap = sumsquare; - checkCudaErrors(cudaMalloc(&pProb_memory, pProb_device.get_size(RefMap.ntotRefMap, param.nTotGridAngles))); + checkCudaErrors(cudaMalloc(&pProb_memory, pProb_device.get_size(RefMap.ntotRefMap, param.nTotGridAngles, param.param_device.writeAngles))); for (int i = 0; i < 2; i++) { checkCudaErrors(cudaStreamCreate(&cudaStream[i])); @@ -466,7 +466,7 @@ int bioem_cuda::deviceStartRun() pProb_device = *pProb_host; pProb_device.ptr = pProb_memory; pProb_device.set_pointers(); - checkCudaErrors(cudaMemcpyAsync(pProb_device.ptr, pProb_host->ptr, pProb_host->get_size(maxRef, param.nTotGridAngles), cudaMemcpyHostToDevice, cudaStream[0])); + checkCudaErrors(cudaMemcpyAsync(pProb_device.ptr, pProb_host->ptr, pProb_host->get_size(maxRef, param.nTotGridAngles, param.param_device.writeAngles), cudaMemcpyHostToDevice, cudaStream[0])); if (FFTAlgo) { @@ -501,7 +501,7 @@ int bioem_cuda::deviceStartRun() int bioem_cuda::deviceFinishRun() { if (GPUAsync) cudaStreamSynchronize(cudaStream[0]); - checkCudaErrors(cudaMemcpyAsync(pProb_host->ptr, pProb_device.ptr, pProb_host->get_size(maxRef, param.nTotGridAngles), cudaMemcpyDeviceToHost, cudaStream[0])); + checkCudaErrors(cudaMemcpyAsync(pProb_host->ptr, pProb_device.ptr, pProb_host->get_size(maxRef, param.nTotGridAngles, param.param_device.writeAngles), cudaMemcpyDeviceToHost, cudaStream[0])); if (FFTAlgo) { diff --git a/include/map.h b/include/map.h index 89b4a31156151f9273fe0ecf547a33ebc345cb6e..1495ec4fde77c6aa99edbe9d0409feeec971b5bd 100644 --- a/include/map.h +++ b/include/map.h @@ -2,6 +2,7 @@ #define BIOEM_MAP_H #include "defs.h" +#include "param.h" #include #include @@ -107,9 +108,11 @@ public: bioem_Probability_map* ptr_map; bioem_Probability_angle* ptr_angle; - static size_t get_size(size_t maps, size_t angles) + static size_t get_size(size_t maps, size_t angles, bool writeAngles) { - return(maps * (angles * sizeof(bioem_Probability_angle) + sizeof(bioem_Probability_map))); + size_t size = sizeof(bioem_Probability_map); + if (writeAngles) size += angles * sizeof(bioem_Probability_angle); + return(maps * size); } void init(size_t maps, size_t angles, bioem& bio); diff --git a/include/param.h b/include/param.h index 2ed175dfe9bc11917f84203c9a5aee3063e619b3..1a471b538f1f8f0527b294d48b3b2c3b5313461a 100644 --- a/include/param.h +++ b/include/param.h @@ -18,6 +18,8 @@ public: int NtotDist; myfloat_t Ntotpi; myfloat_t volu; +// If to write Probabilities of Angles from Model + bool writeAngles; }; class bioem_param @@ -38,8 +40,6 @@ public: mycomplex_t* refCTF; myfloat3_t* CtfParam; -// If to write Probabilities of Angles from Model - bool writeAngles; myfloat_t pixelSize; // Grid Points in Euler angles, assuming uniform sampling d_alpha=d_gamma (in 2pi) & cos(beta)=-1,1 int angleGridPointsAlpha; diff --git a/map.cpp b/map.cpp index 020c765d982538daaee32f0f1e4bf98146d395a2..a65effb1cb9dd53acf95832e593cb652a677bd35 100644 --- a/map.cpp +++ b/map.cpp @@ -276,7 +276,7 @@ void bioem_Probability::init(size_t maps, size_t angles, bioem& bio) { nMaps = maps; nAngles = angles; - ptr = bio.malloc_device_host(get_size(maps, angles)); + ptr = bio.malloc_device_host(get_size(maps, angles, bio.param.param_device.writeAngles)); set_pointers(); } @@ -286,11 +286,14 @@ void bioem_Probability::copyFrom(bioem_Probability* from, bioem& bio) bioem_Probability_map& pProbMapFrom = from->getProbMap(0); memcpy(&pProbMap, &pProbMapFrom, from->nMaps * sizeof(bioem_Probability_map)); - for (int iOrient = 0; iOrient < bio.param.nTotGridAngles; iOrient ++) + if (bio.param.param_device.writeAngles) { - bioem_Probability_angle& pProbAngle = getProbAngle(0, iOrient); - bioem_Probability_angle& pProbAngleFrom = from->getProbAngle(0, iOrient); - memcpy(&pProbAngle, &pProbAngleFrom, from->nMaps * sizeof(bioem_Probability_angle)); + for (int iOrient = 0; iOrient < bio.param.nTotGridAngles; iOrient ++) + { + bioem_Probability_angle& pProbAngle = getProbAngle(0, iOrient); + bioem_Probability_angle& pProbAngleFrom = from->getProbAngle(0, iOrient); + memcpy(&pProbAngle, &pProbAngleFrom, from->nMaps * sizeof(bioem_Probability_angle)); + } } } diff --git a/param.cpp b/param.cpp index 43c26bf1714420fb17cfbe1e91b843282a3a722a..a41a5a1a52f9e86ac5a9e452750b499d1a8e2264 100644 --- a/param.cpp +++ b/param.cpp @@ -209,7 +209,7 @@ int bioem_param::readParameters(const char* fileinput) } else if (strcmp(token, "WRITE_PROB_ANGLES") == 0) { - writeAngles = true; + param_device.writeAngles = true; cout << "Writing Probabilies of each angle \n"; } }