Skip to content
GitLab
Commit 01a19a3d authored by Pilar Cossio-Tejada's avatar Pilar Cossio-Tejada
Browse files

Polishing & including priors without FFTALGO

parent 470535fe
Expand all Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ 01a19a3d
...@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.6) ...@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.6)
project(BioEM) project(BioEM)
###Set up options ###Set up options
option (INCLUDE_CUDA "Build BioEM with CUDA support" ON) option (INCLUDE_CUDA "Build BioEM with CUDA support" OFF)
option (INCLUDE_OPENMP "Build BioEM with OpenMP support" ON) option (INCLUDE_OPENMP "Build BioEM with OpenMP support" ON)
option (INCLUDE_MPI "Build BioEM with MPI support" ON) option (INCLUDE_MPI "Build BioEM with MPI support" ON)
option (PRINT_CMAKE_VARIABLES "List all CMAKE Variables" OFF) option (PRINT_CMAKE_VARIABLES "List all CMAKE Variables" OFF)
...@@ -23,7 +23,7 @@ else() ...@@ -23,7 +23,7 @@ else()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${BIOEM_GCC_FLAGS}") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${BIOEM_GCC_FLAGS}")
endif() endif()
set (BIOEM_SOURCE_FILES "bioem.cpp" "main.cpp" "map.cpp" "model.cpp" "param.cpp" "cmodules/timer.cpp") set (BIOEM_SOURCE_FILES "bioem.cpp" "main.cpp" "map.cpp" "model.cpp" "param.cpp" "timer.cpp")
###Find Required Packages ###Find Required Packages
find_package(PkgConfig) find_package(PkgConfig)
......
License.txt
View file @ 01a19a3d
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images> < BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2016 Pilar Cossio, David Rohr and Gerhard Hummer Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
Volker Lindenstruth and Gerhard Hummer.
This program is free software: you can redistribute it and/or modify This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
......
Manual/BioEM_Manual.tex
View file @ 01a19a3d
This diff is collapsed.
Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision
View file @ 01a19a3d
...@@ -6,14 +6,14 @@ PIXEL_SIZE 1.32 ...@@ -6,14 +6,14 @@ PIXEL_SIZE 1.32
USE_QUATERNIONS USE_QUATERNIONS
##### Constrast transfer integration: ####### ##### Constrast transfer integration: #######
CTF_B_ENV 10.0 500. 5 CTF_B_ENV 2.0 300. 4
CTF_DEFOCUS 0.5 4.5 8 CTF_DEFOCUS 0.5 4.5 8
CTF_AMPLITUDE 0.01 0.601 5 CTF_AMPLITUDE 0.01 0.601 5
SIGMA_PRIOR_B_CTF 50. SIGMA_PRIOR_B_CTF 50.
SIGMA_PRIOR_DEFOCUS 0.4 SIGMA_PRIOR_DEFOCUS 0.4
PRIOR_DEFOCUS_CENTER 2.7 PRIOR_DEFOCUS_CENTER 2.8
##### Center displacement: ####### ##### Center displacement: #######
......
Tutorial_BioEM/Param_Input_Priors
View file @ 01a19a3d
...@@ -17,3 +17,8 @@ DISPLACE_CENTER 10 2 ...@@ -17,3 +17,8 @@ DISPLACE_CENTER 10 2
####### Priors ############### ####### Priors ###############
PRIOR_MODEL 0.005 PRIOR_MODEL 0.005
PRIOR_ANGLES PRIOR_ANGLES
SIGMA_PRIOR_DEFOCUS 0.4
PRIOR_DEFOCUS_CENTER 1.8
SIGMA_PRIOR_B_CTF 10.
Tutorial_BioEM/Param_ProRun
View file @ 01a19a3d
...@@ -6,9 +6,9 @@ PIXEL_SIZE 1.77 ...@@ -6,9 +6,9 @@ PIXEL_SIZE 1.77
USE_QUATERNIONS USE_QUATERNIONS
##### Constrast transfer integration: ####### ##### Constrast transfer integration: #######
CTF_B_ENV 10.0 500. 5. CTF_B_ENV 2.0 300. 4.
CTF_DEFOCUS 0.5 4.5 8. CTF_DEFOCUS 0.5 4.5 8.
CTF_AMPLITUDE 0.01 0.401 4. CTF_AMPLITUDE 0.01 0.601 5.
## Gaussian width of Prior of b-parameter ## Gaussian width of Prior of b-parameter
SIGMA_PRIOR_B_CTF 50. SIGMA_PRIOR_B_CTF 50.
......
bioem.cpp
View file @ 01a19a3d
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images> < BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2016 Pilar Cossio, David Rohr and Gerhard Hummer. Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
Volker Lindenstruth and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany. Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution. See license statement for terms of distribution.
...@@ -35,7 +36,7 @@ ...@@ -35,7 +36,7 @@
#include <fftw3.h> #include <fftw3.h>
#include <math.h> #include <math.h>
#include "cmodules/timer.h" #include "timer.h"
#include "param.h" #include "param.h"
#include "bioem.h" #include "bioem.h"
...@@ -816,8 +817,8 @@ int bioem::run() ...@@ -816,8 +817,8 @@ int bioem::run()
else{outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " [1/A²] ";} else{outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << " [1/A²] ";}
outputProbFile << pProbMap.max.max_prob_cent_x << " [pix] "; outputProbFile << pProbMap.max.max_prob_cent_x << " [pix] ";
outputProbFile << pProbMap.max.max_prob_cent_y << " [pix] " ; outputProbFile << pProbMap.max.max_prob_cent_y << " [pix] " ;
outputProbFile << pProbMap.max.max_prob_norm << " [] " ; if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_norm << " [] " ;}else{outputProbFile << "N.A." << " [] ";}
outputProbFile << pProbMap.max.max_prob_mu << " [] "; if(FFTAlgo){outputProbFile << pProbMap.max.max_prob_mu << " [] ";}else{outputProbFile << "N.A." << " [] ";}
outputProbFile << "\n"; outputProbFile << "\n";
// Writing out CTF parameters if requiered // Writing out CTF parameters if requiered
...@@ -827,8 +828,6 @@ int bioem::run() ...@@ -827,8 +828,6 @@ int bioem::run()
denomi = param.CtfParam[pProbMap.max.max_prob_conv].pos[1] * param.CtfParam[pProbMap.max.max_prob_conv].pos[1] + 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]; param.CtfParam[pProbMap.max.max_prob_conv].pos[2] * param.CtfParam[pProbMap.max.max_prob_conv].pos[2];
outputProbFile << "RefMap: " << iRefMap << " CTFMaxParam: "; outputProbFile << "RefMap: " << iRefMap << " CTFMaxParam: ";
// outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel << " [micro-m]; ";
// outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << denomi ;
outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel*0.0001 << " [micro-m] "; outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel*0.0001 << " [micro-m] ";
outputProbFile << 4*M_PI*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi << " [A²] \n"; outputProbFile << 4*M_PI*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi << " [A²] \n";
} }
...@@ -867,6 +866,7 @@ int bioem::run() ...@@ -867,6 +866,7 @@ int bioem::run()
} }
//************* Writing Cross-Correlations if requiered //************* Writing Cross-Correlations if requiered
//************* This is currently not in the manual *****
if(param.param_device.writeCC){ if(param.param_device.writeCC){
int cc=0; int cc=0;
...@@ -981,6 +981,7 @@ int bioem::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha, ...@@ -981,6 +981,7 @@ int bioem::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha,
//*************************************************************************************** //***************************************************************************************
//***** BioEM routine for comparing reference maps to convoluted maps ***** //***** BioEM routine for comparing reference maps to convoluted maps *****
//***************************************************************************************
if (FFTAlgo) if (FFTAlgo)
{ {
//With FFT Algorithm //With FFT Algorithm
...@@ -997,7 +998,7 @@ int bioem::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha, ...@@ -997,7 +998,7 @@ int bioem::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha,
#pragma omp parallel for schedule(dynamic, 1) #pragma omp parallel for schedule(dynamic, 1)
for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++) for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{ {
compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, conv_map, pProb, param.param_device, RefMap); compareRefMapShifted < -1 > (iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, conv_map, pProb, param.param_device, RefMap);
} }
} }
return(0); return(0);
...@@ -1066,46 +1067,13 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1066,46 +1067,13 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
if(param.doquater){ if(param.doquater){
myfloat_t quater[4]; myfloat_t quater[4];
// myfloat_t cq[4],temp[4];
// myfloat_t a[4];
//quaternion //quaternion
/* quater[1]=param.angles[iMap].pos[1];
quater[2]=param.angles[iMap].pos[2];
quater[3]=param.angles[iMap].pos[0];
quater[0]=param.angles[iMap].quat4;*/
quater[0]=param.angles[iMap].pos[0]; quater[0]=param.angles[iMap].pos[0];
quater[1]=param.angles[iMap].pos[1]; quater[1]=param.angles[iMap].pos[1];
quater[2]=param.angles[iMap].pos[2]; quater[2]=param.angles[iMap].pos[2];
quater[3]=param.angles[iMap].quat4; quater[3]=param.angles[iMap].quat4;
// its conjugate //Rotation Matrix for Quaterions (wikipeda)
/* cq[0]=quater[0];
cq[1]=-quater[1];
cq[2]=-quater[2];
cq[3]=-quater[3];
for(int n = 0; n < Model.nPointsModel; n++)
{
//Rotation with quaterion
// Initial coords
a[0]=0.;
a[1]=Model.points[n].point.pos[0];
a[2]=Model.points[n].point.pos[1];
a[3]=Model.points[n].point.pos[2];
temp[0]=a[0]*cq[0]-a[1]*cq[1]-a[2]*cq[2]-a[3]*cq[3];
temp[1]=a[0]*cq[1]+a[1]*cq[0]+a[2]*cq[3]-a[3]*cq[2];
temp[2]=a[0]*cq[2]-a[1]*cq[3]+a[2]*cq[0]+a[3]*cq[1];
temp[3]=a[0]*cq[3]+a[1]*cq[2]-a[2]*cq[1]+a[3]*cq[0];
// First coordinate is zero Going back to real space
RotatedPointsModel[n].pos[0] = quater[0]*temp[1]+quater[1]*temp[0]+quater[2]*temp[3]-quater[3]*temp[2];
RotatedPointsModel[n].pos[1] = quater[0]*temp[2]-quater[1]*temp[3]+quater[2]*temp[0]+quater[3]*temp[1];
RotatedPointsModel[n].pos[2] = quater[0]*temp[3]+quater[1]*temp[2]-quater[2]*temp[1]+quater[3]*temp[0];
}*/
//Rotation Matrix for Quaterions (wikipeda)
rotmat[0][0] = 1- 2 * quater[1] * quater[1] - 2 * quater[2] * quater[2]; rotmat[0][0] = 1- 2 * quater[1] * quater[1] - 2 * quater[2] * quater[2];
rotmat[1][0] = 2 * ( quater[0] * quater[1] - quater[2] * quater[3]); rotmat[1][0] = 2 * ( quater[0] * quater[1] - quater[2] * quater[3]);
rotmat[2][0] = 2 * ( quater[0] * quater[2] + quater[1] * quater[3]); rotmat[2][0] = 2 * ( quater[0] * quater[2] + quater[1] * quater[3]);
...@@ -1117,10 +1085,6 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1117,10 +1085,6 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
rotmat[2][2] = 1- 2 * quater[0] * quater[0] - 2 * quater[1] * quater[1]; rotmat[2][2] = 1- 2 * quater[0] * quater[0] - 2 * quater[1] * quater[1];
/*Debug for(int i =0; i <3 ; i++){
for(int j =0; j <3 ; j++){cout << "MAT "<< i << " " << j << " " << rotmat[i][j] << "\n";
} } */
} else{ } else{
//*************** Euler Angles**************************** //*************** Euler Angles****************************
...@@ -1249,7 +1213,6 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1249,7 +1213,6 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
localproj[ i * param.param_device.NumberPixels + j] *= ratioDen; localproj[ i * param.param_device.NumberPixels + j] *= ratioDen;
} }
} }
// cout << "HERE " << tempden << " " << Model.NormDen << " \n";
// **** Output Just to check**** // **** Output Just to check****
#ifdef DEBUG #ifdef DEBUG
......
bioem_algorithm.h
View file @ 01a19a3d
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images> < BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer. Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
Volker Lindenstruth and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany. Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution. See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
...@@ -24,7 +25,8 @@ template <int GPUAlgo> ...@@ -24,7 +25,8 @@ 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, bool doAngle, 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)
{ {
//**************** Routine to perform the numerical BioEM intergal *********** // *********** Not using FFT ALGORITHM ******************
//*********** Routine to perform the numerical BioEM intergal ***********
// ******* Summing total Probabilities ************* // ******* Summing total Probabilities *************
...@@ -49,8 +51,8 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef ...@@ -49,8 +51,8 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef
} }
pProbMap.max.max_prob_orient = iOrient; pProbMap.max.max_prob_orient = iOrient;
pProbMap.max.max_prob_conv = iConv; pProbMap.max.max_prob_conv = iConv;
// pProbMap.max.max_prob_norm = - ( -sumC * RefMap.sum_RefMap[iRefMap] + param.Ntotpi * value ) / ( sumC * sumC - sumsquareC * param.Ntotpi); // pProbMap.max.max_prob_norm = - ( -sumC * RefMap.sum_RefMap[iRefMap] + param.Ntotpi * value ) / ( sumC * sumC - sumsquareC * param.Ntotpi);
// pProbMap.max.max_prob_mu = - ( -sumC * value + sumsquareC * RefMap.sum_RefMap[iRefMap] ) / ( sumC * sumC - sumsquareC * param.Ntotpi); // pProbMap.max.max_prob_mu = - ( -sumC * value + sumsquareC * RefMap.sum_RefMap[iRefMap] ) / ( sumC * sumC - sumsquareC * param.Ntotpi);
} }
if (GPUAlgo != 2) pProbMap.Total += exp(logpro - pProbMap.Constoadd); if (GPUAlgo != 2) pProbMap.Total += exp(logpro - pProbMap.Constoadd);
...@@ -69,11 +71,11 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef ...@@ -69,11 +71,11 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef
} }
} }
__device__ static inline myfloat_t calc_logpro(const bioem_param_device& param, const myfloat_t sum, const myfloat_t sumsquare, const myfloat_t crossproMapConv, const myfloat_t sumref, const myfloat_t sumsquareref) __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param, const myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t sum, const myfloat_t sumsquare, const myfloat_t crossproMapConv, const myfloat_t sumref, const myfloat_t sumsquareref)
{ {
// ********** Routine to calculate the log posterior of Eq. of Pmw in SI of JSB paper //*** MAIN ROUTINE TO CALCULATE THE LOGPRO FOR ALL KERNELS*************//
// **** calculate the log posterior of Eq. of Pmw in SI of JSB paper ***//
// Related to Reference calculated Projection // Related to Reference calculated Projection
const myfloat_t ForLogProb = (sumsquare * param.Ntotpi - sum * sum); const myfloat_t ForLogProb = (sumsquare * param.Ntotpi - sum * sum);
...@@ -84,34 +86,32 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param, ...@@ -84,34 +86,32 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
/// ******* 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);
// if(param.debugterm){ myfloat_t logpro = (3 - param.Ntotpi) * 0.5 * log(firstele) + (param.Ntotpi * 0.5 - 2) * log((param.Ntotpi - 2) * ForLogProb);
// printf("Separate cc: %f c: %f oo: %f o: %f co: %f logP: %f\n",sumsquare,sum,sumsquareref,sumref,crossproMapConv, logpro );
// } //*************Adding Gaussian Prior to envelope & Defocus parameter******************
if(not param.tousepsf){
logpro = logpro - env * env / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf -
(pha - param.Priordefcent ) * (pha - param.Priordefcent ) / 2. / param.sigmaPriordefo / param.sigmaPriordefo ;
} else {
myfloat_t envF,phaF;
envF = 4.* M_PI * M_PI * env / ( env * env + pha * pha) ;
phaF = 4.* M_PI * M_PI * pha / ( env * env + pha * pha);
logpro = logpro - envF * envF / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf - (phaF - param.Priordefcent ) * (phaF - param.Priordefcent ) / 2. / param.sigmaPriordefo / param.sigmaPriordefo ;
}
return(logpro); return(logpro);
} }
__device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myfloat_t amp, myfloat_t pha, myfloat_t env, myfloat_t sumC, myfloat_t sumsquareC, myfloat_t 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 amp, myfloat_t pha, myfloat_t env, myfloat_t sumC, myfloat_t sumsquareC, myfloat_t value, int disx, int disy, bioem_Probability& pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
{ {
// MOST IMPORTANT ROUTINE // IMPORTANT ROUTINE Summation of LogProb using FFTALGO
// ******************************************************** // ********************************************************
// *********** Calculates the BioEM probability *********** // *********** Calculates the BioEM probability ***********
// ******************************************************** // ********************************************************
myfloat_t logpro = calc_logpro(param, sumC, sumsquareC, value, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]); myfloat_t logpro = calc_logpro(param, amp, pha, env, sumC, sumsquareC, value, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]);
//Adding Gaussian Prior to envelope parameter
if(not param.tousepsf){
logpro = logpro - env * env / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf - (pha - param.Priordefcent ) * (pha - param.Priordefcent ) / 2. / param.sigmaPriordefo / param.sigmaPriordefo ;
} else {
myfloat_t envF,phaF;
envF = 4.* M_PI * M_PI * env / ( env * env + pha * pha) ;
phaF = 4.* M_PI * M_PI * pha / ( env * env + pha * pha);
logpro = logpro - envF * envF / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf - (phaF - param.Priordefcent ) * (phaF - param.Priordefcent ) / 2. / param.sigmaPriordefo / param.sigmaPriordefo ;
}
//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, param.writeAngles); //update_prob<-1>(logpro, iRefMap, iOrient, iConv, disx, disy, pProb, param.writeAngles);
...@@ -153,6 +153,7 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo ...@@ -153,6 +153,7 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
__device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient, const int iConv, const myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t* lCC, const myfloat_t sumC, const myfloat_t sumsquareC, bioem_Probability& pProb, const bioem_param_device& param, const bioem_RefMap& RefMap) __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient, const int iConv, const myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t* lCC, const myfloat_t sumC, const myfloat_t sumsquareC, bioem_Probability& pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
{ {
//******************** Using FFT algorithm **************************
//******************* Get cross-crollation of Ical to Iobs ******************* //******************* Get cross-crollation of Ical to Iobs *******************
//*********** Routine to get the Cross-Corellation from lCC for the interested center displacement ************* //*********** Routine to get the Cross-Corellation from lCC for the interested center displacement *************
...@@ -180,6 +181,7 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient, ...@@ -180,6 +181,7 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient,
} }
} }
//************ The following if is not in the manual***********
if (param.writeCC) if (param.writeCC)
{ {
// If the Cross-correlation is to be written out and stored using Bayesian analysis // If the Cross-correlation is to be written out and stored using Bayesian analysis
...@@ -226,7 +228,8 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient, ...@@ -226,7 +228,8 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient,
} }
template <int GPUAlgo, class RefT> template <int GPUAlgo, class RefT>
__device__ static inline void compareRefMap(const int iRefMap, const int iOrient, const int iConv, const myfloat_t* 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 myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t sumC,
const myfloat_t sumsquareC, const myfloat_t* 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)
{ {
...@@ -314,7 +317,7 @@ template <int GPUAlgo, class RefT> ...@@ -314,7 +317,7 @@ template <int GPUAlgo, class RefT>
#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, amp, pha, env, sum, sumsquare, crossproMapConv, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]);
} }
else else
{ {
...@@ -485,13 +488,13 @@ template <int GPUAlgo, class RefT> ...@@ -485,13 +488,13 @@ template <int GPUAlgo, class RefT>
} }
template <int GPUAlgo, class RefT> template <int GPUAlgo, class RefT>
__device__ static inline void compareRefMapShifted(const int iRefMap, const int iOrient, const int iConv, const myfloat_t* Mapconv, bioem_Probability& pProb, const bioem_param_device& param, const RefT& RefMap) __device__ static inline void compareRefMapShifted(const int iRefMap, const int iOrient, const int iConv, const myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t sumC, const myfloat_t sumsquareC, const myfloat_t* Mapconv, bioem_Probability& pProb, const bioem_param_device& param, const RefT& RefMap)
{ {
for (int cent_x = -param.maxDisplaceCenter; cent_x <= param.maxDisplaceCenter; cent_x = cent_x + param.GridSpaceCenter) for (int cent_x = -param.maxDisplaceCenter; cent_x <= param.maxDisplaceCenter; cent_x = cent_x + param.GridSpaceCenter)
{ {
for (int cent_y = -param.maxDisplaceCenter; cent_y <= param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter) for (int cent_y = -param.maxDisplaceCenter; cent_y <= param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
{ {
compareRefMap<GPUAlgo>(iRefMap, iOrient, iConv, Mapconv, pProb, param, RefMap, cent_x, cent_y); compareRefMap<GPUAlgo>(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, Mapconv, pProb, param, RefMap, cent_x, cent_y);
} }
} }
} }
......
bioem_cuda.cu
View file @ 01a19a3d
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images> < BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer. Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
Max Planck Institute of Biophysics, Frankfurt, Germany. Volker Lindenstruth and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution. See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
...@@ -80,21 +81,23 @@ bioem_cuda::~bioem_cuda() ...@@ -80,21 +81,23 @@ bioem_cuda::~bioem_cuda()
deviceExit(); deviceExit();
} }
__global__ void compareRefMap_kernel(const int iOrient, const int iConv, const myfloat_t* pMap, bioem_Probability pProb, const bioem_param_device param, const bioem_RefMap_Mod RefMap, const int cent_x, const int cent_y, const int maxRef) __global__ void compareRefMap_kernel(const int iOrient, const int iConv, const myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t sumC,
const myfloat_t sumsquareC, const myfloat_t* pMap, bioem_Probability pProb,
const bioem_param_device param, const bioem_RefMap_Mod RefMap, const int cent_x, const int cent_y, const int maxRef)
{ {
const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX; const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
if (iRefMap < maxRef) if (iRefMap < maxRef)
{ {
compareRefMap<0>(iRefMap, iOrient, iConv, pMap, pProb, param, RefMap, cent_x, cent_y); compareRefMap<0>(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, pMap, pProb, param, RefMap, cent_x, cent_y);
} }
} }
__global__ void compareRefMapShifted_kernel(const int iOrient, const int iConv, const myfloat_t* pMap, bioem_Probability pProb, const bioem_param_device param, const bioem_RefMap_Mod RefMap, const int maxRef) __global__ void compareRefMapShifted_kernel(const int iOrient, const int iConv, const myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t sumC, const myfloat_t sumsquareC, nst myfloat_t* pMap, bioem_Probability pProb, const bioem_param_device param, const bioem_RefMap_Mod RefMap, const int maxRef)
{ {
const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX; const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
if (iRefMap < maxRef) if (iRefMap < maxRef)
{ {
compareRefMapShifted<1>(iRefMap, iOrient, iConv, pMap, pProb, param, RefMap); compareRefMapShifted<1>(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, pMap, pProb, param, RefMap);
} }
} }
...@@ -250,7 +253,7 @@ int bioem_cuda::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t ...@@ -250,7 +253,7 @@ int bioem_cuda::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t
} }
else if (GPUAlgo == 0) //All shifts in one kernel else if (GPUAlgo == 0) //All shifts in one kernel
{ {
compareRefMapShifted_kernel<<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream[0]>>> (iOrient, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, *pRefMap_device_Mod, maxRef); compareRefMapShifted_kernel<<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream[0]>>> (iOrient, iConv, amp, pha, env, sumC, sumsquareC, pConvMap_device[iConv & 1], pProb_device, param.param_device, *pRefMap_device_Mod, maxRef);
} }
else else
{ {
......
include/bioem.h
View file @ 01a19a3d
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
Volker Lindenstruth and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
#ifndef BIOEM_H #ifndef BIOEM_H
#define BIOEM_H #define BIOEM_H
......
include/bioem_cuda.h
View file @ 01a19a3d
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
Volker Lindenstruth and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
#ifndef BIOEM_CUDA_H #ifndef BIOEM_CUDA_H
#define BIOEM_CUDA_H #define BIOEM_CUDA_H
......
include/bioem_cuda_internal.h
View file @ 01a19a3d
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
< BioEM software for Bayesian inference of Electron Microscopy images>
Copyright (C) 2016 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
Volker Lindenstruth and Gerhard Hummer.
Max Planck Institute of Biophysics, Frankfurt, Germany.
See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
#ifndef BIOEM_CUDA_INTERNAL_H #ifndef BIOEM_CUDA_INTERNAL_H
#define BIOEM_CUDA_INTERNAL_H