Home
Manual
Packages
Global Index
Keywords
Quick Reference
|
#if 0
u = double(indgen(100));
v = double(indgen(105))(-,);
data = 20*exp(-4*log(2)*(((u-21.7)/4.1)^2 + ((v-60.33)/4.1)^2)) + 10*(random(numberof(u), numberof(v)) - 0.5);
gg_fit_2d_spike(data);
#endif
require, "yeti_gist_gui.i";
func gg_draw_box (x0, y0, x1, y1, color=, width=)
{
plg, color=color, type=1, width=width, marks=0,
[y0,y1,y1,y0,y0], [x0,x0,x1,x1,x0];
}
func gg_draw_target (xc, yc, dx, dy, color=, width=)
{
x0 = (x1 = xc - 0.5*dx) - dx;
x3 = (x2 = xc + 0.5*dx) + dx;
y0 = (y1 = yc - 0.5*dy) - dy;
y3 = (y2 = yc + 0.5*dy) + dy;
//write, "target:",xc,yc," - ",plsys();
pldj, color=color, type=1, width=width,
[x0,x2,x1,x1,xc,xc,x1,x2], [yc,yc,y1,y2,y0,y2,y1,y1],
[x1,x3,x2,x2,xc,xc,x1,x2], [yc,yc,y1,y2,y1,y3,y2,y2];
}
func gg_draw_roi (xc, yc, dx, dy, size, color=, width=)
{
w = (size - 1)/2;
wx = dx*(w + 0.5);
wy = dy*(w + 0.5);
x1 = xc - wx;
x2 = xc - 0.5*dx;
x3 = xc + 0.5*dx;
x4 = xc + wx;
y1 = yc - wy;
y2 = yc - 0.5*dy;
y3 = yc + 0.5*dy;
y4 = yc + wy;
pldj, color=color, type=1, width=width,
[x1,x4,x4,x1,x2,x3,x3,x2], [y1,y1,y4,y4,y2,y2,y3,y3],
[x4,x4,x1,x1,x3,x3,x2,x2], [y1,y4,y4,y1,y2,y3,y3,y2];
}
func gg_fit_2d_spike (z, x0, y0, x1, y1, background=, size=)
{
local model, param;
dims = dimsof(z);
nx = dims(2);
ny = dims(3);
if (is_void(x0)) x0 = 1.0;
if (is_void(y0)) y0 = 1.0;
if (is_void(x1)) x1 = x0 + (nx - 1);
if (is_void(y1)) y1 = y0 + (ny - 1);
xmin = double(x0);
xmax = double(x1);
ymin = double(y0);
ymax = double(y1);
xstep = (xmax - xmin)/(nx - 1);
ystep = (ymax - ymin)/(ny - 1);
x = xmin + xstep*indgen(0:nx-1);
y = ymin + ystep*indgen(0:ny-1);
x0 -= 0.5*xstep;
x1 += 0.5*xstep;
y0 -= 0.5*ystep;
y1 += 0.5*ystep;
/* Create widgets. */
bw = 77; // button width
bh = 25; // button height
pad = 2; // padding space
bd = 2; // border width
//Button = gg_button(text=string(0), width=bw, height=bh, bd=bd);
Top = gg_toplevel(bd=bd, relief="groove", padx=pad, pady=pad);
Zoom = gg_button(text="zoom", width=bw, height=bh, bd=bd);
Unzoom = gg_button(text="unzoom", width=bw, height=bh, bd=bd);
Verbose = gg_button(text="verbose", width=bw, height=bh, bd=bd);
Redraw = gg_button(text="redraw", width=bw, height=bh, bd=bd);
Ok = gg_button(text="ok", width=bw, height=bh, bd=bd);
Cancel = gg_button(text="cancel", width=bw, height=bh, bd=bd);
Fit = gg_button(text="fit", width=bw, height=bh, bd=bd);
Message = gg_label(text="this is the help...", width=7*bw, height=5*bh,
bd=1, relief=GG_SUNKEN, bg=GG_WHITE, fg=GG_BLUE,
font="helveticaB", fontsize=18, justify="CH");
/* Widgets for "size" selection. */
size_max = min(nx, ny);
size_min = 3;
if (is_void(size)) size = 11;
if (numberof(size)!=1 || dimsof(size)(1)!=0 || structof((size += 0)) != long ||
size%2!=1) error, "bad value for keyword SIZE";
if (size < size_min) error, "value of SIZE keyword too small";
SizeLabel = gg_label(text="size:", width=bw, height=bh,
bd=bd, relief=GG_FLAT, justify="RH");
SizeValue = gg_label(text=swrite(format="%d", size), width=bw, height=bh,
bd=1, relief=GG_SUNKEN, justify="LH",
bg=GG_WHITE, fg=GG_BLACK);
SizePlus = gg_button(text="+", width=12, height=12, bd=1);
SizeMinus = gg_button(text="-", width=12, height=12, bd=1);
/* Widgets for "background" selection. */
Background = gg_label(text="background:", width=bw, height=bh,
bd=bd, relief="flat", justify="RH");
BackgroundNone = gg_button(text="none",
relief=GG_RAISED, fg=GG_DEFAULT_FG,
width=bw, height=bh, bd=bd);
BackgroundConstant = gg_button(text="constant",
relief=GG_RAISED, fg=GG_DEFAULT_FG,
width=bw, height=bh, bd=bd);
BackgroundSlope = gg_button(text="slope",
relief=GG_RAISED, fg=GG_DEFAULT_FG,
width=bw, height=bh, bd=bd);
if (is_void(background)) background = 0;
if (background == 0) selection = BackgroundNone;
else if (background == 1) selection = BackgroundConstant;
else if (background == 2) selection = BackgroundSlope;
else error, "bad value for keyword BACKGROUND";
h_set, selection, relief=GG_SUNKEN, fg=GG_RED;
/* Arrange widgets in a grid. */
gg_manage, Top, 1,1, Zoom;
gg_manage, Top, 2,1, Unzoom;
gg_manage, Top, 3,1, Verbose;
gg_manage, Top, 6,1, Redraw;
gg_manage, Top, 1,2, Background, anchor="e";
gg_manage, Top, 2,2, BackgroundNone;
gg_manage, Top, 3,2, BackgroundConstant;
gg_manage, Top, 4,2, BackgroundSlope;
gg_manage, Top, 1,3, SizeLabel, rowspan=2;
gg_manage, Top, 2,3, SizeValue, rowspan=2;
gg_manage, Top, 3,3, SizePlus, anchor="ws", padx=1, pady=1;
gg_manage, Top, 3,4, SizeMinus, anchor="wn", padx=1, pady=1;
gg_manage, Top, 6,3, Fit, rowspan=2;
gg_manage, Top, 1,5, Message, colspan=6;
gg_manage, Top, 1,6, Ok;
gg_manage, Top, 6,6, Cancel;
fma;
pli, z, x0, y0, x1, y1;
limits;
gg_paint, Top;
/* Event loop. */
npaints = 0;
do_zoom = 0;
verbose = 0;
point = xpoint = ypoint = 0;
for (;;) {
ms = mouse(-1, 0, "");
if (is_void(ms)) continue;
clicked = gg_find_clicked(Top, ms);
if (is_void(clicked)) {
if (do_zoom) {
gg_zoom, ms;
} else {
lim = limits();
if (long(ms(9))==1 &&
(ms(1) - lim(1))*(ms(1) - lim(2)) < 0.0 &&
(ms(2) - lim(3))*(ms(2) - lim(4)) < 0.0 &&
(ms(1) - x0)*(ms(1) - x1) < 0.0 &&
(ms(2) - y0)*(ms(2) - y1) < 0.0) {
if (point) {
/* Erase previous mark. */
gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
}
xpoint = long(floor((ms(1) - xmin)/xstep + 1.5));
ypoint = long(floor((ms(2) - ymin)/ystep + 1.5));
gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
point = 1;
++npaints;
}
}
} else if (clicked == BackgroundNone) {
gg_paint, h_set(BackgroundNone, relief=GG_SUNKEN, fg=GG_RED);
gg_paint, h_set(BackgroundConstant, relief=GG_RAISED, fg=GG_DEFAULT_FG);
gg_paint, h_set(BackgroundSlope, relief=GG_RAISED, fg=GG_DEFAULT_FG);
background = 0;
++npaints;
} else if (clicked == BackgroundConstant) {
gg_paint, h_set(BackgroundNone, relief=GG_RAISED, fg=GG_DEFAULT_FG);
gg_paint, h_set(BackgroundConstant, relief=GG_SUNKEN, fg=GG_RED);
gg_paint, h_set(BackgroundSlope, relief=GG_RAISED, fg=GG_DEFAULT_FG);
background = 1;
++npaints;
} else if (clicked == BackgroundSlope) {
gg_paint, h_set(BackgroundNone, relief=GG_RAISED, fg=GG_DEFAULT_FG);
gg_paint, h_set(BackgroundConstant, relief=GG_RAISED, fg=GG_DEFAULT_FG);
gg_paint, h_set(BackgroundSlope, relief=GG_SUNKEN, fg=GG_RED);
background = 2;
++npaints;
} else if (clicked == SizePlus) {
if (size < size_max) {
if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
size = min(size+2, size_max);
gg_paint, h_set(SizeValue, text=swrite(format="%d", size));
gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
++npaints;
}
} else if (clicked == SizeMinus) {
if (size > size_min) {
if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
size = max(size-2, size_min);
gg_paint, h_set(SizeValue, text=swrite(format="%d", size));
gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
++npaints;
}
} else if (clicked == Fit) {
if (! point) {
gg_paint, h_set(Message, text="You must select a position first:\nunselect \"zoom\" and\nclick mouse near spike in image.");
} else {
/* Region of interest (ROI). */
w = (size - 1)/2;
xrng = max(1, xpoint-w):min(nx, xpoint+w);
yrng = max(1, ypoint-w):min(ny, ypoint+w);
zroi = z(xrng,yrng);
xroi = x(xrng);
yroi = y(yrng)(-,);
param = [max(zroi), w, x(xpoint), y(ypoint)];
if (background > 0) {
grow, param, min(zroi);
if (background == 2) grow, param, 0.0, 0.0;
}
param = fit_2d_spike(zroi, yroi, xroi, /*model, weight=, */
shape="spline",
guess=param, bg=background,
absorption=0, verb=verbose);
xpeak = param(3);
ypeak = param(4);
if (abs(x(xpoint) - xpeak) > 0.5*xstep ||
abs(y(ypoint) - ypeak) > 0.5*ystep) {
gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
xpoint = long(floor((xpeak - xmin)/xstep + 1.5));
ypoint = long(floor((ypeak - ymin)/ystep + 1.5));
gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
msg = "*** WARNING *** box has moved: click fit again.\n";
} else {
msg = "";
}
fwhm = param(2);
msg += swrite(format="peakmax = %g fwhm = %g (%.1f×%.1f)",
param(1), fwhm, fwhm/xstep, fwhm/ystep);
msg += swrite(format="\nxpeak = %g ypeak = %g", xpeak, ypeak);
if (background > 0) {
msg += swrite(format="\nbackground = %g", param(5));
if (background == 2)
msg += swrite(format="+ %g*(x - xpeak) + %g*(y - ypeak)",
param(6), param(7));
}
gg_paint, h_set(Message, text=msg);
}
++npaints;
} else if (clicked == Zoom) {
do_zoom = ! do_zoom;
if (do_zoom) h_set, clicked, relief=GG_SUNKEN, fg=GG_RED;
else h_set, clicked, relief=GG_RAISED, fg=GG_DEFAULT_FG;
gg_paint, clicked;
++npaints;
} else if (clicked == Verbose) {
verbose = ! verbose;
if (verbose) h_set, clicked, relief=GG_SUNKEN, fg=GG_RED;
else h_set, clicked, relief=GG_RAISED, fg=GG_DEFAULT_FG;
gg_paint, clicked;
++npaints;
} else if (clicked == Unzoom || clicked == Redraw) {
fma;
pli, z, x0, y0, x1, y1;
if (clicked == Unzoom) limits;
gg_paint, Top;
if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
npaints = 0;
} else if (clicked == Ok) {
return param;
} else if (clicked == Cancel) {
return;
}
if (npaints > 30) {
fma;
pli, z, x0, y0, x1, y1;
gg_paint, Top;
if (point) gg_draw_roi, x(xpoint), y(ypoint), xstep, ystep, size,
color=GG_XOR;
npaints = 0;
}
}
}
func __fit_spike_gauss (&p, &dp, u) {
/* Gaussian shape: exp(-C1*u*u)
* To have FWHM = 1, use:
* C1 = 4*log(2) ~ 2.77258872223978123766892848583
* C2 = 2*C1 ~ 5.54517744447956247533785697167
*/
p = exp(-2.77258872223978123766892848583*u*u);
dp = -5.54517744447956247533785697167*u*p;
}
/*
* CUBIC B-SPLINE BASIS FUNCTION
* S(x) = (|x|/2 - 1)*x^2 + 2/3 if |x|<=1
* = 1/6*(2 - |x|)^3 if 1<=|x|<=2
* = 0 otherwise
* S(x) is normalised (integral S(x) dx = 1) but we want a function
* F(x) = alpha*S(beta*x) so that:
* f(0) = 1
* f(0.5) = 0.5
* f(-u) = f(u) (i.e. f is symmetric)
* f(u) >= 0 everywhere
* f'(u)>0 for u<0, f'(0)=0, f'(u)<0 for u>0 (i.e. f is "bell-shaped")
* the 3 first properties are needed to properly define the FWHM (Full
* Half Width at Half Maximum).
*
* f(0) = 1 <==> alpha = 3/2
* f(x) = core(|x|) if if |x|<=1
* = wing(|x|) if 1<=|x|<=2
* with:
* core(x) = (beta*x*3/4 - 3/2)*(beta*x)^2 + 1
* wing(x) = (3/2)*(2 - beta*x)^3
* FWHM can be in the "core":
* core(1/2) = 1/2 <=> beta = (4 + 8*cos((arccos(-1/8) - 2*pi)/3))/3
* ~ 1.444703448928752590331004298512891328576
* in this case, beta*1/2 <= 1, so its is possible. FWHM can be in the
* "wing":
* wing(1/2) = 1/2 <=> beta = 4 - 2^(4/3)
* ~ 1.48015790021025367046557878544354329886
* in this case, beta*1/2 < 1 so it is not possible.
*
* Finally BETA = (4 + 8*cos((arccos(-1/8) - 2*pi)/3))/3
* ~ 1.444703448928752590331004298512891328576
*/
func __fit_spike_spline (&p, &dp, u) {
/* To have FWHM = 1, use:
* BETA = (4 + 8*cos((arccos(-1/8) - 2*pi)/3))/3
* ~ 1.444703448928752590331004298512891328576
*/
BETA = 1.444703448928752590331004298512891328576;
s = BETA*sign(u);
u *= s;
p = dp = array(double, dimsof(u));
if (is_array((i = where(u <= 1.0)))) {
ui = u(i);
p(i) = (0.75*ui - 1.5)*ui*ui + 1.0;
dp(i)= (2.25*ui - 3.0)*ui;
}
if (is_array((i = where((u > 1.0)*(u < 2.0))))) {
ui = 2.0 - u(i);
p(i) = 0.25*ui*ui*ui;
dp(i) = -0.75*ui*ui;
}
dp *= s;
}
func __fit_lp (&lmin, &lmax, &lgain)
{
if (is_void(lmin)) lmin = 1e-3;
else if (lmin <= 0.0 || lmin > 1.0) error, "bad value for keyword LMIN";
if (is_void(lmax)) lmax = 1e30;
else if (lmax <= 1.0) error, "bad value for keyword LMAX";
if (is_void(lgain)) lgain = 10.0;
else if (lgain <= 1.0) error, "bad value for keyword LGAIN";
return double(lmin);
}
func fit_2d_spike (z, y, x, &model, weight=, shape=,
guess=, bg=, absorption=, verb=,
maxeval=, maxiter=, conv=,
lmin=, lmax=, lgain=)
/* DOCUMENT
Parameters:
A(1) = peak intensity of spike
A(2) = FWHM of spike
A(3) = abscissa of spike
A(4) = ordinate of spike
A(5) = background level
A(6) = background X-slope
A(7) = background Y-slope
Model:
MODEL(X,Y) = A1*P((X-A3)/A2)*P((Y-A4)/A2) + A5 + A6*(X-A3) + A7*(Y-A4)
where P(U) is the shape function.
SHAPE is the 1-D shape of the spike, it can be "spline" to use a
cubic B-spline function, "gauss" to use a gaussian or any user
defined routine which has the following prototype:
func SHAPE(&p, &dp, u) {gg_fit_2d_spike(data)
p = f(u);
dp = f'(u);
}
where f(u) is a "bell-shaped" function and f'(u) is
its derivative; f(u) should have
with the following properties:
f(0) = 1
f(0.5) = 0.5
f(-u) = f(u) (i.e. f is symmetric)
f(u) >= 0 everywhere
f'(u)>0 for u<0, f'(0)=0, f'(u)<0 for u>0 (i.e. f is "bell-shaped")
the 3 first properties are needed to properly define the FWHM (Full
Half Width at Half Maximum).
Background:
BG = nil/0 no background (BG = 0 everywhere)
BG = 1 uniform background
BG = 2 BG(x,y) = a() + a()*(x-x0) + a()*(y-y0)
Spike:
Support: */
{
local g, h, chi2_best, a_best, delta_chi2, delta_a, p1, dp1, p2, dp2;
if (is_void(conv)) conv = 1e-5;
/* Initializes seetings for Levenberg-Marquardt parameter (smallness of
trust region). */
lp = __fit_lp(lmin, lmax, lgain);
if (! is_func(shape)) {
if (is_void(shape) || shape == "spline") shape = __fit_spike_spline;
else if (shape == "gauss") shape = __fit_spike_gauss;
else error, "bad value for keyword SHAPE";
}
if (! bg) nparams = 4;
else if (bg == 1) nparams = 5;
else if (bg == 2) nparams = 7;
else error, "bad value for keyword BG";
/* ONE is used to make conformable arrays the same dimension list as Z. */
one = array(1.0, dimsof(z));
/* Check WEIGHT, compute W=sqrt(WEIGHT). */
if (is_void(weight)) {
weight = 0n;
} else {
s = structof(weight);
if (s==double || s==float || s==int || s==long || s==short) {
/* Numerical weight. */
if (min(weight) < 0.0) error, "WEIGHT must be non-negative";
if (max(weight) <= 0.0) error, "WEIGHT must not be zero everywhere";
if (dimsof(weight)(1) != 0) {
/* Non-scalar weight. */
w = one*sqrt(weight);gg_fit_2d_spike(data)
weight = 1n;
} else {
/* Scalar weight: same as with no weight. */
weight = 0n;
}
} else error, "bad data type for WEIGHT";
}
if (is_void(guess)) {
/* Starting solution. PEAK is the list of indices in the half topmost
part of the spike (bottommost part for an absorption spike) and is
used to estimate the spike position and FWHM. */
a = array(double, nparams);
zmin = min(z);
zmax = max(z);
if (absorption) {
if (bg) {
a(5) = zmax;
a(1) = zmin - zmax;
peak = where(z <= 0.5*(zmax + zmin));
} else {
a(1) = zmin;
peak = where(z <= 0.5*zmin);
}
} else {
if (bg) {
a(5) = zmin;
a(1) = zmax - zmin;
peak = where(z >= 0.5*(zmax + zmin));
} else {
a(1) = zmax;
peak = where(z >= 0.5*zmax);
}
}
x_peak = (one*x)(peak);
y_peak = (one*y)(peak);
a(2) = 0.5*(abs(x_peak(ptp)) + abs(y_peak(ptp))); /* FWHM */
a(3) = avg(x_peak);
a(4) = avg(y_peak);
peak = x_peak = y_peak = [];
} else {
s = structof(guess);
if ((s==double || s==float || s==long || s==int || s==short || s==char) &&
dimsof(guess)(1)==1 && numberof(guess) == nparams) {
a = double(guess);
} else {
error, "bad GUESS parameters";
}
}
/* Minimum allowed FWHM (taken to be a small fraction of the sampling
to avoid division by zero). */
fwhm_min = 1e-5*sqrt(((max(x)-min(x))*(max(y)-min(y)))/double(numberof(z)));
//write, format="FWHM_min = %g\n", fwhm_min;
/* Model:
* m = a1*p((x-a3)/a2)*p((y-a4)/a2) + a5 + a6*(x-a3) + a7*(y-a4)
*
* Jacobian:
* jac1 = p((x-a3)/a2)*p((y-a4)/a2)
* jac2 = -(a1/a2)*((x-a3)/a2*p'((x-a3)/a2)*p((y-a4)/a2) +
* (y-a4)/a2*p((x-a3)/a2)*p'((y-a4)/a2))
* jac3 = -(a1/a2)*p((y-a4)/a2)*p'((x-a3)/a2) - a6
* jac4 = -(a1/a2)*p((x-a3)/a2)*p'((y-a4)/a2) - a7
* jac5 = 1
* jac6 = x-a3
* jac7 = y-a4
*
* Local quadratic approximation (l = likelihood, a = parameters,
* g = gradient vector with respect to parameters, H = hessian matrix,
* J = jacobian, i.e. paratial derivatives of model with repsect to
* parameters):
* l(a) = sum_x weight(x) * (model(a,x) - data(x))^2
* g(a) = 2 * sum_x weight(x) * (model(a,x) - data(x)) * J(a,x)
* H(a,b) ~ 2 * sum_x weight(x) * J(a,x) * J(b,x)
* ==> l(a+da) ~ l(a) + da'.g(a) + 1/2 da'.H.da
* ==> da_best = -LUsolve(H,g)
* = -LUsolve(1/2 H, 1/2 g)
* trust region control (lp = Levenberg-Marquardt parameter
* = smallness of trust region):
* da(lp) = -LUsolve(1/2 H(lp), 1/2 g)
* with:
* H(lp)(a,b) = H(a,b) where a!=b
* = (1 + lp)*H(a,a) where a=b
*/
warn = string(0);
iter = -1;
eval = 0;
for ( ; ; ) {
/* Apply constraints. */
if (a(2) < fwhm_min) a(2) = fwhm_min;
/* Compute model M, (weighted) residuals R and Chi2. */
s = 1.0/a(2); /* 1/FWHM */
jac6 = x - a(3);
jac7 = y - a(4);
shape, p1, dp1, (u1 = s*jac6);
shape, p2, dp2, (u2 = s*jac7);
if (nparams == 4) m = a(1)*p1*p2;
else if (nparams == 5) m = a(1)*p1*p2 + a(5);
else m = a(1)*p1*p2 + a(5) + a(6)*jac6 + a(7)*jac7;
r = m - z;
if (weight) r *= w;
chi2 = sum(r*r);
++eval;
if (! is_void(maxeval) && eval >= maxeval) {
warn = swrite(format="too many evaluations (%d)", eval);
}
if (eval <= 1 || chi2 < chi2_best) {
/* Solution improved: register current model and current solution as
the best one found so far. */
++iter;
eq_nocopy, model, m;
if (iter) {
delta_chi2 = abs(chi2_best - chi2);
delta_a = abs(a - a_best);
}
chi2_best = chi2;
a_best = a;
if (verb && iter%verb == 0) {
if (eval==1) {
write, format="%s\n%s\n",
"ITER EVAL CHI2 FWHM X0 Y0",
"---- ---- ----------- --------- ------------ ------------";
}
write, format="%4d %5d %11.5e %-9.3g %-+12.5g %-+12.5g\n",
iter, eval, chi2, a(2), a(3), a(4);
}
/* Check for termination. */
if (chi2 <= 0.0 || (iter && delta_chi2 <= conv*abs(chi2))) {
warn = string(0);
break;
}
if (! is_void(maxiter) && iter >= maxiter) {
warn = swrite(format="too many iterations (%d)", iter);
break;
}
if (warn) break; /* too many evaluations */
/* Solution improved: augment trust region. */
lp = max(lp/lgain, lmin);
/* Compute Jacobian:
* jac1 = p1*p2
* jac2 = -(a1/a2)*(u1*dp1*p2 + u2*p1*dp2)
* jac3 = -(a1/a2)*p2*dp1 - a6
* jac4 = -(a1/a2)*p1*dp2 - a7
* jac5 = 1
* jac6 = x-a3
* jac7 = y-a4
*/
jac = array(pointer, nparams);
q = -a(1)/a(2);
if (weight) {
jac(1) = &(w*p1*p2);
jac(2) = &(w*q*(u1*dp1*p2 + u2*p1*dp2));
if (nparams >= 5) jac(5) = &w;
if (nparams == 7) {
jac(3) = &(w*(q*p2*dp1 - a(6)));
jac(4) = &(w*(q*p1*dp2 - a(7)));
jac(6) = &(w*jac6);
jac(7) = &(w*jac7);
} else {
jac(3) = &(w*q*p2*dp1);
jac(4) = &(w*q*p1*dp2);
}
} else {
jac(1) = &(p1*p2);
jac(2) = &(q*(u1*dp1*p2 + u2*p1*dp2));
if (nparams >= 5) jac(5) = &one;
if (nparams == 7) {
jac(3) = &(q*p2*dp1 - a(6));
jac(4) = &(q*p1*dp2 - a(7));
jac(6) = &jac6;
jac(7) = &jac7;
} else {
jac(3) = &(q*p2*dp1);
jac(4) = &(q*p1*dp2);
}
}
/* Compute (half of) Hessian and (half of) gradient to obtain local
quadratic approximation. */
h = array(double, nparams, nparams);
g = array(double, nparams);
for (i=1 ; i<=nparams ; ++i) {
if (weight) {
jac_i = w*(*jac(i));
} else {
local jac_i;
eq_nocopy, jac_i, *jac(i);
}
for (j=1 ; j<i ; ++j) h(i,j) = (h(j,i) = sum(jac_i * *jac(j)));
h(i,i) = sum(jac_i * jac_i);
g(i) = sum(jac_i * r);
}
/* Account for bound constraints on the FWHM. */
if (a(2) <= fwhm_min && g(2) > 0.0) g(2) = 0.0;
} else {
/* Solution not improved, reduce trust region. */
if (lp >= lmax)
warn = "Levenberg-Marquardt parameter overflow (convergence)";
if (warn) break; /* parameter overflow or too many evaluations */
lp = min(lp*lgain, lmax);
}
/* Solve the regularized local quadratic problem. */
tmp = h;
tmp(::nparams+1) *= 1.0 + lp;
a = a_best - LUsolve(tmp, g);
}
if (warn && verb!=0) write, format="WARNING - %s\n", warn;
return a;
}
|