Commit 860031e3 authored by David Rohr's avatar David Rohr

create angular probabilites only if writeAngles set

parent 88c93497
...@@ -92,7 +92,7 @@ int bioem::configure(int ac, char* av[]) ...@@ -92,7 +92,7 @@ int bioem::configure(int ac, char* av[])
// *** 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.param_device.writeAngles = false;
param.dumpMap = false; param.dumpMap = false;
param.loadMap = false; param.loadMap = false;
RefMap.readMRC = false; RefMap.readMRC = false;
...@@ -280,12 +280,15 @@ int bioem::run() ...@@ -280,12 +280,15 @@ int bioem::run()
pProbMap.Total = 0.0; pProbMap.Total = 0.0;
pProbMap.Constoadd = -9999999; pProbMap.Constoadd = -9999999;
pProbMap.max_prob = -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.forAngles = 0.0;
pProbAngle.ConstAngle = -99999999; pProbAngle.ConstAngle = -99999999;
}
} }
} }
// ************************************************************************************** // **************************************************************************************
...@@ -377,7 +380,7 @@ int bioem::run() ...@@ -377,7 +380,7 @@ int bioem::run()
// *** Angular Probability *** // *** Angular Probability ***
ofstream angProbfile; ofstream angProbfile;
if(param.writeAngles) if(param.param_device.writeAngles)
{ {
angProbfile.open ("ANG_PROB"); angProbfile.open ("ANG_PROB");
} }
...@@ -407,7 +410,7 @@ int bioem::run() ...@@ -407,7 +410,7 @@ int bioem::run()
// *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap); // *** For individual files*** //angProbfile.open ("ANG_PROB_"iRefMap);
if(param.writeAngles) if(param.param_device.writeAngles)
{ {
for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++) for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient++)
{ {
...@@ -418,7 +421,7 @@ int bioem::run() ...@@ -418,7 +421,7 @@ int bioem::run()
} }
} }
if(param.writeAngles) if(param.param_device.writeAngles)
{ {
angProbfile.close(); angProbfile.close();
} }
...@@ -540,6 +543,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -540,6 +543,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
} }
// **** Output Just to check**** // **** Output Just to check****
#ifdef PILAR_DEBUG
if(iMap == 10) if(iMap == 10)
{ {
ofstream myexamplemap; ofstream myexamplemap;
...@@ -556,6 +560,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -556,6 +560,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
myexamplemap.close(); myexamplemap.close();
myexampleRot.close(); myexampleRot.close();
} }
#endif
// ***** 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*******
......
...@@ -21,30 +21,30 @@ ...@@ -21,30 +21,30 @@
#endif #endif
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, bool doAngle, 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*****
bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap); bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
if(pProbMap.Constoadd < logpro) if(pProbMap.Constoadd < logpro)
{ {
pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd); pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd);
pProbMap.Constoadd = logpro; pProbMap.Constoadd = logpro;
} }
if (GPUAlgo != 2) pProbMap.Total += exp(logpro - pProbMap.Constoadd);
//Summing probabilities for each orientation if (doAngle)
if(pProbAngle.ConstAngle < logpro)
{ {
pProbAngle.forAngles = pProbAngle.forAngles * exp(-logpro + pProbAngle.ConstAngle); bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
pProbAngle.ConstAngle = logpro;
}
if (GPUAlgo != 2) //Summing probabilities for each orientation
{ if(pProbAngle.ConstAngle < logpro)
pProbMap.Total += exp(logpro - pProbMap.Constoadd); {
pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle); 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 *********** // ********** Getting parameters that maximize the probability ***********
...@@ -90,11 +90,9 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo ...@@ -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]); 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 //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) if(pProbMap.Constoadd < logpro)
{ {
pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd); pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd);
...@@ -102,13 +100,18 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo ...@@ -102,13 +100,18 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
} }
pProbMap.Total += exp(logpro - pProbMap.Constoadd); pProbMap.Total += exp(logpro - pProbMap.Constoadd);
if(pProbAngle.ConstAngle < logpro) if (param.writeAngles)
{ {
pProbAngle.forAngles = pProbAngle.forAngles * exp(-logpro + pProbAngle.ConstAngle); bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
pProbAngle.ConstAngle = logpro;
}
pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle);
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) if(pProbMap.max_prob < logpro)
{ {
pProbMap.max_prob = logpro; pProbMap.max_prob = logpro;
...@@ -291,7 +294,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient ...@@ -291,7 +294,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
if (myShift == 0 && iRefMap < RefMap.ntotRefMap) if (myShift == 0 && iRefMap < RefMap.ntotRefMap)
{ {
const myfloat_t logpro_max = vbuf[myThreadIdxX]; const myfloat_t logpro_max = vbuf[myThreadIdxX];
update_prob<GPUAlgo>(logpro_max, iRefMap, iOrient, iConv, -1, -1, pProb, buf3, bufint); update_prob<GPUAlgo>(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 ...@@ -399,7 +402,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
else else
#endif #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);
} }
} }
......
...@@ -374,7 +374,7 @@ int bioem_cuda::deviceInit() ...@@ -374,7 +374,7 @@ int bioem_cuda::deviceInit()
gpumap->sum_RefMap = sum; gpumap->sum_RefMap = sum;
gpumap->sumsquare_RefMap = sumsquare; 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++) for (int i = 0; i < 2; i++)
{ {
checkCudaErrors(cudaStreamCreate(&cudaStream[i])); checkCudaErrors(cudaStreamCreate(&cudaStream[i]));
...@@ -466,7 +466,7 @@ int bioem_cuda::deviceStartRun() ...@@ -466,7 +466,7 @@ int bioem_cuda::deviceStartRun()
pProb_device = *pProb_host; pProb_device = *pProb_host;
pProb_device.ptr = pProb_memory; pProb_device.ptr = pProb_memory;
pProb_device.set_pointers(); 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) if (FFTAlgo)
{ {
...@@ -501,7 +501,7 @@ int bioem_cuda::deviceStartRun() ...@@ -501,7 +501,7 @@ int bioem_cuda::deviceStartRun()
int bioem_cuda::deviceFinishRun() int bioem_cuda::deviceFinishRun()
{ {
if (GPUAsync) cudaStreamSynchronize(cudaStream[0]); 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) if (FFTAlgo)
{ {
......
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#define BIOEM_MAP_H #define BIOEM_MAP_H
#include "defs.h" #include "defs.h"
#include "param.h"
#include <complex> #include <complex>
#include <math.h> #include <math.h>
...@@ -107,9 +108,11 @@ public: ...@@ -107,9 +108,11 @@ public:
bioem_Probability_map* ptr_map; bioem_Probability_map* ptr_map;
bioem_Probability_angle* ptr_angle; 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); void init(size_t maps, size_t angles, bioem& bio);
......
...@@ -18,6 +18,8 @@ public: ...@@ -18,6 +18,8 @@ public:
int NtotDist; int NtotDist;
myfloat_t Ntotpi; myfloat_t Ntotpi;
myfloat_t volu; myfloat_t volu;
// If to write Probabilities of Angles from Model
bool writeAngles;
}; };
class bioem_param class bioem_param
...@@ -38,8 +40,6 @@ public: ...@@ -38,8 +40,6 @@ public:
mycomplex_t* refCTF; mycomplex_t* refCTF;
myfloat3_t* CtfParam; myfloat3_t* CtfParam;
// If to write Probabilities of Angles from Model
bool writeAngles;
myfloat_t pixelSize; myfloat_t pixelSize;
// Grid Points in Euler angles, assuming uniform sampling d_alpha=d_gamma (in 2pi) & cos(beta)=-1,1 // Grid Points in Euler angles, assuming uniform sampling d_alpha=d_gamma (in 2pi) & cos(beta)=-1,1
int angleGridPointsAlpha; int angleGridPointsAlpha;
......
...@@ -276,7 +276,7 @@ void bioem_Probability::init(size_t maps, size_t angles, bioem& bio) ...@@ -276,7 +276,7 @@ void bioem_Probability::init(size_t maps, size_t angles, bioem& bio)
{ {
nMaps = maps; nMaps = maps;
nAngles = angles; 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(); set_pointers();
} }
...@@ -286,11 +286,14 @@ void bioem_Probability::copyFrom(bioem_Probability* from, bioem& bio) ...@@ -286,11 +286,14 @@ void bioem_Probability::copyFrom(bioem_Probability* from, bioem& bio)
bioem_Probability_map& pProbMapFrom = from->getProbMap(0); bioem_Probability_map& pProbMapFrom = from->getProbMap(0);
memcpy(&pProbMap, &pProbMapFrom, from->nMaps * sizeof(bioem_Probability_map)); 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); for (int iOrient = 0; iOrient < bio.param.nTotGridAngles; iOrient ++)
bioem_Probability_angle& pProbAngleFrom = from->getProbAngle(0, iOrient); {
memcpy(&pProbAngle, &pProbAngleFrom, from->nMaps * sizeof(bioem_Probability_angle)); bioem_Probability_angle& pProbAngle = getProbAngle(0, iOrient);
bioem_Probability_angle& pProbAngleFrom = from->getProbAngle(0, iOrient);
memcpy(&pProbAngle, &pProbAngleFrom, from->nMaps * sizeof(bioem_Probability_angle));
}
} }
} }
......
...@@ -209,7 +209,7 @@ int bioem_param::readParameters(const char* fileinput) ...@@ -209,7 +209,7 @@ int bioem_param::readParameters(const char* fileinput)
} }
else if (strcmp(token, "WRITE_PROB_ANGLES") == 0) else if (strcmp(token, "WRITE_PROB_ANGLES") == 0)
{ {
writeAngles = true; param_device.writeAngles = true;
cout << "Writing Probabilies of each angle \n"; cout << "Writing Probabilies of each angle \n";
} }
} }
......
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