Skip to content

GitLab

Commit 86c07ad3 authored by Johan Richard's avatar Johan Richard
Browse files

Add RBF in source reconstructions

parent 79c11ff3
Hide whitespace changes
Inline Side-by-side
include/fonction.h
View file @ 86c07ad3
......@@ -470,7 +470,7 @@ float einasto_alpha( double r, double rs, int n, double rhos, double kappa_crit)
void read_table_einasto();
void precompsource();
void testpixel(char *imname,char *sname,double *sflux);
void testpixel(char *imname,char *sname,double *sflux,double **simu);
//#ifdef __cplusplus
//};
......
src/bayesapp.c
View file @ 86c07ad3
......@@ -930,7 +930,8 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
int k, ipx; // counters
long int i, j;
double lhood0, chi2;
extern double *sflux;
extern double *sflux,**ero;
//.............................................................................
// USER IMPOSES FINAL TEMPERATURE AND TERMINATION CONDITION !
if ( Common->cool >= 1.0 && M.inverse == 3 )
......@@ -1334,13 +1335,13 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
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
}
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");
......@@ -1355,13 +1356,12 @@ if (M.iclean ==4) /// Added to create multiple testpixel
M.itmax = Common->Nsystem;
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);
testpixel(imname,sname,sflux); //calling testpixel
}
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,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));
......
src/o_pixsource.c
View file @ 86c07ad3
......@@ -14,211 +14,163 @@
/****************************************************************/
void precompsource(){
FILE *fp;
struct point Pi, Ps;
const extern struct pot lens[];
const extern struct g_source S;
extern struct g_mode M; // M.iclean
const extern struct g_observ O;
extern struct g_pixel imFrame,ps,PSF;
extern double **fluxmap;
double **tmpimage,xx ;
extern int **nimage;
const extern struct g_frame F;
FILE *fp;
struct point Pi, Ps;
const extern struct pot lens[];
const extern struct g_source S;
extern struct g_mode M; // M.iclean
const extern struct g_observ O;
extern struct g_pixel imFrame,ps,PSF;
extern double **fluxmap;
double **tmpimage,xx ;
extern int **nimage;
const extern struct g_frame F;
// Setting the imframe parameters
/* imFrame.xmin = F.xmin;
imFrame.xmax = F.xmax;
imFrame.ymin = F.ymin;
imFrame.ymax = F.ymax;
imFrame.nx = M.npixel;
double scale = (imFrame.xmax-imFrame.xmin) / (imFrame.nx-1.);
imFrame.ny = round((imFrame.ymax-imFrame.ymin )/ scale + 1.); // to prevent trunc(ny) behavior
imFrame.pixelx = scale;
imFrame.pixely = scale;
imFrame.ech = 10; */
xx = (double)imFrame.ech;
xx = (double)imFrame.ech;
//Setting the ps parameters
double dlsds = dratio(lens[0].z, M.zclean);
fprintf(stderr,"dlsds=%lf\n",dlsds);
if( strcmp(M.centerfile, "") )
{
ps.pixelx = imFrame.pixelx / ps.ech;
ps.pixely = imFrame.pixelx/ ps.ech;
s_sourcebox(&ps, M.centerfile, dlsds, M.zclean);
}
//Setting the ps parameters
double dlsds = dratio(lens[0].z, M.zclean);
if( strcmp(M.centerfile, "") )
{
ps.pixelx = imFrame.pixelx / ps.ech;
ps.pixely = imFrame.pixelx/ ps.ech;
s_sourcebox(&ps, M.centerfile, dlsds, M.zclean);
}
int ni=imFrame.nx*imFrame.ny;
int ns=ps.nx*ps.ny;
double stepx,stepxs,stepy,stepys;
fluxmap=(double **)alloc_square_double(ni,ns);
double sizeratio=imFrame.pixelx/ps.pixelx; // pixel size ratio between image and source plane
double xs,ys;
tmpimage=(double **)alloc_square_double(imFrame.ny,imFrame.nx);
int is,js,i,j,ii,jj;
double z,rbfsig;
stepx=(imFrame.xmax-imFrame.xmin)/(imFrame.nx-1);
stepy=(imFrame.ymax-imFrame.ymin)/(imFrame.ny-1);
stepxs=(ps.xmax-ps.xmin)/(ps.nx-1);
stepys=(ps.ymax-ps.ymin)/(ps.ny-1);
rbfsig = (double)(O.rbfsig/2.355) ; //Setting the sigma of gaussian radial basis functions
for (i=0;i<ni;i++)
for(j=0;j<ns;j++)
fluxmap[i][j]=0.0;
int ni=imFrame.nx*imFrame.ny;
int ns=ps.nx*ps.ny;
double stepx,stepxs,stepy,stepys;
fluxmap=(double **)alloc_square_double(ni,ns);
double sizeratio=imFrame.pixelx/ps.pixelx; // pixel size ratio between image and source plane
double xs,ys;
fprintf(stderr,"ps.xmax=%lf\t ps.ymax=%lf\n",ps.xmax,ps.ymax);
fprintf(stderr,"ps.xmin=%lf\t ps.ymin=%lf\n",ps.xmin,ps.ymin);
fprintf(stderr,"imframe.nx=%d\t imFrame.ny=%d\n",imFrame.nx,imFrame.ny);
fprintf(stderr,"imFrame.pixelx=%lf\t imFrame.pixely=%lf\n",imFrame.pixelx,imFrame.pixely);
fprintf(stderr,"ps.pixelx=%lf\t ps.pixely=%lf\n",ps.pixelx,ps.pixely);
fprintf(stderr,"imFrame.xmax=%lf\t imFrame.ymax=%lf\n",imFrame.xmax,imFrame.ymax);
fprintf(stderr,"imFrame.xmin=%lf\t imFrame.ymin=%lf\n",imFrame.xmin,imFrame.ymin);
fprintf(stderr, "sizeratio= %lf\n", sizeratio);
fprintf(stderr,"ech=%d\n",imFrame.ech);
tmpimage=(double **)alloc_square_double(imFrame.ny,imFrame.nx);
nimage=(int **)alloc_square_int(imFrame.ny,imFrame.nx);
int is,js,i,j,ii,jj;
double z;
stepx=(imFrame.xmax-imFrame.xmin)/(imFrame.nx-1);
stepy=(imFrame.ymax-imFrame.ymin)/(imFrame.ny-1);
stepxs=(ps.xmax-ps.xmin)/(ps.nx-1);
stepys=(ps.ymax-ps.ymin)/(ps.ny-1);
fprintf(stderr,"stepxs=%.8lf\t stepys=%.8lf\n",stepxs,stepys);
fprintf(stderr,"stepx=%.8lf\t stepy=%.8lf\n",stepx,stepy);
for (i=0;i<imFrame.nx;i++)
for(j=0;j<imFrame.ny;j++)
{
tmpimage[j][i]=0.0;
nimage[j][i]=0;
}
for (i=0;i<ni;i++)
for(j=0;j<ns;j++)
fluxmap[i][j]=0.0;
xs= ps.xmin;
// ys=ps.ymin;
// fp = fopen("fluxmap.dat","w");
for(is=0;is<ps.nx;is++)
{ fprintf(stderr, "is=%d\n",is);
ys=ps.ymin;
for(js=0;js<ps.ny;js++)
{ //Pi.x= imFrame.xmin - stepx/2;
Pi.x= imFrame.xmin - stepx/2 * (1.- 1./xx);
for (i=0;i<imFrame.nx;i++)
for(j=0;j<imFrame.ny;j++)
{
tmpimage[j][i]=0.0;
}
for(i=0;i<imFrame.nx;i++)
{
for (ii = 0; ii < imFrame.ech; ii++)
{ // Pi.y=imFrame.ymin - stepy/2;
Pi.y=imFrame.ymin - stepy/2 * (1.- 1./xx);
for(j=0;j<imFrame.ny;j++)
{
for (jj = 0; jj < imFrame.ech; jj++)
{
e_dpl(&Pi, dlsds,S.zs, &Ps);
z= fmax(fabs(Ps.x-xs)/stepxs,fabs(Ps.y-ys)/stepys);
if(z<0.5)
{
nimage[j][i]++;
// fprintf(stderr,"nimage[%d][%d]=%d\n",j,i,nimage[j][i]);
tmpimage[j][i]+=1.0;//(sizeratio*sizeratio);
}
Pi.y +=stepy/xx;
}
}
Pi.x +=stepx/xx;
}
}
/* convolve tmpimage with PSF */
if (O.setseeing)
d_seeing_omp(tmpimage, imFrame.nx, imFrame.ny, imFrame.pixelx);
/* store tmpimage in flux matrix */
for(i=0;i<imFrame.nx;i++)
{
for(j=0;j<imFrame.ny;j++)
{
fluxmap[i*imFrame.ny+j][is*ps.ny+js]=tmpimage[j][i]/xx/xx;
// if(fluxmap[i*imFrame.ny+j][is*ps.ny+js]>1.e-15)
// fprintf(fp,"fluxmap[%d][%d]=%e\t tmpimage[%d][%d]=%e\n",(i*imFrame.ny+j),(is*ps.ny+js),fluxmap[i*imFrame.ny+j][is*ps.ny+js],j,i,tmpimage[j][i]);
}
}
ys+=stepys;
}
xs+=stepxs;
}
// fclose(fp);
xs= ps.xmin;
for(is=0;is<ps.nx;is++)
{
ys=ps.ymin;
for(js=0;js<ps.ny;js++)
{
Pi.x= imFrame.xmin - stepx/2 * (1.- 1./xx);
for (i=0;i<imFrame.nx;i++)
for(j=0;j<imFrame.ny;j++)
tmpimage[j][i]=0.0;
for(i=0;i<imFrame.nx;i++)
{
for (ii = 0; ii < imFrame.ech; ii++)
{
Pi.y=imFrame.ymin - stepy/2 * (1.- 1./xx);
for(j=0;j<imFrame.ny;j++)
{
for (jj = 0; jj < imFrame.ech; jj++)
{
e_dpl(&Pi, dlsds,S.zs, &Ps);
if (O.setrbfsig)
{
z= ((Ps.x-xs)*(Ps.x-xs)/stepxs/stepxs)+((Ps.y-ys)*(Ps.y-ys)/stepys/stepys);
tmpimage[j][i]+=(1.0/(2*3.1415927*rbfsig*rbfsig))*exp(-z/(2.0*rbfsig*rbfsig));
}
else
{
z= fmax(fabs(Ps.x-xs)/stepxs,fabs(Ps.y-ys)/stepys);
if(z<0.5)
tmpimage[j][i]+=1.0;
}
Pi.y +=stepy/xx;
}
}
Pi.x +=stepx/xx;
}
}
/* convolve tmpimage with PSF */
if (O.setseeing)
d_seeing_omp(tmpimage, imFrame.nx, imFrame.ny, imFrame.pixelx);
wri_fits_abs("npixel_abs.fits", nimage, imFrame.nx, imFrame.ny, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, M.ref_ra, M.ref_dec);
free_square_double(tmpimage,imFrame.ny);
/* store tmpimage in flux matrix */
for(i=0;i<imFrame.nx;i++)
{
for(j=0;j<imFrame.ny;j++)
{
if(M.flux!=0) tmpimage[j][i] = tmpimage[j][i]*((double)ps.ech*ps.ech);
fluxmap[i*imFrame.ny+j][is*ps.ny+js]=tmpimage[j][i]/xx/xx;
}
}
ys+=stepys;
}
xs+=stepxs;
}
free_square_double(tmpimage,imFrame.ny);
}
void testpixel(char *imname,char *sname,double *sflux)
void testpixel(char *imname,char *sname,double *sflux,double **simu)
{
const extern struct g_observ O;
extern struct g_mode M;
extern struct g_pixel ps,imFrame;
extern double **fluxmap;
double xmin, ymin, xmax, ymax, dx, dy;
double **zz,**source;
double f, scale,xx;
int nx, ny;
FILE *fp;
int i,j,k;
extern int **nimage;
double ra, dec, width, height;
xmin = imFrame.xmin;
xmax = imFrame.xmax;
ymin = imFrame.ymin;
ymax = imFrame.ymax;
const extern struct g_observ O;
extern struct g_mode M;
extern struct g_pixel ps,imFrame;
extern double **fluxmap;
double xmin, ymin, xmax, ymax, dx, dy;
double **source;
double f, scale,xx;
int nx, ny;
FILE *fp;
int i,j,k;
extern int **nimage;
double ra, dec, width, height;
xmin = imFrame.xmin;
xmax = imFrame.xmax;
ymin = imFrame.ymin;
ymax = imFrame.ymax;
dx = xmax - xmin;
dy = ymax - ymin;
nx = imFrame.nx;// default: assign np to nx and adjust ny
// ny = imFrame.ny;
scale = dx / (nx-1);
ny = round(dy / scale + 1.); // to prevent trunc(ny) behaviour
dx = xmax - xmin;
dy = ymax - ymin;
nx = imFrame.nx;// default: assign np to nx and adjust ny
scale = dx / (nx-1);
ny = round(dy / scale + 1.); // to prevent trunc(ny) behaviour
dx = (nx - 1) * scale;
dy = (ny - 1) * scale;
xmax = xmin + dx;
ymax = ymin + dy;
int ni=imFrame.nx*imFrame.ny;
int ns=ps.nx*ps.ny;
fprintf(stderr, "\timage (%d,%d) s=%.3lf [%.3lf,%.3lf] [%.3lf,%.3lf]\n",nx, ny, scale, xmin, xmax, ymin, ymax);
dx = (nx - 1) * scale;
dy = (ny - 1) * scale;
xmax = xmin + dx;
ymax = ymin + dy;
int ni=imFrame.nx*imFrame.ny;
int ns=ps.nx*ps.ny;
zz = (double **) alloc_square_double(ny, nx);
fprintf(stderr,"nx=%d\t ny=%d\n",nx,ny);
char filename[128];
sprintf(filename, "params_%s.dat", imname);
// fprintf(stderr,"filename = %s\n",filename);
// fp = fopen(filename,"w");
for(i=0;i<imFrame.nx;i++)
{
for(j=0;j<imFrame.ny;j++)
{ zz[j][i]=0.0;
// xx = (double)nimage[j][i];
for(k=0;k<ns;k++)
{ //fprintf(fp,"fluxmap[%d][%d]=%e\t sflux[%d]=%e\t",i*ny+j,k,fluxmap[i*ny+j][k],k,sflux[k]);
zz[j][i]+=fluxmap[i*ny+j][k]*sflux[k];
}
// fprintf(fp,"zz[%d][%d]=%e\t fluxmap[%d][%d]=%e\n",j,i,zz[j][i]);
// if (nimage[j][i]>1)
// zz[j][i]=zz[j][i]/xx;
}
}
// fclose(fp);
wrf_fits_abs(imname, zz, nx, ny, xmin, xmax, ymin, ymax, M.ref_ra, M.ref_dec);
source = (double **) alloc_square_double(ps.ny, ps.nx);
// Also write the corresponding source file using the ps parameters
for(i=0;i<ps.nx;i++)
{
for(j=0;j<ps.ny;j++)
source[j][i]=sflux[i*ps.ny+j];
}
// Nprintf(stderr,"source[%d][%d]=%e\n",9,19,source[9][19]);
wrf_fits_abs(sname, source, ps.nx, ps.ny, ps.xmin, ps.xmax, ps.ymin, ps.ymax, M.ref_ra, M.ref_dec);
free_square_double(zz, ny);
free_square_double(source, ps.ny);
NPRINTF(stderr, "\timage (%d,%d) s=%.3lf [%.3lf,%.3lf] [%.3lf,%.3lf]\n",nx, ny, scale, xmin, xmax, ymin, ymax);
source = (double **) alloc_square_double(ps.ny, ps.nx);
for(i=0;i<ps.nx;i++)
for(j=0;j<ps.ny;j++)
source[j][i]=sflux[i*ps.ny+j];
if (O.setrbfsig)
d_rbf(source, ps.nx, ps.ny, 1.);
for(i=0;i<imFrame.nx;i++)
for(j=0;j<imFrame.ny;j++)
{
simu[j][i]=0.0;
for(k=0;k<ns;k++)
simu[j][i]+=fluxmap[i*ny+j][k]*sflux[k];
}
wrf_fits_abs(imname, simu, nx, ny, xmin, xmax, ymin, ymax, M.ref_ra, M.ref_dec);
wrf_fits_abs(sname, source, ps.nx, ps.ny, ps.xmin, ps.xmax, ps.ymin, ps.ymax, ra, dec);
free_square_double(source, ps.ny);
}
src/readBayesModels.c
View file @ 86c07ad3
......@@ -199,7 +199,7 @@ void setBayesModel( long int iVal, long int nVal, double **array)
extern int cblock[NPAMAX];
extern int vfblock[NPAMAX];
extern struct sigposStr sigposAs;
extern double *np_b0,*sflux;
extern double *np_b0,*sflux,**ero;
extern long int narclet;
extern struct g_pixel ps;
......@@ -302,11 +302,11 @@ void setBayesModel( long int iVal, long int nVal, double **array)
// cleanset = 4 case
if (M.iclean == 4)
{
int ns = ps.nx*ps.ny;
imethod = -2;
for ( ilens = 0; ilens < ns; ilens++ )
sflux[ilens] = method(nVal, iParam++, array);
testpixel("testpixel_method-2.fits","sflux_method-2.fits",sflux);
int ns = ps.nx*ps.ny;
imethod = -2;
for ( ilens = 0; ilens < ns; ilens++ )
sflux[ilens] = method(nVal, iParam++, array);
testpixel("testpixel_method-2.fits","sflux_method-2.fits",sflux,ero);
}
// Set the cosmological parameters
......
Markdown is supported
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