Skip to content
GitLab
Commit 4cc0cc1d authored by Johan Richard's avatar Johan Richard
Browse files

Add optimisation within contours for shapemodels

parent b5cf3051
Hide whitespace changes
Inline Side-by-side
include/fonction.h
View file @ 4cc0cc1d
......@@ -221,6 +221,7 @@ int inverse(const struct point gsource[][NGGMAX], struct point *P, struct bi
void isoima(struct ellipse *es, struct ellipse *ampli, struct ellipse *ei);
double iter(double phi0, double coeur, double ct);
void keep_cl(double **image, struct pixlist *pl, int *npl);
void keep_cl_extended(struct pixlist *pl, int npl, struct pixlist *pl_ext, int *npl_ext); // mod bclement shapemodel_contour_ext
int lire(struct galaxie *liste, char *name);
complex lncpx(complex c);
void local(int type, char localfile[]);
......@@ -284,6 +285,7 @@ void o_mag_m(int n, struct galaxie *arclet);
double o_min_loc(double y0);
double o_min_slope(double y0);
void o_pixel(double **zz, int nx, int ny, double scale, double xmin, double xmax, double ymin, double ymax, struct galaxie *source, double **dpl_x, double **dpl_y);
void o_pixel_contour(double **zz, int nx, int ny, double scale, double xmin, double xmax, double ymin, double ymax, struct galaxie *source, struct pixlist *pl, struct pixlist *pl_dpl, int npl); // mod bclement shapemodel_contour
double o_prep();
void o_prep_mult(int ntmult, struct galaxie *mult);
void o_print(FILE *OUT, double chi0);
......
include/structure.h
View file @ 4cc0cc1d
......@@ -817,6 +817,8 @@ struct pixlist
{
double i;
double j;
int ii; // mod bclement shapemodel_contour
int jj; // mod bclement shapemodel_contour
double flux;
};
......
src/checkpar.c
View file @ 4cc0cc1d
......@@ -8,20 +8,14 @@
/****************************************************************/
/* nom: checkpar */
/* auteur: Johan Richard */
/* date: 12/12/17 */
/* date: 25/06/19 */
/* place: Lyon */
/* */
/*
* Check the consistency of parameters and constraints in the .par
*
* List of points to check (non-exhaustive):
* - number of optimised z_mult vs number of multiple images without redshifts
* (already implemented but message should be more explicit on the non-matched systems).
* - check if two images in one or several systems have exactly the same location (should not
* normally happen, but maybe could be override by the user?)
* - redshift differences in the same system ID.
* number of lenses (nlens) vs number of potentials (Warning, re-adjustment, confirmation by the user...)
* check if some wcs coordinates are completely off the rest (fixed limit ~ 5 arcmin? or field limits?, give a warning)
* Sources with sigx + eps => sigy adjusted
*
****************************************************************
*/
......@@ -30,14 +24,31 @@ void checkpar()
const extern struct pot lens[];
const extern struct g_grille G;
const extern struct g_mode M;
extern struct g_source S;
extern struct galaxie source[NFMAX];
FILE *OUT;
int i;
double aa, bb, Cx, Cy;
//Double check galaxy sources
for (i=0; i<S.ns;i++)
{
if(source[i].eps>0)
{
if(source[i].E.a>0)
{
if(source[i].E.b>0)
fprintf(stderr, "WARNING: source %d sig_y parameter readjusted according to eps\n", i);
source[i].E.b = source[i].E.a * (1. - source[i].eps) / (1. + source[i].eps);
}
else
{
fprintf(stderr, "ERROR: in source %d eps parameter is set without sigx\n", i);
exit(-1);
}
}
}
return;
......
src/d_seeing_omp.c
View file @ 4cc0cc1d
......@@ -471,7 +471,7 @@ void d_seeing_omp(double **im, int ny, int nx, double scale)
if (nx != nx_init || ny != ny_init || scale != scale_init)
{
fprintf(stderr, "Error | d_seeing_omp.c/d_seeing_omp | nx, ny or scale have been changed from first run");
exit(EXIT_FAILURE);
//exit(EXIT_FAILURE);
}
if (O.setseeing != O_psf_f.setseeing ||
O.seeing != O_psf_f.seeing ||
......@@ -480,7 +480,7 @@ void d_seeing_omp(double **im, int ny, int nx, double scale)
O.seeing_angle != O_psf_f.seeing_angle )
{
fprintf(stderr, "Error | d_seeing_omp.c/d_seeing_omp | O structure have been changed from first run");
exit(EXIT_FAILURE);
//exit(EXIT_FAILURE);
}
}
......
src/e_pixel.c
View file @ 4cc0cc1d
......@@ -71,9 +71,8 @@ void e_pixel(int np, char *iname, char *sname, struct galaxie *source)
o_pixel(zz, nx, ny, scale, xmin, xmax, ymin, ymax, source, NULL, NULL);
if (O.setseeing)
d_seeing_omp(zz, nx, ny, scale);
d_seeing_omp(zz, nx, ny, scale); // mod bclement shapemodel_opt
ny_sav = ny;
if (O.setbin)
......@@ -126,6 +125,8 @@ void e_pixel(int np, char *iname, char *sname, struct galaxie *source)
zz[j][k] = 0.;
for( l = 0; l < S.ns; l++ )
zz[j][k] += d_profil(A.x, A.y, &source[l]);
}
}
......@@ -151,29 +152,54 @@ void o_pixel(double **zz, int nx, int ny, double scale, double xmin,
const extern struct g_mode M;
const extern struct g_source S;
extern struct pot lens[];
int j, k, l;
struct point pi, ps;
for (j = 0; j < ny; j++)
{
for (k = 0; k < nx; k++)
{
zz[j][k] = 0.0;
zz[j][k] = 0.0;
}
}
//#pragma omp parallel for private(ps,pi,j,k,l)
// mod bclement shapemodel_opt
// dpl_x and dpl_y are only assigned if lens model is fixed and pre-computed for a zs_ref=100 in o_global()
// TODO : save computation time of this integral by storing the value dlsds_ref somewhere
double dlsds_ref = 1;
if(dpl_x)
{
double zs_ref = 100.0;
dlsds_ref = dratio(lens[0].z, zs_ref);
}
//#pragma omp parallel for private(ps,pi,j,k,l)
for (l = 0 ; l < S.ns ; l++ )
{
// TODO : can use collapse(2) or collapse(3) and put pragma before the source loop
#pragma omp parallel for schedule(static,1) private(ps,pi,j,k)
for (j = 0; j < ny; j++)
{
pi.y = ymin + j * scale;
pi.y = ymin + j * scale;
for (k = 0; k < nx; k++)
{
pi.x = xmin + k * scale;
e_dpl(&pi,source[l].dr,source[l].z,&ps);
// mod bclement shapemodel_opt
if(dpl_x)
{
double fdr;
fdr = source[l].dr / dlsds_ref;
ps.x = pi.x - dpl_x[j][k]*fdr;
ps.y = pi.y - dpl_y[j][k]*fdr;
}
else
{
e_dpl(&pi,source[l].dr,source[l].z,&ps);
}
double f = d_profil(ps.x, ps.y, &source[l]);
/* if ((f>sqrt(O.SKY/O.gain)) || (f>source[l].I0/100.))
zz[j][k]+=d_integrer(source[l],pi.x,pi.y,scale/2.,f,0);
......@@ -186,6 +212,66 @@ void o_pixel(double **zz, int nx, int ny, double scale, double xmin,
}
// mod bclement shapemodel_contour
/****************************************************************/
/* nom: o_pixel_contour */
/* auteur: Benjamin Clement */
/* date: 13/02/19 */
/* place: Lyon/Geneva */
/****************************************************************
* Simulate an image from a source inside a user-defined list of contours
* Called by o_chi() and o_chires()
*/
void o_pixel_contour(double **zz, int nx, int ny, double scale, double xmin,
double xmax, double ymin, double ymax, struct galaxie *source, struct pixlist *pl, struct pixlist *pl_dpl, int npl)
{
const extern struct g_mode M;
const extern struct g_source S;
const extern struct g_observ O;
extern struct pot lens[];
int j, k, l, i, n;
struct point pi, ps;
// TODO : only zero the array inside and around contours
for (j = 0; j < ny; j++)
{
for (k = 0; k < nx; k++)
{
zz[j][k] = 0.0;
}
}
for (l = 0 ; l < S.ns ; l++ )
{
// TODO : can use collapse(2) or collapse(3) and put pragma before the source loop ?
#pragma omp parallel for schedule(static) private(ps,pi,j,k)
for(n = 0; n < npl; n++)
{
pi.y = pl[n].i;
pi.x = pl[n].j;
j = pl[n].ii;
k = pl[n].jj;
// compute the source position from pre-computed dpl in o_global
if(pl_dpl)
{
ps.x = pi.x - pl_dpl[n].i*source[l].dr;
ps.y = pi.y - pl_dpl[n].j*source[l].dr;
}
else // recompute dpl
{
e_dpl(&pi,source[l].dr,source[l].z,&ps);
}
double f = d_profil(ps.x, ps.y, &source[l]);
zz[j][k] += f;
}
}
}
/****************************************************************/
/* nom: pixelsource */
/* auteur: Johan Richard */
......
src/keep_cl.c
View file @ 4cc0cc1d
......@@ -121,6 +121,8 @@ void keep_cl(double **image, struct pixlist *pl, int *npl)
{
pl[kk].i = P.y;
pl[kk].j = P.x;
pl[kk].ii = ii; // mod bclement shapemodel_contour
pl[kk].jj = jj; // mod bclement shapemodel_contour
pl[kk++].flux = image[ii][jj];
clean[ii][jj] = image[ii][jj];
}
......
src/o_chi.c
View file @ 4cc0cc1d
......@@ -91,6 +91,12 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
// extern struct pixlist plo[];
extern double **imo, **wo, **ero, **soo, lhood_wo;
// mod bclement shapemodel_contour
extern struct pixlist *imo_pl;
extern struct pixlist *imo_pl_dpl;
extern int imo_npl;
int n;
/* local variables */
struct point A, B, D;
struct ellipse ampli;
......@@ -183,10 +189,10 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
fprintf(stderr, "Error | imFrame.pixelx(%f) != imFrame.pixely(%f) but we assume that they are equal\n", imFrame.pixelx, imFrame.pixely);
exit(EXIT_FAILURE);
}
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);
// mod bclement shapemodel_contour
//o_pixel(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, map_axx, map_ayy);
o_pixel_contour(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, imo_pl, imo_pl_dpl, imo_npl); // mod bclement shapemodel_contour
if (O.setseeing)
d_seeing_omp(ero, imFrame.nx, imFrame.ny, imFrame.pixelx);
......@@ -195,27 +201,38 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
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];
}
// mod bclement shapemodel_contour
// compute the chi2 only inside the contours
//for (i = 0; i < imFrame.ny; i++ )
//for (j = 0; j < imFrame.nx; j++ )
for(n = 0; n < imo_npl; n++)
{
i = imo_pl[n].ii;
j = imo_pl[n].jj;
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.) );
}
}
// mod bclement shapemodel_contour
// compute the chi2 only inside the contours
//for (i = 0; i < imFrame.ny; i++ )
// for (j = 0; j < imFrame.nx; j++ )
for(n = 0; n < imo_npl; n++)
{
i = imo_pl[n].ii;
j = imo_pl[n].jj;
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.) );
}
}
// restore relative source positions
if( I.n_mult > 0 )
......
src/o_chires.c
View file @ 4cc0cc1d
......@@ -91,6 +91,12 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
// extern struct pixlist plo[];
extern double **imo, **wo, **ero, **soo, lhood_wo;
// mod bclement shapemodel_contour
extern struct pixlist *imo_pl;
extern struct pixlist *imo_pl_dpl;
extern int imo_npl;
int n;
/* local variables */
struct point A, B, D;
struct ellipse ampli;
......@@ -183,10 +189,10 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
fprintf(stderr, "Error | imFrame.pixelx(%f) != imFrame.pixely(%f) but we assume that they are equal\n", imFrame.pixelx, imFrame.pixely);
exit(EXIT_FAILURE);
}
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);
// mod bclement shapemodel_contour
//o_pixel(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, map_axx, map_ayy);
o_pixel_contour(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, imo_pl, imo_pl_dpl, imo_npl); // mod bclement shapemodel_contour
if (O.setseeing)
d_seeing_omp(ero, imFrame.nx, imFrame.ny, imFrame.pixelx);
......@@ -195,27 +201,38 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
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];
}
// mod bclement shapemodel_contour
// compute the chi2 only inside the contours
//for (i = 0; i < imFrame.ny; i++ )
//for (j = 0; j < imFrame.nx; j++ )
for(n = 0; n < imo_npl; n++)
{
i = imo_pl[n].ii;
j = imo_pl[n].jj;
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.) );
}
}
// mod bclement shapemodel_contour
// compute the chi2 only inside the contours
//for (i = 0; i < imFrame.ny; i++ )
// for (j = 0; j < imFrame.nx; j++ )
for(n = 0; n < imo_npl; n++)
{
i = imo_pl[n].ii;
j = imo_pl[n].jj;
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.) );
}
}
// restore relative source positions
if( I.n_mult > 0 )
......
src/o_global.c
View file @ 4cc0cc1d
......@@ -52,6 +52,11 @@ int nwarn; // warning emitted during image plane chi2
double *np_b0; // non parametric accelerator
//int nplo; // number of points for criticinv.c DEPRECATED
// mod bclement shapemodel_contour
struct pixlist *imo_pl; // list of points inside the polygons for shapemodel
struct pixlist *imo_pl_dpl; // list of dpl inside the polygons for shapemodel
int imo_npl;
/****************************************************************/
/* nom: global */
/* auteur: Jean-Paul Kneib */
......@@ -272,6 +277,8 @@ void o_global_free()
if (M.iclean == 2 )
{
free_square_double(ero, imFrame.ny);
free(imo_pl); // mod bclement shapemodel_contour
free(imo_pl_dpl); // mod bclement shapemodel_contour
}
else if (M.iclean == 3 )
{
......@@ -460,6 +467,47 @@ void readConstraints()
// allocate temporary array to compute chi2 in o_chi()
ero = (double **) alloc_square_double(imFrame.ny, imFrame.nx);
}
// mod bclement shapemodel_contour
// Initialize and build the pixlist defining the contours
if(M.iclean == 2)
{
/* Allocate the maximum number of points to send to source plane */
imo_pl = (struct pixlist*)malloc((unsigned int)
imFrame.nx * imFrame.ny * sizeof(struct pixlist));
/* keep only the relevant points */
if (imFrame.ncont == 0)
{
NPRINTF(stderr, "ERROR: the number of contours is not defined. Set the keyword ncont in the parameter file.\n",imFrame.ncont);
exit(-1);
}
else
{
keep_cl(imo, imo_pl, &imo_npl);
}
// TODO : remove duplicate pixels in list when overlapping contours occur
// This would cause a bug o_pixel_contour()
// TODO possibly reallocate imo_pl given imo_npl
imo_pl = (struct pixlist*)realloc(imo_pl, (unsigned int)
imo_npl * sizeof(struct pixlist));
// redo the pre-computation of the loglikelihood normalization factor
lhood_wo = 0.;
int n;
for( n = 0; n < imo_npl; n++ )
{
j = imo_pl[n].ii;
i = imo_pl[n].jj;
if( wo[j][i] > 0 )
lhood_wo -= log( 2.*M_PI*wo[j][i] );
}
}
if(M.iclean==3)
{
// allocate temporary array to compute chi2 in o_chi()
......@@ -482,22 +530,41 @@ void readConstraints()
if( G.no_lens == 0 && pot_nopt == 0 )
{
map_axx = (double **) alloc_square_double(imFrame.ny, imFrame.nx);
map_ayy = (double **) alloc_square_double(imFrame.ny, imFrame.nx);
//mod bclement shapemodel_contour
struct point pi, ps;
double zs_ref = 100.;
double dlsds_ref = dratio(lens[0].z, zs_ref);
for (i = 0; i < imFrame.ny; i++)
{
pi.y = imFrame.ymin + i * imFrame.pixelx;
for (j = 0; j < imFrame.nx; j++)
{
pi.x = imFrame.xmin + j * imFrame.pixelx;
e_dpl(&pi, dlsds_ref, zs_ref, &ps);
map_axx[i][j] = pi.x - ps.x;
map_ayy[i][j] = pi.y - ps.y;
}
}
if(M.iclean == 2) // compute normalized dpl only inside contour
{
imo_pl_dpl = (struct pixlist*)malloc((unsigned int)
imo_npl * sizeof(struct pixlist));
int n;
for(n = 0; n < imo_npl; n++)
{
pi.y = imo_pl[n].i;
pi.x = imo_pl[n].j;
e_dpl(&pi, dlsds_ref, zs_ref, &ps);
imo_pl_dpl[n].i = (pi.x - ps.x)/dlsds_ref;
imo_pl_dpl[n].j = (pi.y - ps.y)/dlsds_ref;
}
}
else
{
map_axx = (double **) alloc_square_double(imFrame.ny, imFrame.nx);
map_ayy = (double **) alloc_square_double(imFrame.ny, imFrame.nx);
for (i = 0; i < imFrame.ny; i++)
{
pi.y = imFrame.ymin + i * imFrame.pixelx;
for (j = 0; j < imFrame.nx; j++)
{
pi.x = imFrame.xmin + j * imFrame.pixelx;
e_dpl(&pi, dlsds_ref, zs_ref, &ps);
map_axx[i][j] = pi.x - ps.x;
map_ayy[i][j] = pi.y - ps.y;
}
}
}
}
// read and prepare sources for o_chi()
......@@ -749,6 +816,9 @@ int getNConstraints()
extern struct g_image I;
extern struct g_dyn Dy; //TV Oct 2011
extern long int narclet;
extern struct pixlist *imo_pl; // mod bclement shapemodel_contour
extern double **wo; // mod bclement shapemodel_contour
extern int imo_npl; // mod bclement shapemodel_contour
int constraints, i;
int opi; // number of observable per images
......@@ -795,8 +865,20 @@ int getNConstraints()
// FITS image reconstruction
if ( M.iclean == 2 )
{
extern struct g_pixel imFrame;
constraints += imFrame.nx * imFrame.ny;
//extern struct g_pixel imFrame;
//constraints += imFrame.nx * imFrame.ny;
// mod bclement shapemodel_contour
// Number of constraints is now the number
// of non-zero weight pixels inside the contours
int n,j;
for( n = 0; n < imo_npl; n++ )
{
j = imo_pl[n].ii;
i = imo_pl[n].jj;
if( wo[j][i] > 0 )
constraints++;
}
}
//CHECK IF M.iclean = 4
......
src/o_prep_mult.c
View file @ 4cc0cc1d
......@@ -77,8 +77,8 @@ void o_prep_mult(int ntmult, struct galaxie *mult)
{
if(j>=NIMAX-2)
{
fprintf(stderr, "ERROR: too many images (>%d) in system %s\n",NIMAX-1,multi[k][0].n);
fprintf(stderr,"Please increase the NIMAX parameter in dimension.h\n");
NPRINTF(stderr, "ERROR: too many images (>%d) in system %s\n",NIMAX-1,multi[k][0].n);
NPRINTF(stderr,"Please increase the NIMAX parameter in dimension.h\n");
exit(-1);
}
multi[k][++j] = mult[i];
......@@ -98,7 +98,7 @@ void o_prep_mult(int ntmult, struct galaxie *mult)
if ( k >= NFMAX && i < ntmult )
{
fprintf(stderr, "ERROR: too many systems in %s (maximum %d)\n",
NPRINTF(stderr, "ERROR: too many systems in %s (maximum %d)\n",
I.multfile, NFMAX);
exit(-1);
}
......@@ -120,7 +120,7 @@ void o_prep_mult(int ntmult, struct galaxie *mult)
// no single images without redshift
if ( I.mult[i] == 1 && multi[i][0].z == 0 )
{
fprintf(stderr, "ERROR: no multiple images in file %s\n",
NPRINTF(stderr, "ERROR: no multiple images in file %s\n",
I.multfile);
exit(-1);
};
......@@ -165,7 +165,7 @@ void o_prep_mult(int ntmult, struct galaxie *mult)
//checks differences in redshifts in a same system
if(multi[i][j].z!=multi[i][jcheck].z)