Commit 43ef25ff authored by Johan Richard's avatar Johan Richard
Browse files

Added chi2 computation for pixsource

parent 2eff55a3
......@@ -168,6 +168,8 @@ void follow(struct point A, struct point B, struct point O, double dl0s, do
struct ellipse formeli(double a, double b, double c);
void frprmn( double *p, int n, double ftol, int *iter, double *fret, double (*func)(double*), void (*dfunc)(double*, double*) );
void fr_sq_point(struct point **square, int nbr_lin);
void free_cubic_double(double ***cube,int nbr_lin, int nbr_col);
void free_cubic_int(int ***cube,int nbr_lin, int nbr_col);
void f_shape2(long int *istart, struct galaxie *liste, char *name);
void f_shape3(long int *istart, struct galaxie *liste, char *name);
void f_shape4(long int *istart, struct galaxie *liste, char *name);
......
......@@ -242,6 +242,73 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
}
}
else if (M.iclean == 4)
{ fprintf(stderr,"Enter o_chi");
// Add up barycenter position to source positions
extern struct g_source S;
extern struct galaxie source[NFMAX];
extern struct galaxie multi[NFMAX][NIMAX];
const extern struct g_image I;
struct point Bs[NFMAX]; // list of source barycenters
struct point Ps[NIMAX];
char str[IDSIZE], *pch; // get image identifier
// Compute image
extern struct g_pixel imFrame;
extern struct galaxie source[NFMAX];
const extern struct g_observ O;
double dx = imFrame.xmax - imFrame.xmin;
double dy = imFrame.ymax - imFrame.ymin;
extern double *sflux;
int verbose_save = M.verbose;
M.verbose = 0;
//2Eric: You assume that imFrame.pixelx == imFrame.pixely ?
if (fabs(( imFrame.pixelx - imFrame.pixely)/imFrame.pixelx) > 1e-4)
{
fprintf(stderr, "Error | imFrame.pixelx(%f) != imFrame.pixely(%f) but we assume that they are equal\n", imFrame.pixelx, imFrame.pixely);
exit(EXIT_FAILURE);
}
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);
// if (O.setseeing)
// d_seeing_omp(ero, imFrame.nx, imFrame.ny, imFrame.pixelx);
//if (O.bruit)
// d_bruiter_omp(ero, imFrame.nx, imFrame.ny);
M.verbose = verbose_save;
if (wo != NULL)
{
for (i = 0; i < imFrame.ny; i++ )
for (j = 0; j < imFrame.nx; j++ )
{
ero[i][j] -= imo[i][j];
chi_im += ero[i][j] * ero[i][j] * wo[i][j];
}
lhood0 = lhood_wo; // read from lhood_wo variable pre-computed in o_global()
}
else
{
for (i = 0; i < imFrame.ny; i++ )
for (j = 0; j < imFrame.nx; j++ )
{
ero[i][j] -= imo[i][j];
chi_im += ero[i][j] * ero[i][j] / fabs(imo[i][j] + 1.);
lhood0 += log( 2.*M_PI*fabs(imo[i][j] + 1.) );
}
}
}
else
/* chi_im is the error associated to the transformation image -> source*/
......@@ -865,9 +932,9 @@ static double chi2SglImage( struct galaxie *pima, struct point *ps )
if ( nimages > 1 )
{
if ( fabs(I.forme) == 10 )
if ( abs(I.forme) == 10 )
I2x = I2y = pima->E.a * pima->E.a;
else if ( fabs(I.forme) == 11 )
else if ( abs(I.forme) == 11 )
{
I2x = pima->E.a * cos(pima->E.theta);
I2y = pima->E.b * cos(pima->E.theta);
......
......@@ -242,6 +242,74 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
}
}
else if (M.iclean == 4)
{ fprintf(stderr,"Enter o_chi\n");
// Add up barycenter position to source positions
extern struct g_source S;
extern struct galaxie source[NFMAX];
extern struct galaxie multi[NFMAX][NIMAX];
const extern struct g_image I;
struct point Bs[NFMAX]; // list of source barycenters
struct point Ps[NIMAX];
char str[IDSIZE], *pch; // get image identifier
// Compute image
extern struct g_pixel imFrame;
extern struct galaxie source[NFMAX];
const extern struct g_observ O;
double dx = imFrame.xmax - imFrame.xmin;
double dy = imFrame.ymax - imFrame.ymin;
extern double *sflux;
int verbose_save = M.verbose;
M.verbose = 0;
//2Eric: You assume that imFrame.pixelx == imFrame.pixely ?
if (fabs(( imFrame.pixelx - imFrame.pixely)/imFrame.pixelx) > 1e-4)
{
fprintf(stderr, "Error | imFrame.pixelx(%f) != imFrame.pixely(%f) but we assume that they are equal\n", imFrame.pixelx, imFrame.pixely);
exit(EXIT_FAILURE);
}
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);
// if (O.setseeing)
// d_seeing_omp(ero, imFrame.nx, imFrame.ny, imFrame.pixelx);
//if (O.bruit)
// d_bruiter_omp(ero, imFrame.nx, imFrame.ny);
M.verbose = verbose_save;
if (wo != NULL)
{
for (i = 0; i < imFrame.ny; i++ )
for (j = 0; j < imFrame.nx; j++ )
{
ero[i][j] -= imo[i][j];
chi_im += ero[i][j] * ero[i][j] * wo[i][j];
}
lhood0 = lhood_wo; // read from lhood_wo variable pre-computed in o_global()
}
else
{
for (i = 0; i < imFrame.ny; i++ )
for (j = 0; j < imFrame.nx; j++ )
{
ero[i][j] -= imo[i][j];
chi_im += ero[i][j] * ero[i][j] / fabs(imo[i][j] + 1.);
lhood0 += log( 2.*M_PI*fabs(imo[i][j] + 1.) );
}
}
}
else
/* chi_im is the error associated to the transformation image -> source*/
......@@ -865,9 +933,9 @@ static double chi2SglImage( struct galaxie *pima, struct point *ps )
if ( nimages > 1 )
{
if ( fabs(I.forme) == 10 )
if ( abs(I.forme) == 10 )
I2x = I2y = pima->E.a * pima->E.a;
else if ( fabs(I.forme) == 11 )
else if ( abs(I.forme) == 11 )
{
I2x = pima->E.a * cos(pima->E.theta);
I2y = pima->E.b * cos(pima->E.theta);
......
......@@ -32,10 +32,11 @@ double chi_im, chip, chix, chiy, chia, chil, chis, chi_vel,chi_mass; //I adde
//double amplifi[NFMAX][NIMAX];
double **imo, **wo, **soo, **ero, lhood_wo;
double ***cubeo,***ercubeo;
double **fluxmap;
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 nshmap;
int optim_z;
int block_s[NLMAX][NPAMAX];
......@@ -264,6 +265,7 @@ void o_global_free()
extern struct g_cube cubeFrame;
extern double **map_axx, **map_ayy;
int i;
// free the square maps
if (M.iclean == 2 )
......@@ -275,15 +277,20 @@ void o_global_free()
free_cubic_double(cubeo, cubeFrame.ny, cubeFrame.nx);
free_cubic_double(ercubeo, cubeFrame.ny, cubeFrame.nx);
}
else if (M.iclean == 4 )
{
free_square_double(ero, imFrame.ny);
free_square_double(fluxmap, imFrame.ny*imFrame.nx);
}
else
{
free_square_double(soo, ps.ny);
free_square_int(imuo, ps.ny);
}
free_square_double(imo, imFrame.ny);
free_square_double(wo, wFrame.ny);
free_square_double(fluxmap, imFrame.ny*imFrame.nx);
int pot_nopt = 0;
extern struct g_pot P[NPOTFILE];
......@@ -339,8 +346,6 @@ void o_global_free()
}
// Reset source counter to zero
// Comment out this line because of a bug with cleanlens 2 (http://projets.lam.fr/issues/2783)
// This reset is not justified
//S.ns = 0;
}
}
......@@ -358,7 +363,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 galaxie arclet[];
extern long int narclet;
......@@ -374,11 +379,8 @@ void readConstraints()
extern struct g_pixel wFrame;
extern struct g_cube cubeFrame;
extern struct galaxie source[NFMAX];
const extern struct g_pixel ps;
extern double **map_axx, **map_ayy;
//Read imframe
//Read imframe
if ( imFrame.format !=0 )
{ NPRINTF(stderr, "READ: imframe\n");
imo = (double **) readimage(&imFrame);
......@@ -406,8 +408,17 @@ void readConstraints()
// fprintf(stderr, "ERROR: O value found in %s. Cannot compute chi2.\n", imFrame.pixfile);
// exit(-1);
// }
//check for iclean = 4
if (M.iclean==4)
{//imtosou(M.zclean, ps.pixfile);
precompsource();
fprintf(stderr,"okprecomp\n");
}
if(wFrame.format != 0)
{
wo = (double **) readimage(&wFrame);
......@@ -428,6 +439,7 @@ void readConstraints()
fprintf(stderr, "ERROR: dimensions mismatch in %s and %s", wFrame.pixfile, imFrame.pixfile);
exit(-1);
}
// Pre-compute the loglikelihood normalization factor
lhood_wo = 0.;
for( i = 0; i < wFrame.nx; i++ )
......@@ -436,10 +448,11 @@ void readConstraints()
lhood_wo -= log( 2.*M_PI*wo[i][j] );
}
if ( M.iclean == 2 || M.iclean ==3 )
if ( M.iclean == 2 || M.iclean ==3 || M.iclean ==4)
{
if(M.iclean==2)
if(M.iclean==2|| M.iclean ==4)
{
// allocate temporary array to compute chi2 in o_chi()
ero = (double **) alloc_square_double(imFrame.ny, imFrame.nx);
......@@ -744,6 +757,9 @@ int getNConstraints()
if ( I.forme == 3 ) opi = 5; // +flux + shape
constraints = 0;
// Multiple images constraints
for ( i = 0; i < I.n_mult; i ++)
if ( I.mult[i] > 1 )
......@@ -782,6 +798,13 @@ int getNConstraints()
constraints += imFrame.nx * imFrame.ny;
}
//CHECK IF M.iclean = 4
if ( M.iclean == 4 )
{
extern struct g_pixel imFrame;
constraints += imFrame.nx * imFrame.ny;
}
return(constraints);
}
......@@ -793,14 +816,33 @@ int getNParameters()
extern struct g_pot P[NPOTFILE];
extern struct g_image I;
extern struct g_source S;
extern struct g_mode M;
extern int block[][NPAMAX];
extern int cblock[], sblock[NFMAX][NPAMAX];
extern struct sigposStr sigposAs;
const extern struct g_pixel ps;
long int parameters, i, j;
int ipx;
parameters = 0;
//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.;
}
}
for ( i = 0; i < G.no_lens; i ++ )
for ( ipx = CX; ipx <= PMASS; ipx ++ )
if ( block[i][ipx] != 0 )
......@@ -852,10 +894,10 @@ int getNParameters()
parameters++;
}
if ( I.dsigell != -1. ) parameters++;
return(parameters);
}
/* Initialise the np_grad and np_grad2 global variables
* used after in e_grad() and e_grad2().
*
......
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
#include<lt.h>
/****************************************************************/
/* nom: o_pixsource */
/* auteur: Johan Richard */
......@@ -14,79 +14,131 @@
/****************************************************************/
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 ;
int ni=imFrame.nx*imFrame.ny;
int ns=ps.nx*ps.ny;
double stepx,stepxs,stepy,stepys;
double ra, dec, width, height;
fluxmap=(double **)alloc_square_double(ni,ns);
double dlsds = dratio(lens[0].z, M.zclean);
double sizeratio=imFrame.pixelx/ps.pixelx; // pixel size ratio between image and source plane
double xi,yi,xs,ys,xs2,ys2;
double **tmpimage=(double **)alloc_square_double(imFrame.nx,imFrame.ny);
double xs,ys;
wcsfull(ps.wcsinfo, &ra, &dec, &width, &height);
ps.xmin -= 3600.*cos(M.ref_dec * DTR) * (ra - M.ref_ra);
ps.xmax -= 3600.*cos(M.ref_dec * DTR) * (ra - M.ref_ra);
ps.ymin -= -3600.*(dec - M.ref_dec);
ps.ymax -= -3600.*(dec - M.ref_dec);
/* fprintf(stderr,"ps.nx=%d\t ps.ny=%d\n",ps.nx,ps.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,"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, "sizeratio= %lf\n", sizeratio);
tmpimage=(double **)alloc_square_double(imFrame.ny,imFrame.nx);
int is,js,i,j;
for(is=0;is<ps.nx;is++)
{
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;
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);
// xs = is * (ps.xmax - ps.xmin) / (ps.nx - 1) + ps.xmin;
ys=ps.ymin;
for(js=0;js<ps.ny;js++)
{
xs = is * (ps.xmax - ps.xmin) / (ps.nx - 1) + ps.xmin;
ys = js * (ps.ymax - ps.ymin) / (ps.ny - 1) + ps.ymin;
{//fprintf(stderr, "js=%d\t ys=%lf\n",js,ys);
// ys = js * (ps.ymax - ps.ymin) / (ps.ny - 1) + ps.ymin;
Pi.x= imFrame.xmin;
for(i=0;i<imFrame.nx;i++)
{
Pi.y=imFrame.ymin;
// Pi.x = i * (imFrame.xmax - imFrame.xmin) / (imFrame.nx - 1) + imFrame.xmin;
for(j=0;j<imFrame.ny;j++)
{
xi = i * (imFrame.xmax - imFrame.xmin) / (imFrame.nx - 1) + imFrame.xmin;
yi = j * (imFrame.ymax - imFrame.ymin) / (imFrame.ny - 1) + imFrame.ymin;
//e_dpl(xi,yi, dlsds, xs2,ys2);
if((xs2-xs)*(xs2-xs)+(ys2-ys)*(ys2-ys)<0.25) tmpimage[i][j]=sizeratio;
else tmpimage[i][j]=0.0;
// fprintf(fp, "is=%d\t xs=%lf\t js=%d\t ys=%lf\t i=%d\t Pi.x=%lf\t j=%d\t Pi.y=%lf\n",is,xs,js,ys,i,Pi.x,j,Pi.y);
// Pi.y = j * (imFrame.ymax - imFrame.ymin) / (imFrame.ny - 1) + imFrame.ymin;
e_dpl(&Pi, dlsds,S.zs, &Ps);
z= (Ps.x-xs)*(Ps.x-xs)/stepxs/stepxs+(Ps.y-ys)*(Ps.y-ys)/stepys/stepys;
// fprintf(fp, "Pi.x=%lf\t Pi.y=%lf\t Ps.x=%lf \t Ps.y=%lf \t xs=%lf \t ys=%lf\n",Pi.x,Pi.y,Ps.x,Ps.y,xs,ys);
if(z<0.25)
{
tmpimage[j][i]=sizeratio;
// fprintf(stderr, "sizeratio\n");
}
else
tmpimage[j][i]=0.0;
Pi.y +=stepy;
}
Pi.x +=stepx;
}
//convolve tmpimage with PSF
/* convolve tmpimage with PSF */
if (O.setseeing)
d_seeing_omp(tmpimage, imFrame.nx, imFrame.ny, imFrame.pixelx);
//store tmpimage in flux matrix
/* 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[i][j];
}
fluxmap[i*imFrame.ny+j][is*ps.ny+js]=tmpimage[j][i];
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);
free_square_double(tmpimage,imFrame.ny);
}
void testpixel(double *sflux)
void pixelsource()
{
const extern struct g_observ O;
extern struct g_mode M;
extern struct g_pixel ps,imFrame;
extern double **fluxmap;
const extern struct g_frame F;
extern struct g_mode M; //M.pixfile N.npixel
double xmin, ymin, xmax, ymax, dx, dy;
double f, scale;
int nx, ny, ny_sav;
double **zz;
int i,j,k;
/*
xmin = imFrame.xmin;
xmax = imFrame.xmax;
ymin = imFrame.ymin;
ymax = imFrame.ymax;
xmin = F.xmin;
xmax = F.xmax;
ymin = F.ymin;
ymax = F.ymax;
dx = xmax - xmin;
dy = ymax - ymin;
nx = imFrame.nx; // default: assign np to nx and adjust ny
nx = M.npixel; // default: assign np to nx and adjust ny
scale = dx / (nx-1);
ny = dy / scale + 1; // warning: trunc(ny) behavior
......@@ -100,6 +152,8 @@ void testpixel(double *sflux)
zz = (double **) alloc_square_double(ny, nx);
ps=
int ns=ps.nx*ps.ny;
for(i=0;i<nx;i++)
{
......@@ -116,11 +170,81 @@ void testpixel(double *sflux)
if (O.setseeing)
d_seeing(zz, nx, ny, scale);
ny_sav = ny;
if (O.setbin)
d_binning(zz, &nx, &ny, O.bin);
if (O.bruit)
d_bruiter(zz, nx, ny);
NPRINTF(stderr, "COMP: image of arcs and arclets --> %s\n", iname);
if ( M.iref > 0 )
wrf_fits_abs("testpixel.fits", zz, nx, ny, xmin, xmax, ymin, ymax, M.ref_ra, M.ref_dec);
wrf_fits_abs(iname, zz, nx, ny, xmin, xmax, ymin, ymax, M.ref_ra, M.ref_dec);
else
wrf_fits("testpixel.fits", zz, nx, ny, xmin, xmax, ymin, ymax);
wrf_fits(iname, zz, nx, ny, xmin, xmax, ymin, ymax);
free_square_double(zz, ny_sav);
*/
}
void testpixel(double **zz,double *sflux)
{
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 f, scale;
int nx, ny;
FILE *fp;
int i,j,k;
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 = dy / scale + 1; // warning: trunc(ny) behavior
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);
zz = (double **) alloc_square_double(ny, nx);
fp = fopen("final_image1.dat","w");
for(i=0;i<imFrame.nx;i++)
{
for(j=0;j<imFrame.ny;j++)