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

Prior on defocus

parent af0a4c90
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
See license statement for terms of distribution. See license statement for terms of distribution.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
#ifdef WITH_MPI
#include <mpi.h> #include <mpi.h>
#define MPI_CHK(expr) \ #define MPI_CHK(expr) \
...@@ -14,9 +14,13 @@ ...@@ -14,9 +14,13 @@
{ \ { \
fprintf(stderr, "Error in MPI function %s: %d\n", __FILE__, __LINE__); \ fprintf(stderr, "Error in MPI function %s: %d\n", __FILE__, __LINE__); \
} }
#endif
#include <fstream> #include <fstream>
#include <boost/program_options.hpp> #include <boost/program_options.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <iostream> #include <iostream>
#include <algorithm> #include <algorithm>
#include <iterator> #include <iterator>
...@@ -67,6 +71,7 @@ const int num_colors = sizeof(colors)/sizeof(colors[0]); ...@@ -67,6 +71,7 @@ const int num_colors = sizeof(colors)/sizeof(colors[0]);
using namespace boost; using namespace boost;
namespace po = boost::program_options; namespace po = boost::program_options;
namespace bran= boost::random;
using namespace std; using namespace std;
...@@ -104,6 +109,8 @@ int bioem::configure(int ac, char* av[]) ...@@ -104,6 +109,8 @@ int bioem::configure(int ac, char* av[])
HighResTimer timer; HighResTimer timer;
std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap; std::string infile, modelfile, mapfile,Inputanglefile,Inputbestmap;
if (mpi_rank == 0) if (mpi_rank == 0)
{ {
...@@ -309,7 +316,7 @@ int bioem::configure(int ac, char* av[]) ...@@ -309,7 +316,7 @@ int bioem::configure(int ac, char* av[])
if (mpi_rank != 0)param.angles = (myfloat3_t*) mallocchk(param.nTotGridAngles * sizeof (myfloat3_t)); if (mpi_rank != 0)param.angles = (myfloat3_t*) mallocchk(param.nTotGridAngles * sizeof (myfloat3_t));
MPI_Bcast(param.angles, param.nTotGridAngles * sizeof (myfloat3_t),MPI_BYTE, 0, MPI_COMM_WORLD); MPI_Bcast(param.angles, param.nTotGridAngles * sizeof (myfloat3_t),MPI_BYTE, 0, MPI_COMM_WORLD);
#ifdef PILAR_DEBUG #ifdef DEBUG
for(int n=0;n<param.nTotGridAngles;n++){ for(int n=0;n<param.nTotGridAngles;n++){
cout << "CHECK: Angle orient " << mpi_rank << " "<< n << " " << param.angles[n].pos[0] << " " << param.angles[n].pos[1] << " " << param.angles[n].pos[2] << " " << param.angles[n].quat4 << " " << "\n";} cout << "CHECK: Angle orient " << mpi_rank << " "<< n << " " << param.angles[n].pos[0] << " " << param.angles[n].pos[1] << " " << param.angles[n].pos[2] << " " << param.angles[n].quat4 << " " << "\n";}
...@@ -614,7 +621,7 @@ int bioem::run() ...@@ -614,7 +621,7 @@ int bioem::run()
for (int i = 0;i < RefMap.ntotRefMap;i++) for (int i = 0;i < RefMap.ntotRefMap;i++)
{ {
bioem_Probability_map& pProbMap = pProb.getProbMap(i); bioem_Probability_map& pProbMap = pProb.getProbMap(i);
#ifdef PILAR_DEBUG #ifdef DEBUG
cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd << "\n"; cout << "Reduction " << mpi_rank << " Map " << i << " Prob " << pProbMap.Total << " Const " << pProbMap.Constoadd << "\n";
#endif #endif
tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]); tmp1[i] = pProbMap.Total * exp(pProbMap.Constoadd - tmp2[i]);
...@@ -1021,7 +1028,7 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_ ...@@ -1021,7 +1028,7 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_
doRefMapFFT(iRefMap, iOrient, iConv, amp, pha, env, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap); doRefMapFFT(iRefMap, iOrient, iConv, amp, pha, env, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
#ifdef PILAR_DEBUG #ifdef DEBUG
if (param.param_device.writeCC) if (param.param_device.writeCC)
{ int cc=0; { int cc=0;
for (int cent_x = 0; cent_x < param.param_device.NumberPixels ; 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)
...@@ -1105,7 +1112,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1105,7 +1112,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
gam = param.angles[iMap].pos[2]; gam = param.angles[iMap].pos[2];
//*** To see how things are going: //*** To see how things are going:
#ifdef PILAR_DEBUG #ifdef DEBUG
cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n";
#endif #endif
// ********** Creat Rotation with pre-defiend grid of orientations********** // ********** Creat Rotation with pre-defiend grid of orientations**********
...@@ -1151,7 +1158,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1151,7 +1158,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
for(int n = 0; n < Model.nPointsModel; n++) for(int n = 0; n < Model.nPointsModel; n++)
{ {
if(Model.points[n].radius < param.pixelSize){ if(Model.points[n].radius <= param.pixelSize){
// cout << "Radius less than Pixel size: use keyword NO_PROJECT_RADIUS in inputfile\n"; // cout << "Radius less than Pixel size: use keyword NO_PROJECT_RADIUS in inputfile\n";
i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f); i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f); j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
...@@ -1207,7 +1214,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1207,7 +1214,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
// **** Output Just to check**** // **** Output Just to check****
#ifdef PILAR_DEBUG #ifdef DEBUG
// if(iMap == 0) // if(iMap == 0)
{ {
ofstream myexamplemap; ofstream myexamplemap;
...@@ -1317,7 +1324,14 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -1317,7 +1324,14 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
if (mpi_rank == 0 && param.printModel) if (mpi_rank == 0 && param.printModel)
{ {
MTRand mtr; // MTRand mtr;
bran::mt19937 gen;
//Generating random seed so the maps do not have correlated Noise
gen.seed(static_cast<unsigned int>(std::time(0)));
//Uniform Noise: bran::uniform_int_distribution<> dist(1, 6);
//Gaussian noise
bran::normal_distribution <> distn(0.0,param.stnoise);
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);
...@@ -1334,7 +1348,8 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj ...@@ -1334,7 +1348,8 @@ int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj
if(!param.withnoise){ if(!param.withnoise){
myexamplemap << "\nMAP " << k+param.ddx << " " << j+param.ddy << " " << Mapconv[k * param.param_device.NumberPixels + j] / norm2 *param.bestnorm +param.bestoff ; myexamplemap << "\nMAP " << k+param.ddx << " " << j+param.ddy << " " << Mapconv[k * param.param_device.NumberPixels + j] / norm2 *param.bestnorm +param.bestoff ;
} else { } else {
myexamplemap << "\nMAP " << k+param.ddx << " " << j+param.ddy << " " << Mapconv[k * param.param_device.NumberPixels + j] / norm2 *param.bestnorm +param.bestoff+ mtr.randNorm(0.0,param.stnoise); myexamplemap << "\nMAP " << k+param.ddx << " " << j+param.ddy << " " << Mapconv[k * param.param_device.NumberPixels + j] / norm2 *param.bestnorm +param.bestoff+distn(gen);
// cout << distn(gen) << "CHECK\n";
} }
} }
myexamplemap << " \n"; myexamplemap << " \n";
......
...@@ -105,11 +105,12 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo ...@@ -105,11 +105,12 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
//Adding Gaussian Prior to envelope parameter //Adding Gaussian Prior to envelope parameter
if(not param.tousepsf){ if(not param.tousepsf){
logpro = logpro - env * env / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf ; logpro = logpro - env * env / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf - (pha - param.Priordefcent ) * (pha - param.Priordefcent ) / 2. / param.sigmaPriordefo / param.sigmaPriordefo ;
} else { } else {
myfloat_t envF; myfloat_t envF,phaF;
envF = 4.* M_PI * M_PI * env / ( env * env + pha * pha) ; envF = 4.* M_PI * M_PI * env / ( env * env + pha * pha) ;
logpro = logpro - envF * envF / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf ; 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
......
...@@ -22,6 +22,9 @@ public: ...@@ -22,6 +22,9 @@ public:
myfloat_t Ntotpi; myfloat_t Ntotpi;
myfloat_t volu; myfloat_t volu;
myfloat_t sigmaPriorbctf; myfloat_t sigmaPriorbctf;
myfloat_t sigmaPriordefo;
myfloat_t Priordefcent;
// If to write Probabilities of Angles from Model // If to write Probabilities of Angles from Model
bool writeAngles; bool writeAngles;
......
#ifdef WITH_MPI
#include <mpi.h> #include <mpi.h>
#define MPI_CHK(expr) \ #define MPI_CHK(expr) \
...@@ -5,6 +6,7 @@ ...@@ -5,6 +6,7 @@
{ \ { \
fprintf(stderr, "Error in MPI function %s: %d\n", __FILE__, __LINE__); \ fprintf(stderr, "Error in MPI function %s: %d\n", __FILE__, __LINE__); \
} }
#endif
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
......
...@@ -83,7 +83,7 @@ int bioem_param::readParameters(const char* fileinput) ...@@ -83,7 +83,7 @@ int bioem_param::readParameters(const char* fileinput)
writeCTF=false; writeCTF=false;
elecwavel=0.019866; elecwavel=0.019866;
ignoreCCoff=false; ignoreCCoff=false;
doquater= false; doquater=false;
nocentermass=false; nocentermass=false;
printrotmod=false; printrotmod=false;
readquatlist=false; readquatlist=false;
...@@ -98,7 +98,8 @@ int bioem_param::readParameters(const char* fileinput) ...@@ -98,7 +98,8 @@ int bioem_param::readParameters(const char* fileinput)
shiftX=0; shiftX=0;
shiftY=0; shiftY=0;
param_device.sigmaPriorbctf=100.; param_device.sigmaPriorbctf=100.;
param_device.sigmaPriordefo=1.0;
param_device.Priordefcent=3.0;
ifstream input(fileinput); ifstream input(fileinput);
if (!input.good()) if (!input.good())
...@@ -390,8 +391,21 @@ int bioem_param::readParameters(const char* fileinput) ...@@ -390,8 +391,21 @@ int bioem_param::readParameters(const char* fileinput)
{ {
token = strtok(NULL, " "); token = strtok(NULL, " ");
param_device.sigmaPriorbctf=atof(token); param_device.sigmaPriorbctf=atof(token);
cout << "Chainging Gaussian width in Prior of Envelope b parameter" << param_device.sigmaPriorbctf << "\n"; cout << "Chainging Gaussian width in Prior of Envelope b parameter: " << param_device.sigmaPriorbctf << "\n";
} }
else if (strcmp(token, "SIGMA_PRIOR_DEFOCUS") == 0)
{
token = strtok(NULL, " ");
param_device.sigmaPriordefo=atof(token);
cout << "Gaussian Width in Prior of defocus parameter: " << param_device.sigmaPriordefo << "\n";
}
else if (strcmp(token, "PRIOR_DEFOCUS_CENTER") == 0)
{
token = strtok(NULL, " ");
param_device.Priordefcent=atof(token);
cout << "Gaussian Center in Prior of defocus parameter: " << param_device.Priordefcent << "\n";
}
else if (strcmp(token, "IGNORE_POINTSOUT") == 0) else if (strcmp(token, "IGNORE_POINTSOUT") == 0)
{ {
ignorepointsout=true; ignorepointsout=true;
...@@ -437,6 +451,8 @@ int bioem_param::readParameters(const char* fileinput) ...@@ -437,6 +451,8 @@ int bioem_param::readParameters(const char* fileinput)
//Asigning values of envelope according to b-envelope (not b-factor) //Asigning values of envelope according to b-envelope (not b-factor)
startGridEnvelop = startBfactor ;// 2.f; startGridEnvelop = startBfactor ;// 2.f;
endGridEnvelop = endBfactor ; // / 2.f; endGridEnvelop = endBfactor ; // / 2.f;
param_device.Priordefcent *= M_PI * 2.f * 10000 * elecwavel ;
param_device.sigmaPriordefo *= M_PI * 2.f * 10000 * elecwavel ;
} }
if(elecwavel==0.019688)cout << "Using default electron wave length: 0.019688 (A) of 300kV microscope\n"; if(elecwavel==0.019688)cout << "Using default electron wave length: 0.019688 (A) of 300kV microscope\n";
...@@ -469,12 +485,11 @@ int bioem_param::forprintBest(const char* fileinput) ...@@ -469,12 +485,11 @@ int bioem_param::forprintBest(const char* fileinput)
writeCTF=false; writeCTF=false;
elecwavel=0.019866; elecwavel=0.019866;
ignoreCCoff=false; ignoreCCoff=false;
doquater= false; doquater=false;
nocentermass=false; nocentermass=false;
printrotmod=false; printrotmod=false;
readquatlist=false; readquatlist=false;
doaaradius=true; doaaradius=true;
doquater=true;
shiftX=0; shiftX=0;
shiftY=0; shiftY=0;
...@@ -890,7 +905,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN ...@@ -890,7 +905,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
angles[n].pos[1] = (myfloat_t) b; angles[n].pos[1] = (myfloat_t) b;
angles[n].pos[2] = (myfloat_t) g; angles[n].pos[2] = (myfloat_t) g;
angles[n].quat4 =0.0;//just to be sure */ angles[n].quat4 =0.0;//just to be sure */
#ifdef PILAR_DEBUG #ifdef DEBUG
if(yespriorAngles) cout << "check orient: " << n << " " << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior: " << angprior[n]<< "\n"; if(yespriorAngles) cout << "check orient: " << n << " " << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior: " << angprior[n]<< "\n";
#endif #endif
} }
...@@ -1049,7 +1064,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN ...@@ -1049,7 +1064,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
sscanf (tmpVals, "%f", &pp); sscanf (tmpVals, "%f", &pp);
if(pp <0.0000001)cout << "Sure you're input is correct? Very small prior.\n"; if(pp <0.0000001)cout << "Sure you're input is correct? Very small prior.\n";
angprior[n] = pp; angprior[n] = pp;
#ifdef PILAR_DEBUG #ifdef DEBUG
if(yespriorAngles) cout << "check orient: " << n << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior: " << angprior[n] << "\n"; if(yespriorAngles) cout << "check orient: " << n << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior: " << angprior[n] << "\n";
#endif #endif
} }
......
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