Tutorial Corr

Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision

@@ -2,16 +2,13 @@
 NUMBER_PIXELS 224
 PIXEL_SIZE 1.32
 
-###### Sincs its a full atom structure #####
-NO_PROJECT_RADIUS
-
 ##### Quaterion grid points: #######
 USE_QUATERNIONS
 
 ##### Constrast transfer integration: #######
-CTF_B_FACTOR 100.0 300.5 8
+CTF_B_ENV 100.0 300.5 8
 CTF_DEFOCUS 1.0 6.0 5
-CTF_AMPLITUDE 0.1 0.1 1
+CTF_AMPLITUDE 0.01 0.5 3
 
 ##### Center displacement: #######
 DISPLACE_CENTER 10 2

Tutorial_BioEM/MODEL_COMPARISION/subtract_LogP.sh

+paste $1 $2 | grep LogProb: | awk 'BEGIN{cum=0;dif=0}{if($1=="RefMap:"){dif=$4-$10;cum+=dif;printf"MapNum: %6d Mod1: %10.5f Mod2: %10.5f Dif: %10.5f Cum: %10.5f\n",$2,$4,$10,dif,cum}}'

bioem_algorithm.h

@@ -85,13 +85,23 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
 	return(logpro);
 }
 
-__device__ static inline void calProb(int iRefMap, int iOrient, int iConv, 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)
 {
 	// ********************************************************
 	// *********** Calculates the BioEM probability ***********
 	// ********************************************************
 
-	const myfloat_t logpro = calc_logpro(param, sumC, sumsquareC, value, RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]);
+	myfloat_t logpro = calc_logpro(param, 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 ; 
+	} else {
+        myfloat_t envF;
+        envF = 4.* M_PI * M_PI * env / ( env * env + pha * pha) ; 
+	logpro = logpro - envF * envF / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf ;
+	}
 
 
 	//GCC is too stupid to inline properly, so the code is copied here
@@ -115,9 +125,12 @@ __device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myflo
 	pProbMap.Total += exp(logpro - pProbMap.Constoadd);
 
 	if(param.debugterm){
-//         	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);
+         	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);
         }
 
+//       printf("Separate Ptot: %f Const: %f logProb %f %f %f \n",amp,pha,env,logpro,param.sigmaPriorbctf);
+
+//      cout << amp,pha,env,logpro,sumC;
 	if (param.writeAngles)
 	{
 		bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
@@ -135,28 +148,28 @@ __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* 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)
 {
 	for (int cent_x = 0; cent_x <= param.maxDisplaceCenter; cent_x = cent_x + param.GridSpaceCenter)
 	{
 		for (int cent_y = 0; cent_y <= param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y, pProb, param, RefMap);
+			calProb(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y, pProb, param, RefMap);
 		}
 		for (int cent_y = param.NumberPixels - param.maxDisplaceCenter; cent_y < param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y - param.NumberPixels, pProb, param, RefMap);
+			calProb(iRefMap, iOrient, iConv,amp, pha, env, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y - param.NumberPixels, pProb, param, RefMap);
 		}
 	}
 	for (int cent_x = param.NumberPixels - param.maxDisplaceCenter; cent_x < param.NumberPixels; cent_x = cent_x + param.GridSpaceCenter)
 	{
 		for (int cent_y = 0; cent_y < param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			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, pProb, param, RefMap);
+			calProb(iRefMap, iOrient, iConv,amp, pha, env, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x - param.NumberPixels, cent_y, pProb, param, RefMap);
 		}
 		for (int cent_y = param.NumberPixels - param.maxDisplaceCenter; cent_y <= param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			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);
+			calProb(iRefMap, iOrient, iConv, amp, pha, env, 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!!!");