2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/emap_CEvalMAP.cpp $
5 $LastChangedRevision: 85 $
6 $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
11 ===========================================================================
12 Evaluation of the MAP algorithms (the results are written to "user1")
13 ===========================================================================
27 void emap_class::EvalMAP(REAL* recon, INTEGER iter)
29 //------ local variables
30 INTEGER i, j, iadr, nb, np, nr, numb, list[4096], usnr;
31 REAL snorm, weight[4096], diff1;
32 double sum1, sum2, y, diffnbs, totalsum, sclf, smx, pictk, ssum1, ssum2;
36 //------ log-likelihood
38 sclf = MAPEM.sclfactor;
39 for (i = 0; i < GeoPar.prjnum; i++)
42 for (j = 0; j < GeoPar.usrays; j++)
46 nr = usnr + GeoPar.fusray;
48 y = MAPEM.prj[(np) * GeoPar.usrays + usnr] * sclf;
53 wray(np, nr, list, weight, &numb, &snorm);
57 Blob.bwray(np, nr, list, weight, &numb, &snorm);
62 for (nb = 0; nb < numb; nb++)
66 smx += pictk * weight[nb];
73 ssum1 += y * log(smx) - smx;
77 fprintf(output, "\n In EvalMAP smx = %14.7f y = %14.7f", smx,
79 fprintf(output, "\n smx must be greater than 0.0");
86 sum2 = 0.0; // moved out of the `if' block for initialization. Lajos, Feb 10, 2005
88 { // when blob, gamma is 0.0, so this step is skept. hstau 6/25/2003
90 for (i = GeoPar.nelem; i < GeoPar.area - 2 * GeoPar.nelem + 1; i +=
94 for (j = i + 1; j < i + GeoPar.nelem - 2; j++)
98 * (recon[j - 1 - GeoPar.nelem]
99 + recon[j - GeoPar.nelem]
100 + recon[j + 1 - GeoPar.nelem]
101 + recon[j - 1] + recon[j + 1]
102 + recon[j - 1 + GeoPar.nelem]
103 + recon[j + GeoPar.nelem]
104 + recon[j + 1 + GeoPar.nelem]);
105 diffnbs = diff1 * MAPEM.sclfactor;
106 ssum2 += diffnbs * diffnbs;
111 sum2 = -0.5 * MAPEM.gamma1 * sum2;
114 totalsum = sum1 + sum2;
116 //------ write results
119 if (MAPEM.gamma1 == 0.0)
121 fprintf(MAPUser1, "\n iter log-likelihood");
125 ///if(MAPEM.gamma1 != 0.0) {
127 "\n iter log-likelihood log-posterior");
131 if (MAPEM.gamma1 == 0.0)
133 fprintf(MAPUser1, "\n %5i %20.9f", iter, sum1);
137 ///if(MAPEM.gamma1 != 0.0) {
138 fprintf(MAPUser1, "\n %5i %20.9f %20.9f", iter, sum1, totalsum);