Commit 505dc0a3 authored by Johan Richard's avatar Johan Richard
Browse files

New Newton implementation for multiple images search

parent d45d2319
......@@ -443,6 +443,7 @@ void e_transform(struct triplet *I, double dlsds, double zs, struct triplet *
void interieur(const struct triplet *E, struct triplet *I);
int unlens_bc(const struct point *Ps, struct point Bs, struct galaxie *multi, struct point *multib, int n, int n_famille, int *det_stop);
int unlens_bc_single(const struct point Ps, struct point Bs, struct galaxie *multi, struct point *vmultib, int i);
int unlens_bc_single_newton(const struct point Ps, struct point Bs, struct galaxie *multi, struct point *vmultib, int i); // newton mul search (bclement)
void updatecut(int i);
void update_emass(int i);
void update_epot(int i, double *epot);
......
......@@ -377,6 +377,7 @@ struct g_image
double Aymax;
double Amean;
double Adisp;
int newton; // newton mul search (bclement)
};
struct g_source
......
......@@ -17,14 +17,14 @@
// BCLEMENT : chi2 proximity bug fix
/* qsort struct comparision function (x field) */
static int struct_cmp_by_x(const void *a, const void *b)
{
static int struct_cmp_by_x(const void *a, const void *b)
{
struct point *da = (struct point *)a;
struct point *db = (struct point *)b;
if (da->x < db->x) return -1;
if (da->x > db->x) return 1;
return 0;
}
}
static double chi2SglImage( struct galaxie *pima, struct point *ps );
static int chi2_img( double *chi2, double *lh0 );
......@@ -1571,8 +1571,17 @@ static int chi2_img( double *chi2, double *lh0 )
else //MULTIPLE IMAGES system
{
/*return in vmultib[task_idx] an image for current arclet */
int ok = unlens_bc_single(vPs[task_idx], vBs[i], &multi[i][j],
&vmultib[task_idx], i);
int ok;
if(I.newton == 1) // newton mul search (bclement)
{
ok = unlens_bc_single_newton(vPs[task_idx], vBs[i], &multi[i][j],
&vmultib[task_idx], i);
}
else
{
ok = unlens_bc_single(vPs[task_idx], vBs[i], &multi[i][j],
&vmultib[task_idx], i);
}
//if !ok then we cannot find this image
#ifdef CHIRES
......
......@@ -17,15 +17,14 @@
// BCLEMENT : chi2 proximity bug fix
/* qsort struct comparision function (x field) */
static int struct_cmp_by_x(const void *a, const void *b)
{
static int struct_cmp_by_x(const void *a, const void *b)
{
struct point *da = (struct point *)a;
struct point *db = (struct point *)b;
if (da->x < db->x) return -1;
if (da->x > db->x) return 1;
return 0;
}
}
static double chi2SglImage( struct galaxie *pima, struct point *ps );
static int chi2_img( double *chi2, double *lh0 );
......@@ -1572,8 +1571,17 @@ static int chi2_img( double *chi2, double *lh0 )
else //MULTIPLE IMAGES system
{
/*return in vmultib[task_idx] an image for current arclet */
int ok = unlens_bc_single(vPs[task_idx], vBs[i], &multi[i][j],
&vmultib[task_idx], i);
int ok;
if(I.newton == 1) // newton mul search (bclement)
{
ok = unlens_bc_single_newton(vPs[task_idx], vBs[i], &multi[i][j],
&vmultib[task_idx], i);
}
else
{
ok = unlens_bc_single(vPs[task_idx], vBs[i], &multi[i][j],
&vmultib[task_idx], i);
}
//if !ok then we cannot find this image
#ifdef CHIRES
......
......@@ -216,6 +216,12 @@ void r_image(FILE *IN, FILE *OUT)
sscanf(third, "%lf", &I.Aseeing);
fprintf(OUT, "\t%s\t\t%lf\n", second, I.Aseeing);
}
else if (!strcmp(second, "newton")) // newton mul search (bclement)
{
sscanf(third, "%d", &I.newton);
fprintf(OUT, "\t%s\t\t%d\n", second, I.newton);
}
// Read the next line
fmot(IN, second);
......
......@@ -66,6 +66,7 @@ void set_default()
I.zarclet = -1.;
I.n_mult = 0;
I.multmode = 0;
I.newton = 0; // newton mul search (bclement)
// z_a_limit
zalim.bk = 0;
......
......@@ -51,6 +51,8 @@ int unlens_bc_single(struct point Psi,
struct point *multibi,
int n_famille);
static int newton_method(struct point im, struct point vBs, double dl0s, double dlsds, double conv, double z, double distmin, struct point *imfound); // newton mul search (bclement)
//det_stop: if det_stop than we should stop imediately
// It can happen in openmp mode then chi2_img
// have already decided return -1.
......@@ -214,3 +216,113 @@ int unlens_bc_single(struct point Psi,
return 0; //do not count the current image
}
}
//function for single arclet using Newton's Method (bclement)
int unlens_bc_single_newton(struct point Psi,
struct point Bs,
struct galaxie *multii,
struct point *multibi,
int n_famille)
{
const extern double distmin[NFMAX];
int error;
int ok;
double dl0s,dos,dlsds,conv;
dl0s = multii->dl0s;
dos = multii->dos;
dlsds = multii->dr;
conv = 1./ dos;
struct point imfound;
error = newton_method(multii->C, Bs, dl0s, dlsds, conv, multii->z, distmin[n_famille], &imfound);
if ( error == 0 )
{
*multibi = imfound;
ok = 1;
}
else // Newton method didn't converge
{
*multibi = multii->C;
//NPRINTF(stderr,"WARNING: could not find the searched image\n");
// NPRINTF(stderr,"Ai=%.3lf Aig=%.3lf D:%lf i:%d j:%d it:%d\n",1./e_amp(multii->C,dlsds),1./e_amp(CG,dlsds),D,n_famille,i,it);
ok = 0;
}
return ok;
}
static int newton_method(struct point im, struct point vBs, double dl0s, double dlsds, double conv, double z, double distmin, struct point *imfound)
{
struct point dpl;
struct point p;
struct point p0;
struct point d;
struct matrix J;
struct matrix Jinv;
struct matrix grad2;
int it = 0;
int maxit = 200;
dpl = e_grad(&im, dlsds, z); // already computed this when computing the barycentre
dpl.x *= dlsds;
dpl.y *= dlsds;
p.x = im.x - vBs.x - dpl.x; //f(x,y)
p.y = im.y - vBs.y - dpl.y; //g(x,y)
p0.x = p.x;
p0.y = p.y;
while(sqrt(p.x*p.x + p.y*p.y) > 1E-8)
{
it++;
if(it == maxit)
{
return 1;
}
// we want to find p.x = 0 && p.y = 0 using Newton's Method:
// dx = -Jf^-1(x)f(x)
// Jf is the Jacobian of p(x,y)
// Jf = | df(x,y)/dx df(x,y)/dy) | = | (1 - d2phixx) (-d2phiyx) |
// | dg(x,y)/dx dg(x,y)/dy) | | (-d2phiyx) (1 - d2phiyy) |
//
// Jf = | (1 - g2.a) (-g2.b) |
// | (-g2.b) (1 - g2.c) |
grad2 = e_grad2(&im, dl0s, z);
J.a = 1 - conv * grad2.a;
J.b = - conv * grad2.b;
J.c = 1 - conv * grad2.c;
J.d = - conv * grad2.b;
// inverse J and multiply by -1
double det = 1/(J.a*J.c - J.b*J.d);
Jinv.a = -J.c*det;
Jinv.b = J.b*det;
Jinv.c = -J.a*det;
Jinv.d = J.d*det;
// multiply by p
d.x = Jinv.a*p.x + Jinv.b*p.y;
d.y = Jinv.d*p.x + Jinv.c*p.y;
// adaptive first jump
// if jump too large, make it smaller according to parameter distmin
double f = 1;
if(it == 1 && sqrt(d.x*d.x+d.y*d.y) >= distmin/2. )
f = sqrt(d.x*d.x+d.y*d.y)/(distmin/2.);
im.x = im.x + d.x/f;
im.y = im.y + d.y/f;
dpl = e_grad(&im, dlsds, z);
dpl.x *= dlsds;
dpl.y *= dlsds;
p.x = im.x - vBs.x - dpl.x; //f(x,y)
p.y = im.y - vBs.y - dpl.y; //g(x,y)
}
*imfound = im;
return 0;
}
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