Commit 82a148c6 authored by Johan Richard's avatar Johan Richard
Browse files

Source reconstruction Massinf

parent 71a3c881
......@@ -468,7 +468,7 @@ float einasto_alpha( double r, double rs, int n, double rhos, double kappa_crit)
void read_table_einasto();
void precompsource();
void testpixel();
void testpixel(char *imname,char *sname,double *sflux);
//#ifdef __cplusplus
//};
......
......@@ -318,11 +318,10 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
// MassInf optimization
if( Common->Valency )
{
{ // fprintf(stderr,"Entering UserBuild\n");
// Weak-Lensing is done in UserFoot() under the linear approximation
*Lhood = Object->Lhood;
return 1;
// Compute the likelihood for the SL
extern struct g_grille G;
double *np_b0 = (double *)calloc(G.nlens - G.nmsgrid, sizeof(double));
......@@ -339,8 +338,8 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
for( i = 0; i < I.n_mult; i++ )
{
vB[i].x = vB[i].y = 0.;
for( j = 0; j < I.mult[i]; j++ )
{
for( j = 0; j < I.mult[i]; j++ ) /* Not sure what is happening */
{ /* in these loops */
vB[i].x += multi[i][j].C.x;
vB[i].y += multi[i][j].C.y;
}
......@@ -401,7 +400,7 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
//----------------------------------------------------------------------
// Case with WL only but with moving clumps (Common->Valency == 0)
// Copy lens[G.nmsgrid], Natoms times
// Copy lens[G.nmsgrid], Natoms times ---- Which case is this ??
//----------------------------------------------------------------------
extern struct g_image I;
extern struct g_grille G;
......@@ -428,7 +427,6 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
int ipx, i;
int valid = 1;
*Lhood = 0.;
// Create Mock data
double *Mock = (double *) calloc(Ndata, sizeof(double));
......@@ -457,8 +455,8 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
Mock[narclet + k] += 2. * grad2.b * arclet[k].dr;
}
}
// Accumulate logLikelihood from grid
// Accumulate logLikelihood from grid ?? What does this do?
double Lnorm = 0.0; // Likelihood normalisation
double C = 0.0; // chisquared
for( k = 0; k < Ndata; k++ )
......@@ -495,7 +493,7 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
// Check for single atom
if ( Natoms != 1 )
return valid;
if ( M.restart != 0 )
{
// Read Lhood from restart file if it exists
......@@ -512,11 +510,11 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
return valid;
}
}
// Copy by value of the cube (cube is modified in rescaleCube_1Atom)
for( j = 0; j < Common->Ndim; j++ )
tmpCube[j] = Cubes[0][j];
/* Set the clumps parameters from the cube and eventually reject the whole object*/
valid = rescaleCube_1Atom(tmpCube, Common->Ndim);
......@@ -525,7 +523,7 @@ static int UserBuild( // O +ve = OK, 0 = DO NOT USE, -ve = error
/* Compute the likelihood for this object */
if ( valid )
{
{
UserCommon->Nchi2++;
*Lhood = o_lhood(&error);
......@@ -621,20 +619,33 @@ int UserFoot(
extern double **zbitstmp1;
extern struct g_grille G;
extern struct g_image I;
extern struct g_pixel ps;
extern struct g_mode M;
long int ilens;
int i;
int nDim = Common->Ndim;
ilens = (int)((G.nlens - G.nmsgrid + I.n_mult) * Cube[nDim-1]);
if( ilens < G.nlens - G.nmsgrid )
{
FILE *fp;
fp= fopen("userfoot.dat","a");
if (M.iclean == 4)
{
ilens = (int)(ps.nx* ps.ny * Cube[nDim-1]);
// fprintf(stderr,"Entering userfoot for cleanset = 4\n");
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");
ilens = (int)((G.nlens - G.nmsgrid + I.n_mult) * Cube[nDim-1]);
if( ilens < G.nlens - G.nmsgrid ) //this means no multfile involved
{
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)));
}
else
{
{ //this means multfile involved
int j;
int isrc = ilens - (G.nlens - G.nmsgrid);
......@@ -656,10 +667,12 @@ int UserFoot(
zbits[j+I.mult[isrc]] = 1.;
nimage++;
}
}
}
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;
}
......@@ -863,9 +876,8 @@ static double getParVal(int i, int ipx)
}
static void addStat(double val, UserCommonStr *UserCommon, int param, long int samp)
{
{ //fprintf(stderr,"entered addstat \n");
double var, var1, mean;
UserCommon->err[param] += val * val;
UserCommon->avg[param] += val;
mean = UserCommon->avg[param] / samp;
......@@ -900,27 +912,25 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
CommonStr* Common, // I Full general information
ObjectStr* Objects) // I Full ensemble of sample objects [ENSEMBLE]
{
// Track NaN errors
if ( Common->cool != Common->cool )
{
fprintf(stderr, "ERROR: Common->cool=NaN found in bayesapp.c:UserMonitor().\n");
exit(2);
}
// Variables used everywhere
const double Sqrt2Pi = 2.50662827463100050240;
extern struct g_mode M;
extern struct g_grille G;
extern struct g_pixel ps;
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;
//.............................................................................
// USER IMPOSES FINAL TEMPERATURE AND TERMINATION CONDITION !
if ( Common->cool >= 1.0 && M.inverse == 3 )
......@@ -940,7 +950,14 @@ 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)
{
fp1 = fopen("lenstoolpar.dat","a");
}
// fprintf(fp1," #Iter \t #object \t #natoms\t #nzeros\n");
if ( Common->cool < 1. )
bayes = fopen( burninName, "a");
else
......@@ -953,7 +970,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
extern double *np_b0;
if( Common->Valency )
{
{
extern struct g_image I;
extern double **zbitstmp1;
int Ndata = Common->Ndata;
......@@ -963,8 +980,31 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
// Initialize arrays
double *Mock = (double *) calloc(Ndata, sizeof(double));
memset(np_b0, 0, (unsigned long int)((G.nlens - G.nmsgrid) * 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)));
// sprintf(imname, "testpixel_iter%d_obj%d.fits",Common->Nsystem,k);
// sprintf(sname, "sflux_iter%d_obj%d.fits", Common->Nsystem,k);
fprintf(fp1,"%d %d %d\n", Common->Nsystem,k,Objects[k].Natoms);
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];
}
// testpixel(imname,sname,sflux); //calling testpixel
}
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++ )
{
......@@ -977,7 +1017,6 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
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];
}
......@@ -1000,7 +1039,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
}
}
}
}
// Compute Likelihood from SL
chi2 = 0.; lhood0 = 0.;
......@@ -1008,7 +1047,9 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
for( j = 0; j < Ndata; j++ )
{
chi2 += (Mock[j] - Data[j]) * Acc[j] * Acc[j] * (Mock[j] - Data[j]);
lhood0 += log(Sqrt2Pi / Acc[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]);
}
// Compute Likelihood (not to include otherwise there is no match
......@@ -1017,7 +1058,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
free(Mock);
}
else if( G.nmsgrid != G.nlens )
{
{
// CONDITION NOT USED BECAUSE EVERYTHING IS DONE WITH MASSINF NOW
long int ilens;
......@@ -1031,7 +1072,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
int error = o_chi_lhood0(&chi2, &lhood0, np_b0); // TODO: catch error
}
else
{
{
// Reset the tmpCube array
memset(tmpCube, 0, (unsigned long int)(UserCommon->Ndim * sizeof(double)));
......@@ -1053,17 +1094,17 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
if( error )
continue;
}
// fprintf(stderr,"UserCommon->sum_2= %lf\n",UserCommon->sum);
// Appended till here
double o_lhood_rez = (-0.5*(chi2 + lhood0));
fprintf( bayes, "%ld %lf ", UserCommon->Nsample + 1, o_lhood_rez);
// SAVE: save lens parameters
int iparam = 0;
long int samp = UserCommon->Nsample * Common->ENSEMBLE + k + 1;
extern int block[][NPAMAX];
double val = 0.;
for ( i = 0; i < G.no_lens; i++ )
for ( i = 0; i < G.no_lens; i++ ) // doesnt go in this one for WL grid
for ( ipx = CX; ipx <= PMASS; ipx++ )
if ( block[i][ipx] != 0 )
{
......@@ -1074,11 +1115,31 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
int nzeros= 0;
for ( i = 0; i < G.nlens - G.nmsgrid; i++ )
{
{// goes in this one for WL grid
val = np_b0[i];
if ( fabs(val) > 1E-6 )
{// goes in this one for WL grid
if( nzeros > 0 )
fprintf(bayes, "%dx0.0 %lf ", nzeros, val);
else
fprintf(bayes, "%lf ", val);
nzeros = 0;
}
else
nzeros++;
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
......@@ -1089,14 +1150,14 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
nzeros++;
addStat(val, UserCommon, iparam++, samp);
}
}
if ( nzeros > 0 )
fprintf(bayes, "%dx0.0 ", nzeros);
// SAVE: source positions for grid model (for DEBUG purposes for now)
extern struct g_image I;
if ( Common->Valency && I.n_mult != 0)
{
if ( Common->Valency && I.n_mult != 0) //doesnt go here in WL grid
{
double *Ps = (double *)calloc(I.n_mult * 2, sizeof(double));
long int ilens;
for ( i = 0; i < Objects[k].Natoms; i++ )
......@@ -1180,7 +1241,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
extern int vfblock[NPAMAX];
for ( ipx = VFCX; ipx <= VFSIGMA; ipx++ )
if ( vfblock[ipx] != 0 )
{
{
val = getParVal(0, ipx);
if( ipx == VFTHETA ) val *= RTD;
if( ipx == VFI ) val *= RTD;
......@@ -1242,7 +1303,7 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
// SAVE: the nuisance parameters
extern struct sigposStr sigposAs;
if ( sigposAs.bk != 0 )
{
{
for ( i = 0; i < I.n_mult; i++ )
for ( j = 0; j < I.mult[i]; j++ )
{
......@@ -1269,33 +1330,52 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
}
fclose(bayes);
fclose(fp1);
//.............................................................................
// 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); //calling testpixel
}
printf(" \r");
fprintf(stdout, // (stderr flushes output immediately)
"Burn-in : %10.6lf %5d %12.3lf %10.3lf %d/%.0lf %.0lfchi2/s\r",
Common->cool, Nrem, chi2, Common->Evidence, minusone, nchi2, chi2rate);
}
else
{
{ char imname[128],sname[128];
if ( M.itmax == 0 )
M.itmax = Common->Nsystem;
Nrem = M.itmax - UserCommon->Nsample; // or any other setting
// fprintf(fp1,"nrem in samp=%d \n",Common->Nsystem);
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);
testpixel(imname,sname,sflux); //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));
}
fclose(fp1);
fflush(stdout);
minusone = 0;
nchi2 = 0.;
......
......@@ -119,10 +119,10 @@ void e_pixel(int np, char *iname, char *sname, struct galaxie *source)
int j, k, l;
for ( j = 0; j < ps.ny; j++ )
{
A.y = ps.ymin + ((double)j + 1) * ps.pixely;
A.y = ps.ymin + ((double)j ) * ps.pixely;
for ( k = 0; k < ps.nx; k++ )
{
A.x = ps.xmin + ((double)k + 1) * ps.pixelx;
A.x = ps.xmin + ((double)k ) * ps.pixelx;
zz[j][k] = 0.;
for( l = 0; l < S.ns; l++ )
zz[j][k] += d_profil(A.x, A.y, &source[l]);
......@@ -240,7 +240,7 @@ void pixelsource(int np, char *iname)
dy = ymax - ymin;
nx = np; // default: assign np to nx and adjust ny
scale = dx / (nx-1);
ny = dy / scale + 1; // warning: trunc(ny) behavior
ny = round(dy / scale + 1.); // to prevent trunc(ny) behavior
// f = dx / scale;
// nx = (int) f; // BUG in C
// f = dy / scale;
......
......@@ -255,7 +255,7 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
char str[IDSIZE], *pch; // get image identifier
// Compute image
extern struct g_pixel imFrame;
extern struct g_pixel imFrame,ps;
extern struct galaxie source[NFMAX];
const extern struct g_observ O;
......@@ -272,7 +272,7 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
exit(EXIT_FAILURE);
}
testpixel(ero,sflux);
//testpixel(ero,sflux);
// o_pixel(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, map_axx, map_ayy);
//wrf_fits_abs("simu.fits", ero, imFrame.nx, imFrame.ny, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, M.ref_ra, M.ref_dec);
......
......@@ -255,7 +255,7 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
char str[IDSIZE], *pch; // get image identifier
// Compute image
extern struct g_pixel imFrame;
extern struct g_pixel imFrame,ps;
extern struct galaxie source[NFMAX];
const extern struct g_observ O;
......@@ -272,7 +272,7 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
exit(EXIT_FAILURE);
}
testpixel(ero,sflux);
//testpixel(ero,sflux);
// o_pixel(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, map_axx, map_ayy);
//wrf_fits_abs("simu.fits", ero, imFrame.nx, imFrame.ny, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, M.ref_ra, M.ref_dec);
......
......@@ -37,6 +37,7 @@ double drima; /* distance ratio between the image to inverse and the main lens
double **sm_p; /* optimization spline map : initial save */
double *sflux;
double **fluxmap;
int **nimage;
int nshmap;
int optim_z;
int block_s[NLMAX][NPAMAX];
......@@ -84,7 +85,7 @@ void o_global()
/*-------------------------------------------------------*/
/* optimisation */
NPRINTF(stderr, "OPTIMISATION");
if ( G.nmsgrid != G.nlens )
if (( G.nmsgrid != G.nlens ) || (M.iclean == 4))
NPRINTF(stderr, "(NON PARAM)");
NPRINTF(stderr, " :\n");
......@@ -259,7 +260,7 @@ void o_global_free()
if (M.iclean != 0)
{
extern struct g_source S;
const extern struct g_pixel ps;
extern struct g_pixel ps;
extern struct g_pixel imFrame;
extern struct g_pixel wFrame;
extern struct g_cube cubeFrame;
......@@ -279,8 +280,10 @@ void o_global_free()
}
else if (M.iclean == 4 )
{
free_square_double(ero, imFrame.ny);
free(sflux);
// free_square_double(ero, imFrame.ny);
free_square_double(fluxmap, imFrame.ny*imFrame.nx);
free_square_int(nimage,imFrame.ny);
}
else
{
......@@ -363,7 +366,7 @@ void readConstraints()
extern struct cline cl[];
extern struct sigposStr sigposAs;
extern struct galaxie multi[NFMAX][NIMAX];
const extern struct g_pixel ps;
extern struct g_pixel ps;
extern struct galaxie arclet[];
extern long int narclet;
......@@ -411,11 +414,11 @@ void readConstraints()
//check for iclean = 4
if (M.iclean==4)
{//imtosou(M.zclean, ps.pixfile);
{
precompsource();
fprintf(stderr,"okprecomp\n");
fprintf(stderr,"precomp generated\n");
}
......@@ -686,7 +689,7 @@ void readConstraints()
}
// Prepare Accelerated potentials
if ( G.nmsgrid != G.nlens )
if (( G.nmsgrid != G.nlens ) || (M.iclean == 4))
{
prep_non_param();
......@@ -820,7 +823,7 @@ int getNParameters()
extern int block[][NPAMAX];
extern int cblock[], sblock[NFMAX][NPAMAX];
extern struct sigposStr sigposAs;
const extern struct g_pixel ps;
extern struct g_pixel ps;
long int parameters, i, j;
int ipx;
......@@ -828,18 +831,7 @@ int getNParameters()
//check M.iclean =4 case first
if (M.iclean == 4)
{ //Initializing sflux
int ni= ps.nx*ps.ny;
sflux = (double*)alloc_vector_double(ni);
for (i=0; i<ni;i++)
{parameters++;
if (fabs(i-(0.5*ps.nx*ps.ny))< 10)
{
sflux[i]=1.;
}
else
sflux[i]=0.;
}
parameters = ps.nx*ps.ny;
}
......@@ -914,6 +906,8 @@ void prep_non_param()
extern long int narclet;
extern struct galaxie arclet[];
extern struct pot lmax[];
extern struct g_pixel ps,imFrame;
extern struct g_mode M;
struct galaxie *image;
long int i, j, k, l, n, nimages;
......@@ -924,55 +918,73 @@ void prep_non_param()
long int gcount, lcount; // global and local counters
struct ellipse ampli;
// Number of clumps
n = G.nlens - G.nmsgrid;
if (M.iclean == 4 )
{
fprintf(stderr,"entered cleanset = 4 loop in prep_non_param\n");
n = ps.nx*ps.ny;
sflux = (double *) calloc((unsigned) n, sizeof(double));
double **source;
ps.header= imFrame.header;
ps.format = imFrame.format; //Editing the format since its default value is 0
// NPRINTF(stderr,"ps->format = %d",ps.format);
// source = (double **) readimage(&ps);
for(i=0;i<ps.nx;i++) //Initializing sflux
{
for(j=0;j<ps.ny;j++)
sflux[i*ps.ny+j] = 0.;
}
}
else
{
// Number of clumps
n = G.nlens - G.nmsgrid;
// Create a global array for the lens.b0
np_b0 = (double *) calloc((unsigned) n, sizeof(double));