Commit 82245305 authored by Johan Richard's avatar Johan Richard
Browse files

Corrected output and identation in bayesapp/cleanset 4

parent 7a35098e
......@@ -624,8 +624,6 @@ int UserFoot(
long int ilens;
int i;
int nDim = Common->Ndim;
FILE *fp;
//fp= fopen("userfoot.dat","a");
if (M.iclean == 4)
{
ilens = (int)(ps.nx* ps.ny * Cube[nDim-1]);
......@@ -633,7 +631,6 @@ int UserFoot(
nbits[0] = Common->Ndata; nbits[1] = 0; nbits[2] = 0;
memcpy(ibits, ibitstmp, (unsigned int)(Common->Ndata * sizeof(int)));;
memcpy(zbits, zbitstmp1[ilens], (unsigned int)(Common->Ndata * sizeof(double)));
//fprintf(fp,"%d \t %lf \t %lf\n",Common->Nsystem,Cube[0],Cube[1]);
}
else
{ //fprintf(stderr,"Entered else in userfoot\n");
......@@ -671,7 +668,6 @@ int UserFoot(
}
}
increase_nchi2();
//fclose(fp);
// fprintf(stderr,"cube=%lf\t zbits=%lf \t nbits=%d\n",Cube[ilens],zbits[ilens],nbits[0]);
return (Cube && Common && ibits && zbits && nbits) ? 1 : 1;
}
......@@ -882,7 +878,7 @@ static void addStat(double val, UserCommonStr *UserCommon, int param, long int s
UserCommon->avg[param] += val;
mean = UserCommon->avg[param] / samp;
var1 = UserCommon->err[param] / samp - mean * mean;
var = abs(var1); //TV
var = fabs(var1); //TV
UserCommon->sum += var / mean / mean;
}
......@@ -926,12 +922,11 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
UserCommonStr* UserCommon = (UserCommonStr*)Common->UserCommon;
extern double *tmpCube;
char bayesName[30],burninName[30];
FILE *bayes;
FILE *bayes,*fp1;
int k, ipx; // counters
long int i, j;
double lhood0, chi2;
extern double *sflux,**ero;
//.............................................................................
// USER IMPOSES FINAL TEMPERATURE AND TERMINATION CONDITION !
if ( Common->cool >= 1.0 && M.inverse == 3 )
......@@ -951,12 +946,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
#endif
sprintf(bayesName, "%s.dat", bayesName);
sprintf(burninName, "%s.dat", burninName);
// char imname[128],sname[128];
if(M.iclean==4)
{
}
if ( Common->cool < 1. )
bayes = fopen( burninName, "a");
else
......@@ -979,78 +969,73 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
// Initialize arrays
double *Mock = (double *) calloc(Ndata, sizeof(double));
if (M.iclean == 4 )
{ // fprintf(stderr, "Entering User Monitor\n");
extern double **ero; // To call testpixel at each stage
//
int ns = ps.nx*ps.ny ;
memset(sflux, 0, (unsigned long int)( ns* sizeof(double)));
for( i = 0; i < Objects[k].Natoms; i++ )
if (M.iclean == 4 )
{
int ns = ps.nx*ps.ny ;
memset(sflux, 0, (unsigned long int)( ns* sizeof(double)));
for( i = 0; i < Objects[k].Natoms; i++ )
{
ilens = (int)(ns * Cubes[i][Common->Ndim - 1]);
sflux[ilens] += Cubes[i][Common->Ndim];
for( j = 0; j < Ndata; j++ )
Mock[j] += Cubes[i][Common->Ndim] * zbitstmp1[ilens][j];
for( j = 0; j < Ndata; j++ )
Mock[j] += Cubes[i][Common->Ndim] * zbitstmp1[ilens][j];
}
}
else
{
memset(np_b0, 0, (unsigned long int)((G.nlens - G.nmsgrid) * sizeof(double)));
// MassInf optimization (commented: 1 optimized pot.)
for( i = 0; i < Objects[k].Natoms; i++ )
{
// Not implemented in UserBuild()
// for( j = 0; j < Common->Ndim - 2; j ++ )
// tmpCube[j] += Cubes[i][j]; // TODO: rescale between 0..1 with Unif prior
}
else
{
memset(np_b0, 0, (unsigned long int)((G.nlens - G.nmsgrid) * sizeof(double)));
// MassInf optimization (commented: 1 optimized pot.)
for( i = 0; i < Objects[k].Natoms; i++ )
{
// Not implemented in UserBuild()
// for( j = 0; j < Common->Ndim - 2; j ++ )
// tmpCube[j] += Cubes[i][j]; // TODO: rescale between 0..1 with Unif prior
// Cubes[][Common->Ndim - 1] is coordinate
ilens = (int)((G.nlens - G.nmsgrid + I.n_mult) * Cubes[i][Common->Ndim - 1]);
if( ilens < G.nlens - G.nmsgrid )
{
np_b0[ilens] += Cubes[i][Common->Ndim];
for( j = 0; j < Ndata; j++ )
Mock[j] += Cubes[i][Common->Ndim] * zbitstmp1[ilens][j];
}
else
{
int isrc = ilens - (G.nlens - G.nmsgrid);
int nimage = 0;
for( j = 0; j < isrc; j ++ )
nimage += I.mult[j];
// Cubes[][Common->Ndim - 1] is coordinate
ilens = (int)((G.nlens - G.nmsgrid + I.n_mult) * Cubes[i][Common->Ndim - 1]);
if( ilens < G.nlens - G.nmsgrid )
{
np_b0[ilens] += Cubes[i][Common->Ndim];
for( j = 0; j < Ndata; j++ )
Mock[j] += Cubes[i][Common->Ndim] * zbitstmp1[ilens][j];
}
else
{
int isrc = ilens - (G.nlens - G.nmsgrid);
int nimage = 0;
for( j = 0; j < isrc; j ++ )
nimage += I.mult[j];
int nmult = nimage;
for( j = isrc; j < I.n_mult; j++ )
nmult += I.mult[j];
int nmult = nimage;
for( j = isrc; j < I.n_mult; j++ )
nmult += I.mult[j];
for( j = 0; j < I.mult[isrc]; j++ )
{
Mock[nimage] += Cubes[i][Common->Ndim + 1];
Mock[nimage + nmult] += Cubes[i][Common->Ndim + 2];
nimage++;
}
}
}
}
// Compute Likelihood from SL
chi2 = 0.; lhood0 = 0.;
for( j = 0; j < I.mult[isrc]; j++ )
{
Mock[nimage] += Cubes[i][Common->Ndim + 1];
Mock[nimage + nmult] += Cubes[i][Common->Ndim + 2];
nimage++;
}
}
}
}
// Compute Likelihood from SL
chi2 = 0.; lhood0 = 0.;
// Accumulate logLikelihood from grid
for( j = 0; j < Ndata; j++ )
{
chi2 += (Mock[j] - Data[j]) * Acc[j] * Acc[j] * (Mock[j] - Data[j]);
lhood0 += log(Sqrt2Pi / Acc[j]);
// fprintf(stderr,"chi2=%e \t lhood0=%e\n",chi2,lhood0);
// fprintf(stderr,"Mock[%d]=%e \t Data[%d]=%e\n",j,Mock[j],j,Data[j]);
}
// Accumulate logLikelihood from grid
for( j = 0; j < Ndata; j++ )
{
chi2 += (Mock[j] - Data[j]) * Acc[j] * Acc[j] * (Mock[j] - Data[j]);
lhood0 += log(Sqrt2Pi / Acc[j]);
}
// Compute Likelihood (not to include otherwise there is no match
// when reading with readBayesModel.c)
//int error = o_chi_lhood0(&chi2, &lhood0, np_b0); // TODO: catch error
free(Mock);
}
// Compute Likelihood (not to include otherwise there is no match
// when reading with readBayesModel.c)
//int error = o_chi_lhood0(&chi2, &lhood0, np_b0); // TODO: catch error
free(Mock);
}
else if( G.nmsgrid != G.nlens )
{
// CONDITION NOT USED BECAUSE EVERYTHING IS DONE WITH MASSINF NOW
......@@ -1126,25 +1111,26 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
addStat(val, UserCommon, iparam++, samp);
}
if (M.iclean == 4)
{ nzeros= 0;
val = 0.;
int ns = ps.nx*ps.ny;
for ( i = 0; i < ns; i++ )
{
val = sflux[i];
if ( fabs(val) > 1E-6 )
{
if( nzeros > 0 )
fprintf(bayes, "%dx0.0 %lf ", nzeros, val);
else
fprintf(bayes, "%lf ", val);
nzeros = 0;
}
else
nzeros++;
{
nzeros= 0;
val = 0.;
int ns = ps.nx*ps.ny;
for ( i = 0; i < ns; i++ )
{
val = sflux[i];
if ( fabs(val) > 1E-6 )
{
if( nzeros > 0 )
fprintf(bayes, "%dx0.0 %lf ", nzeros, val);
else
fprintf(bayes, "%lf ", val);
nzeros = 0;
}
else
nzeros++;
addStat(val, UserCommon, iparam++, samp);
}
addStat(val, UserCommon, iparam++, samp);
}
}
if ( nzeros > 0 )
fprintf(bayes, "%dx0.0 ", nzeros);
......@@ -1324,26 +1310,17 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
}
fclose(bayes);
//.............................................................................
// Print the status line of STDOUT
int Nrem;
if( Common->cool < 1. )
{ // fprintf(stderr,"chi2=%e\n",chi2);
{
double chi2rate = difftime(end, start);
chi2rate = chi2rate > 0 ? chi2rate : -1;
// fprintf(stderr,"chi2rate=%lf\n",chi2rate);
chi2rate = nchi2 / chi2rate;
Nrem = Common->Nsystem; // # to go, as (2 * Nsample++) catches up with Nsystem++
if (M.iclean ==4) /// Added to create multiple testpixel
{
char imname[128],sname[128];
sprintf(imname, "testpixel_burn%ld.fits", Nrem);
sprintf(sname, "sflux_burn%ld.fits", Nrem);
testpixel(imname,sname,sflux,ero); //calling testpixel
}
printf(" \r");
fprintf(stdout, // (stderr flushes output immediately)
"Burn-in : %10.6lf %5d %12.3lf %10.3lf %d/%.0lf %.0lfchi2/s\r",
......@@ -1357,15 +1334,16 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
Nrem = M.itmax - UserCommon->Nsample; // or any other setting
if (M.iclean ==4) /// Added to create multiple testpixel
{
sprintf(imname, "testpixel_samp%ld.fits", UserCommon->Nsample);
sprintf(sname, "sflux_samp%ld.fits", UserCommon->Nsample);
{
sprintf(imname, "testpixel_samp%ld.fits", UserCommon->Nsample);
sprintf(sname, "sflux_samp%ld.fits", UserCommon->Nsample);
testpixel(imname,sname,sflux,ero); //calling testpixel
}
}
UserCommon->sum /= Common->ENSEMBLE * Common->Ndim;
fprintf(stdout, "Sampling : %10.6lf %7ld/%d %12.3lf %10.3lf[CTRL-C to interrupt]\r",
Common->cool, UserCommon->Nsample + 1, M.itmax + 1, chi2, sqrt(UserCommon->sum));
}
}
fflush(stdout);
minusone = 0;
nchi2 = 0.;
......
Supports Markdown
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