Commit 8ae18d67 authored by Pilar Cossio's avatar Pilar Cossio

Corrected and debugged Sep.29.15

parent 14dd2863
##### Micrograph Parameters #######
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
##### Euler angle grid points: #######
GRIDPOINTS_ALPHA 10
GRIDPOINTS_BETA 5
##### Constrast transfer integration: #######
CTF_B_FACTOR 100.0 300.5 8
CTF_DEFOCUS 1.0 6.0 5
CTF_AMPLITUDE 0.1 0.1 1
##### Center displacement: #######
DISPLACE_CENTER 10 2
#### Write probability as function of Orientations ####
WRITE_CROSSCOR
CROSSCOR_NOTBAYESIAN
FLIPPED
##### Micrograph Parameters #######
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
##### Euler angle grid points: #######
GRIDPOINTS_ALPHA 3
GRIDPOINTS_BETA 2
##### PSF integration: #######
USE_PSF
PSF_ENVELOPE 0.0002 0.0008 3
PSF_PHASE 0.002 0.008 3
PSF_AMPLITUDE 0.5 1.0 1
##### Center displacement: #######
DISPLACE_CENTER 10 2
##### Micrograph Parameters #######
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
##### Euler angle grid points: #######
GRIDPOINTS_ALPHA 2
GRIDPOINTS_BETA 2
##### Constrast transfer integration: #######
CTF_B_FACTOR 100.0 300.5 2
CTF_DEFOCUS 1.0 6.0 5
CTF_AMPLITUDE 0.1 0.1 1
##### Center displacement: #######
DISPLACE_CENTER 10 2
####### Priors ###############
PRIOR_MODEL 0.005
PRIOR_ANGLES
##### Micrograph Parameters #######
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
USE_QUATERNIONS
##### Constrast transfer integration: #######
CTF_B_FACTOR 100.0 300.5 2
CTF_DEFOCUS 1.0 6.0 5
CTF_AMPLITUDE 0.1 0.1 1
##### Center displacement: #######
DISPLACE_CENTER 10 2
####### Priors ###############
PRIOR_MODEL 0.005
PRIOR_ANGLES
##### Micrograph Parameters #######
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
##### Quaterion grid points: #######
USE_QUATERNIONS
GRIDPOINTS_QUATERNION 5
##### Constrast transfer integration: #######
CTF_B_FACTOR 100.0 300.5 8
CTF_DEFOCUS 1.0 6.0 5
CTF_AMPLITUDE 0.1 0.1 1
##### Center displacement: #######
DISPLACE_CENTER 10 2
##### Micrograph Parameters #######
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
##### Euler angle grid points: #######
GRIDPOINTS_ALPHA 10
GRIDPOINTS_BETA 5
##### Constrast transfer integration: #######
CTF_B_FACTOR 100.0 300.5 8
CTF_DEFOCUS 1.0 6.0 5
CTF_AMPLITUDE 0.1 0.1 1
##### Center displacement: #######
DISPLACE_CENTER 10 2
#### Write probability as function of Orientations ####
WRITE_PROB_ANGLES
This diff is collapsed.
......@@ -44,6 +44,8 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef
}
pProbMap.max.max_prob_orient = iOrient;
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_mu = - ( -sumC * value + sumsquareC * RefMap.sum_RefMap[iRefMap] ) / ( sumC * sumC - sumsquareC * param.Ntotpi);
}
if (GPUAlgo != 2) pProbMap.Total += exp(logpro - pProbMap.Constoadd);
......@@ -51,6 +53,7 @@ __device__ static inline void update_prob(const myfloat_t logpro, const int iRef
{
bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
// printf("Here\n");
//Summing probabilities for each orientation
if(pProbAngle.ConstAngle < logpro)
{
......@@ -76,9 +79,9 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
const myfloat_t logpro = (3 - param.Ntotpi) * 0.5 * log(firstele) + (param.Ntotpi * 0.5 - 2) * log((param.Ntotpi - 2) * ForLogProb);
if(param.debugterm){
// if(param.debugterm){
// printf("Separate cc: %f c: %f oo: %f o: %f co: %f logP: %f\n",sumsquare,sum,sumsquareref,sumref,crossproMapConv, logpro );
}
// }
return(logpro);
}
......@@ -90,9 +93,12 @@ __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]);
//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);
bioem_Probability_map& pProbMap = pProb.getProbMap(iRefMap);
if(pProbMap.Constoadd < logpro)
{
pProbMap.Total = pProbMap.Total * exp(-logpro + pProbMap.Constoadd);
......@@ -105,17 +111,18 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
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_mu = - ( -sumC * value + sumsquareC * RefMap.sum_RefMap[iRefMap] ) / ( sumC * sumC - sumsquareC * param.Ntotpi);
}
pProbMap.Total += exp(logpro - pProbMap.Constoadd);
if(param.debugterm){
// printf("Separate Ptot: %f Const: %f logProb %f \n",pProbMap.Total,pProbMap.Constoadd,logpro);
printf("Separate Ptot: %f Const: %f logProb %f %f %f \n",pProbMap.Total,pProbMap.Constoadd,logpro,pProbMap.max.max_prob_norm,pProbMap.max.max_prob_mu);
}
if (param.writeAngles)
{
bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
// if(iConv==0 && disx==0 && disy==0)
// if(iOrient==1)printf("Separate Ptot: %f Const: %f logProb %f param: %d %d %d \n",logpro,pProbAngle.ConstAngle,pProbAngle.forAngles,disx,disx,iOrient);
if(pProbAngle.ConstAngle < logpro)
{
......@@ -123,6 +130,7 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
pProbAngle.ConstAngle = logpro;
}
pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle);
// if(iOrient==5)printf("After separate Ptot: %f Const: %f logProb %f \n",logpro,pProbAngle.ConstAngle,pProbAngle.forAngles);
}
}
......@@ -151,6 +159,7 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient,
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x - param.NumberPixels, cent_y - param.NumberPixels, pProb, param, RefMap);
}
}
// printf("HERE!!!");
if (param.writeCC)
{
......@@ -186,7 +195,7 @@ __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient,
pProbCC.forCC += exp(ttmp - pProbCC.ConstCC);
}
// printf("Separate Ptot: %f Const: %f logProb %f \n",pProbCC.forCC,pProbCC.ConstCC,ttmp);
// printf("Separate %d %d Ptot: %f Const: %f logProb %f \n",cent_x,cent_y,pProbCC.forCC,pProbCC.ConstCC,ttmp);
cc++;
}
......
......@@ -50,6 +50,7 @@ struct myfloat3_t
{
myfloat_t pos[BIOEM_FLOAT_3_PHYSICAL_SIZE];
myfloat_t quat4;
// myfloat_t prior;
};
#ifdef BIOEM_GPUCODE
......
......@@ -100,6 +100,7 @@ class bioem_Probability_angle
public:
myfloat_t forAngles;
myfloat_t ConstAngle;
myfloat_t priorang;
};
class bioem_Probability_cc
......
......@@ -38,6 +38,7 @@ public:
bioem_param();
~bioem_param();
int readParameters(const char* fileinput,const char* fileangles);
int CalculateGridsParam();
int CalculateRefCTF();
......@@ -50,7 +51,10 @@ public:
bool printrotmod;
bool readquatlist;
bool showrotatemod;
bool notsqure;
bool notnormmap;
bool usepsf;
myfloat_t elecwavel;
bioem_param_device param_device;
......@@ -63,6 +67,11 @@ public:
size_t getCtfParamCount() {return nTotCTFs;}
myfloat_t pixelSize;
// Priors
myfloat_t priorMod;
bool yespriorAngles;
myfloat_t* angprior;
// Grid Points in Euler angles, assuming uniform sampling d_alpha=d_gamma (in 2pi) & cos(beta)=-1,1
int angleGridPointsAlpha;
int angleGridPointsBeta;
......@@ -81,6 +90,13 @@ public:
int numberGridPointsDisplaceCenter;
// Grid sampling for the convolution kernel
// CTF
myfloat_t startBfactor, endBfactor;
int numberBfactor;
myfloat_t startDefocus, endDefocus;
int numberDefocus;
//ENVELOPE
myfloat_t startGridEnvelop;
myfloat_t endGridEnvelop;
......
......@@ -220,7 +220,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
cout << ".";
ntotRefMap = nummap + 1;
maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
cout << "Particle Maps read from Standard File \n";
cout << "Particle Maps read from Standard File: " << ntotRefMap << "\n";
}
......@@ -456,7 +456,7 @@ int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
if(nr!=param.param_device.NumberPixels || nc!=param.param_device.NumberPixels )
{
cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << nc << ", j " << nr << ")" << "\n";
exit(1);
if(!param.notsqure) exit(1);
}
if (ntotRefMap == 0)
......@@ -501,21 +501,26 @@ int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
}
else
{
if(!param.notsqure){
maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = (myfloat_t) currfloat;
st += currfloat;
st2 += currfloat*currfloat;
} else {
if( i > 595 && i < 675 && j > 1250 && j< 1330 && nmap >230 && nmap <310)cout << "map1: " << i << " "<< j << " " << nmap << " " << currfloat <<"\n";
}
}
}
if(param.notsqure)exit(1);
//Normaling maps to zero mean and unit standard deviation
if(!param.notnormmap){
st /= float(nr*nc);
st2 = sqrt(st2 / float(nr * nc) - st * st);
for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){
maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] / st2 - st/st2;
//(nmap+ntotRefMap==300)
//cout <<"MAP:: " << i << " " << j << " " << maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] << "\n";
}
}
}
ntotRefMap += ns ;
// cout << ntotRefMap << "\n";
}
......
This diff is collapsed.
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