Commit 7a35098e authored by Johan Richard's avatar Johan Richard
Browse files

Corrected cleanset 4 parameters in o_run_bayes and identations

parent 4ee8a49c
......@@ -79,7 +79,8 @@ double o_run_bayes()
// Initialize Bayesys for MassInf. WL only with a Grid
if( G.nmsgrid != G.nlens )
{NPRINTF(stderr,"Stating G.nmsgrid != G.nlens\n ");
{
NPRINTF(stderr,"Stating G.nmsgrid != G.nlens\n ");
if( I.stat == 8 || I.n_mult != 0 )
{
extern long int narclet;
......@@ -93,7 +94,7 @@ double o_run_bayes()
long int ilens;
double e;
int k;
NPRINTF(stderr,"I.n_mult = %d\n",I.n_mult);
NPRINTF(stderr,"I.n_mult = %d\n",I.n_mult);
int Ndata = narclet;
for ( i = 0; i < I.n_mult; i++ ) //Adding data from multfile if any
Ndata += I.mult[i];
......@@ -105,7 +106,7 @@ double o_run_bayes()
double *zbitstmp0 = (double *)calloc(Ndata, sizeof(double));
zbitstmp1 = alloc_square_double(G.nlens - G.nmsgrid + 2 * I.n_mult, Ndata);
ibitstmp = (int *)calloc(Ndata, sizeof(int));
NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
// Contribution from the grid potentials
for( ilens = 0; ilens < G.nlens - G.nmsgrid; ilens++ )
......@@ -167,31 +168,12 @@ NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
for( i = 0; i < narclet; i++ )
{
grad2 = e_grad2_pot(&arclet[i].C, ilens);
grad2 = e_grad2_pot(&arclet[i].C, ilens);
zbitstmp0[2 * nmult + i] += (grad2.a - grad2.c) * arclet[i].dr;
zbitstmp0[2 * nmult + narclet + i] += 2. * grad2.b * arclet[i].dr;
}
}
// Shift source positions to (F.xmin, F.ymin)
// nimage = 0;
// double xmin = 1000.;
// double ymin = 1000.;
// for( i = 0; i < I.n_mult; i++ )
// for( j = 0; j < I.mult[i]; j++ )
// {
// if( multi[i][j].C.x < xmin ) xmin = multi[i][j].C.x;
// if( multi[i][j].C.y < ymin ) ymin = multi[i][j].C.y;
// }
//
// for( i = 0; i < I.n_mult; i++ )
// for( j = 0; j < I.mult[i]; j++ )
// {
// zbitstmp0[nimage] += xmin;
// zbitstmp0[nmult + nimage] += ymin;
// nimage++;
// }
// Indexes for the source positions
nimage = 0;
for ( i = 0; i < I.n_mult; i++ )
......@@ -247,24 +229,26 @@ NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
Common->Alpha = -0.02 * (G.nlens - G.nmsgrid) - 2 * I.n_mult - 1; // +ve for Poisson, -ve for geometric
if( nDim < G.nlens )
{ NPRINTF(stderr,"Stating nDim<G.nlens\n ");
{
NPRINTF(stderr,"Stating nDim<G.nlens\n ");
// MassInf parameters
Common->Ndim = 1; // Coordinate in the grid
Common->Valency = 1; // On/off switch for MassInf/Grid optim (Default: off)
Common->MassInf = 1;
Common->MassInf = 1;
if( I.n_mult > 0 && G.no_lens > 0 )
{
Common->Valency = 3; // Add 2 valencies for the SL source positions X and Y
Common->Method = -3; // Remove lifeStory2 algorithm for conflict of atom valencies
Common->MassInf = 2; // Positive/Negative prior on the flux
}
Common->ProbON = 0.7;
Common->FluxUnit0 = -10.; // 1 for pmass, 10 for b0
nDim = G.nlens - G.nmsgrid; // To initialize statistics arrays err, tmpCube and avg
Common->ProbON = 0.7;
Common->FluxUnit0 = -10.; // 1 for pmass, 10 for b0
nDim = G.nlens - G.nmsgrid; // To initialize statistics arrays err, tmpCube and avg
}
else
{ NPRINTF(stderr,"Stating nDim = G.nlens\n");
{
NPRINTF(stderr,"Stating nDim = G.nlens\n");
// Normal optimization but with Natoms (see UserBuild)
Common->Ndim = nDim; // Number of parameters of the optimized clumps
Common->Valency = 0; // On/off switch for MassInf/Grid optim (Default: off)
......@@ -273,7 +257,8 @@ NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
}
else
{ NPRINTF(stderr,"Stating G.nmsgrid = G.nlens\n ");
{
NPRINTF(stderr,"Stating G.nmsgrid = G.nlens\n ");
// Number of Atoms = 1 Standard optimization
Common->MinAtoms = 1; // >= 1
......@@ -283,70 +268,69 @@ NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
Common->Ndim = nDim; // Number of parameters of the optimized clumps
Common->Valency = 0; // On/off switch for MassInf/Grid optim (Default: off)
// Entering pixelized source plane reconstruction
if (M.iclean == 4)
{
NPRINTF(stderr,"Entering into cleanset = 4 loop in o_runbayes\n ");
int Ndata = imFrame.nx*imFrame.ny;
NPRINTF(stderr,"NDATA = %d\n",Ndata);
int ns= ps.nx*ps.ny;
// Ndata = 2 * Ndata; //Each data point has two dimensions x and y
double *Data = (double *)malloc((unsigned int)(Ndata * sizeof(double)));
double *Acc = (double *)malloc((unsigned int)(Ndata * sizeof(double)));
double **im;
extern double **imo,**wo;
extern double **fluxmap;
zbitstmp1 = alloc_square_double(ns, Ndata);
ibitstmp = (int *)calloc(Ndata, sizeof(int));
// im = (double **) readimage(&imFrame);
// wo = (double **) readimage(&wFrame);
//Data
for( j = 0; j < Ndata; j++ )
ibitstmp[j] = j; //Filling up ibits matrix
for(i=0;i<imFrame.nx;i++)
{
for(j=0;j<imFrame.ny;j++)
{
Data[i*imFrame.ny+j] = imo[j][i];
if(wFrame.format != 0)
Acc[i*imFrame.ny+j]=sqrt(wo[j][i]);
else
Acc[i*imFrame.ny+j]=20.0;
}
}
// free_square_double(im, imFrame.ny);
for( i = 0; i < ns; i++ )
{
for( j = 0; j < Ndata; j++ )
{ zbitstmp1[i][j]=fluxmap[j][i];
// NPRINTF(stderr,"zbitstmp[%d][%d]=%lf\n",i,j,zbitstmp1[i][j]);
}
}
// Entering pixelized source plane reconstruction
if (M.iclean == 4)
{
NPRINTF(stderr,"Entering into cleanset = 4 loop in o_runbayes\n ");
int Ndata = imFrame.nx*imFrame.ny;
NPRINTF(stderr,"NDATA = %d\n",Ndata);
int ns= ps.nx*ps.ny;
// Ndata = 2 * Ndata; //Each data point has two dimensions x and y
double *Data = (double *)malloc((unsigned int)(Ndata * sizeof(double)));
double *Acc = (double *)malloc((unsigned int)(Ndata * sizeof(double)));
extern double **imo,**wo;
extern double **fluxmap;
zbitstmp1 = alloc_square_double(ns, Ndata); //What do they denote?
ibitstmp = (int *)calloc(Ndata, sizeof(int));
//Data
for( j = 0; j < Ndata; j++ )
ibitstmp[j] = j; //Filling up ibits matrix
for(i=0;i<imFrame.nx;i++)
{
for(j=0;j<imFrame.ny;j++)
{
Data[i*imFrame.ny+j] = imo[j][i];
if(wFrame.format != 0)
Acc[i*imFrame.ny+j]=sqrt(wo[j][i]);
else
Acc[i*imFrame.ny+j]=1. ; //sqrt(fabs(imo[j][i]));
}
}
// free_square_double(im, imFrame.ny);
for( i = 0; i < ns; i++ )
{
for( j = 0; j < Ndata; j++ )
{
zbitstmp1[i][j]=fluxmap[j][i];
}
}
Common->Ndata = Ndata;
Common->Data = Data;
Common->Acc = Acc;
Common->Ndata = Ndata;
Common->Data = Data;
Common->Acc = Acc;
// UserBuild uses mock data
for( i = 0; i < Common->ENSEMBLE; i++ )
Objects[i].Mock = (double*)malloc(Ndata * sizeof(double));
// UserBuild uses mock data
for( i = 0; i < Common->ENSEMBLE; i++ )
Objects[i].Mock = (double*)malloc(Ndata * sizeof(double));
// Number of Atoms free
Common->MinAtoms = 1; // >= 1
Common->MaxAtoms = 0; // >= MinAtoms, or 0 = infinity(not implemented)
Common->Alpha = -0.02 * (ps.nx*ps.ny)- 1;
// MassInf parameters
Common->Ndim = 1; // Coordinate in the grid
Common->Valency = 1; // On/off switch for MassInf/Grid optim (Default: off)
Common->MassInf = 1;
Common->ProbON = 0.7;
Common->FluxUnit0 = -10.; // 1 for pmass, 10 for b0
// Number of Atoms free
Common->MinAtoms = 1; // >= 1
Common->MaxAtoms = 0; // >= MinAtoms, or 0 = infinity(not implemented)
Common->Alpha = -0.02 * (ps.nx*ps.ny)- 1;
// MassInf parameters
Common->Ndim = 1; // Coordinate in the grid
Common->Valency = 1; // On/off switch for MassInf/Grid optim (Default: off)
Common->MassInf = 1;
Common->ProbON = 0.7;
Common->FluxUnit0 = -20.; // 1 for pmass, 10 for b0
}
}
}
}
// Initialize some statistics
Common->UserCommon = (void*)UserCommon;
UserCommon->Nsample = 0;
......@@ -354,8 +338,8 @@ NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
UserCommon->Nchi2 = 0;
UserCommon->sum = 0.;
UserCommon->Ndim = nDim;
;
NPRINTF(stderr,"UserCommon->Ndim = %d\n",UserCommon->Ndim);
NPRINTF(stderr,"UserCommon->Ndim = %d\n",UserCommon->Ndim);
tmpCube = (double *)malloc((unsigned int) nDim * sizeof(double));
UserCommon->err = (double *)malloc((unsigned int) nDim * sizeof(double));
UserCommon->avg = (double *)malloc((unsigned int) nDim * sizeof(double));
......@@ -370,6 +354,7 @@ NPRINTF(stderr,"Filling up of zbits and ibits matrices\n");
optInterrupt = 0;
// Initialize mutex for nchi2 counts
init_mutex();
//Run bayesys
int code = BayeSys3(Common, Objects);
printf("\nBayeSys3 return code %d\n", code);
......@@ -457,13 +442,8 @@ int bayesHeader()
#endif
fprintf( bayes, "#Nsample\n");
//if ( sigposAs.bk != 0 || I.dsigell != -1 ||
// (I.stat == 6 && I.statmode == 1) )
fprintf( bayes, "#ln(Lhood)\n");
//else
// fprintf( bayes, "#Chi2\n");
fprintf( bayes, "#ln(Lhood)\n");
// fprintf( bayes, "#Evidence\n");
// Number of lens parameters in the atom
for ( i = 0; i < G.no_lens; i++ )
......@@ -481,13 +461,13 @@ int bayesHeader()
if( block[i][B0] != 0 )
nDim++;
}
if (M.iclean == 4)
{
for ( i = 0; i < (ps.nx*ps.ny); i++ )
fprintf( bayes, "#%d : %s\n", i, "flux" );
NPRINTF(stderr,"Entered bayesheader() in o_run_bayes\n");
nDim = ps.nx*ps.ny;
}
if (M.iclean == 4)
{
for ( i = 0; i < (ps.nx*ps.ny); i++ )
fprintf( bayes, "#%d : %s\n", i, "flux" );
NPRINTF(stderr,"Entered bayesheader() in o_run_bayes\n");
nDim = ps.nx*ps.ny;
}
// source positions of multiple images with grid optimization
if( G.nmsgrid != G.nlens && nDim < G.nlens - G.nmsgrid )
......
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