 # document.write(pagetitle) Home Download Tutorial 1 Tutorial 2 Tutorial 3 Tutorial 4 Tutorial 5
```;+
; NAME:
;               ULY_FIT_LIN
; PURPOSE:
;               Fit the linear coefficients and return the best fit
;
; USAGE:
;                      PAR_LOSVD=par_losvd, CMP=cmp,
;                      SPECTRUM=SignalLog,
;                      GOODPIXELS=goodPixels,
;                      VOFF=voff,
;                      POLPEN=polpen [,/LUM_WEIGHT] [,/QUIET])
;
; DESCRIPTION:
;    Evaluate the different components of the model to fit and convolve
;    them with the LOSVD for a given set of the free parameters.
;    Determine the optimal weight of each component and the coefficients
;    of the multiplicative and additive continuum (linear fit).
;    Return the best fitting model.
;
;    This function is called at each step of the Levenberg-Marquart
;    minimization performed in ULY_FIT to compute the values of the
;    non-linear parameters of the components and LOSVD.
;
;    For performance reasons, the large data arrays are passed through
;    commons (not as arguments).
;
; KEYWORDS:
;   POLPEN (input)
;      Bias level of the multiplicative polynomial. This keyword can be
;      used to reduce the effect of insignificant terms of the multiplicative
;      polynomial.
;      Terms whose coefficients are smaller, in absolute value, than <polpen>
;      times their statistical error are amorted by the factor
;      (abs(coef)/(<polpen>*err))^2.
;      This feature is active only if MDEGREE>0. <polpen>=2 seems a
;      reasonable choice.
;      By default,  no bias is applied.
;   PAR_LOSVD (input)
;      Parameters of the LOSVD, array of 0 to 6 numbers ordered as:
;      [cz, sigma, h3, h4, h5, h6], cz and sigma are in pixels.
;   SPECTRUM  (input)
;      Structure containing the spectrum to analyse and its associated error.
;      It is used to determine the weight of each pixel (.err) and to
;      fit the multiplicative continuum (.data)
;   GOODPIXELS (input)
;      List of pixels to be used in the fit.
;   VOFF (input)
;      Velocity offset between the spectrum to analyse and the model. km/s.
;   CMP       (Input & output)
;      Array of components to fit (see documentation in ULY_FIT)
;      In output CMP.WEIGHT is updated with the weight of each component
;      computed by this routine.
;   MPOLY  (Input/output)
;      multiplicative continuum Legendre polynomial
;      structure {lmdegree, mpolcoefs, poly}
;      .LMDEGREE (integer), Input, maximum degree of legendre polynomials.
;         Unchanged in output.
;      .MPOLCOEFS (double 1D array), LMDEGREE+1 terms.
;         Output: updated values of the coefficients
;         The input values of the coefficients are ignored.
;      .POLY (1D array)
;         Input: polynomial determined at a previous iterration
;         Output: updated value of the polynomial (consistent with .MPOLCOEFS)
;      Usually MPOLY is initialized by ULY_FIT, but if it is NULL in input,
;      it is initialized to contain a degree 0 polynomial (scaling factor).
;   /LUM_WEIGHT
;      Instruct to compute the luminosity weights and the errors on the
;      weights.
;   /QUIET
;      Set this keyword to run silently.
;
; AUTHOR:
;               Mina Koleva, Philippe Prugniel
;               with credit to Michele Cappellari & Igor Chilingarian
;               for initial development
; HISTORY:
;    The principle of solving the linear dependences within the Levenberg-Marquart
;    function evaluation is borrowed from M. Cappellari in his ppxf package.
;    v.1.0      From early development by M. Cappelleri & I. Chilingarian
;    v.2.0      Separated from the main program, Mina Koleva Sep 2005
;    v.2.3      Changing the interpolation (improve performance): mk 2007/03/01
;    v.3.0      Interpolating between different Mg/Fe: mk 2007/04/17
;    v.4.0      Ph. Prugniel 2007/12/14: linear fit of mult pol. Cache Legendre
;    v.5.0      Ph. Prugniel, 2008/05/27: general n-components fit
;    v.6        Ph. Prugniel, 2011/02: Implement Tikhonov regularisation
;    v.7        Ph. Prugniel, 2015/09: Major debugging and improvements
;-
; CATEGORY:    ULY
;------------------------------------------------------------------------------
;
; uly_fit_fractsolve:
;  This routine determines the weight of the different 'components' of the
;  model: one or several components and 'npoly' terms of the additive
;  polynomial.
;  'a' can be a vector (if there is a single template and no polynomial)
;  or an array whose first 'npoly' lines are the polynomial basis, and
;  the next line(s) the templates.
;  'b' is the 'observed' spectrum
;  bvls_state is a I/O array used to store the state of BVLS_PS, for
;  the warm start
;
;  The function returns an array containing the weight of each component.
;  It is called by 'uly_fit_lin_weight'.

function uly_fit_fractsolve, a, b, npoly, BOUNDS=bounds, STATE=bvls_state
compile_opt idl2, hidden
s = size(a)

if s eq 1 then begin                   ; Fit a single cmp, no add polynomial
soluz = total(a*b,/DOUBLE)/total(a^2,/DOUBLE)
endif else if s eq npoly+1 then begin ; Fitting a single cmp (no bound)
soluz = la_least_squares(transpose(a), b, /DOUBLE)
endif else begin                         ; Fitting multiple templates (bounds)
bnd = replicate(1d,2,s) ; faster than dblarr
;    bnd = dblarr(2,s, /NOZERO)
mx = (machar(/DOUBLE)).xmax
if npoly gt 0 then begin
bnd[0,0:npoly-1] = -mx ; No bounds on additive polynomials
bnd[1,0:npoly-1] = mx  ; No bounds on additive polynomials
endif
bnd[*,npoly:*] = double(bounds)
bvls_ps, a, b, bnd, soluz, bvls_state, ITMAX=15*s, STATUS=ierr
;    bvls_lh, a, b, bnd, soluz, ITMAX=15*s, STATUS=ierr
if ierr ne 0 then message, 'BVLS returned with error #' + strtrim(ierr,2)
endelse

return, soluz
end

;-----------------------------------------------------------------------------
; Determination (or refining) of the multiplicative Legendre polynomial
function uly_fit_lin_mulpol, bestfit, CMP=cmp, MPOLY=mpoly,           \$
SPECTRUM=SignalLog,                      \$
GOODPIXELS=goodPixels,                   \$
QUIET=quiet,                             \$
STATUS=status

status = 0
npix = (size(*SignalLog.data, /DIM))        ; number of wavelength bins
galaxy = double((*SignalLog.data)[goodPixels]) ; signal for non-masked points
noise = double((*SignalLog.err)[goodPixels])   ; error for non-masked points

if mpoly.lmdegree eq 0 then return, bestfit

; The polynomial may get crazy in the unconstrained regions, and this
; may affect the edges of these regions...
; It is critical when we want to fit spectral segment separated by a large
; range without any data.
; To prevent this, we may fit the polynomial everywhere, replacing the
; missing data with values in a correct range, and giving them low weight.
; like:
;    galtmp = galaxy*0 + 1e
;    galtmp[goodPixels] = galaxy[goodPixels]
;    meas_err = galaxy*0 + (mean(noise[goodPixels]) * 100)
;    meas_err[goodPixels] = noise[goodPixels]
; and then fit as:
;    coefs_pol = mregress(cmul, galtmp, MEASURE_ERRORS=meas_err, INVERSE=inverse, STATUS=status)
; This should work in most places, and replacing the '1' by an
; interpolation of the data should make this solution quite general.
; Though, for the moment, we prefer to keep the conservative solution.
; So, the problem still exists.

bestfitp = bestfit / mpoly.poly
cmul = mpoly.leg_array * rebin(bestfitp,n_elements(bestfitp),mpoly.lmdegree+1)

coefs_pol = mregress(cmul[goodPixels,*], galaxy, MEASURE_ERRORS=noise, STATUS=status_mregress)
if status_mregress eq 1 then status = 1
; note that la_least_squares or bvls_ps are about 50% slower than mregress
;  coefs_pol = la_least_squares(transpose(cmul[goodPixels,*]), galaxy, /DOUBLE, STATUS=status)
;  coefs_pol = transpose(coefs_pol)

if n_elements(polpen) ne 0 then begin
;  penalization of unsignificant terms in the pol.
penal_pol = abs(coefs_pol)/(polpen*inverse.sigma)
pen = 1+[where (penal_pol[1:*] lt 1)]
coefs_pol[pen] *= penal_pol[pen]^2
endif

mpoly.poly = mpoly.leg_array[*,0:mpoly.lmdegree] # transpose(coefs_pol)

minmaxpol = minmax(mpoly.poly[goodPixels])
if minmaxpol lt 0 then begin
lmd = mpoly.lmdegree
while min(mpoly.poly) le 0 and lmd gt 0 do begin
lmd -= 1
coefp = mregress(cmul[goodPixels,0:lmd], galaxy, \$
MEASURE_ERRORS=noise,       \$
STATUS=status_mregress)
if status_mregress eq 1 then status = 1
if n_elements(polpen) ne 0 then begin
penal_pol = abs(coefp) / (polpen*inv.sigma)
pen = 1+[where (penal_pol[1:*] lt 1)]
coefs_pol[pen] *= penal_pol[pen]^2
endif
mpoly.poly = (mpoly.leg_array[*,0:lmd] # transpose(coefp))
endwhile

if not keyword_set(quiet) then \$
print, 'WARNING: the polynomial became negative, min:'+strtrim(minmaxpol,2),' its degree is reduced to:', lmd

coefs_pol = 0*coefs_pol + coefp

if min(mpoly.poly) le 0 then if not keyword_set(quiet) then begin
message, /INFO, 'polynomial is negative'
status = 1
return, 0
endif
endif

bestfitp *= mpoly.poly
if keyword_set(lum_weight) and n_elements(cmp.weight) gt 0 then begin
cmp.e_weight /= coefs_pol
cmp.l_weight /= coefs_pol
endif
mpoly.mpolcoefs = coefs_pol

return, bestfitp
end

;------------------------------------------------------------------------------
; Determination of the weight of each component
function uly_fit_lin_weight, SignalLog, goodpixels, tmp, mpoly, \$
LUM_WEIGHT=lum_weight,             \$
CMP=cmp,                           \$
STATE=bvls_state,                  \$
QUIET=quiet,                       \$
STATUS=status
; Input arguments:
;   SignalLog  : observation (spect structure)
;   goodpixels : list of pixels to use
;   tmp        : templates
;   mpoly      : structure containing the multiplicative polynomial
; Keywords (input)
;   lum_weight : keyword switch: if set, compute e_weight and l_weight
; Keyword (output)
;   cmp        : result returned as cmp.weight (and .e_weight .l_weight if /LUM_WEIGHT)
;   status     : error status: 0= no error
; RETURN
;   bestfit    : array with the best matching model

n_cmp = n_elements(cmp)                        ; number of components
npix = (size(*SignalLog.data, /DIM))        ; number of wavelength bins
n_cmp = n_elements(cmp)                        ; number of components
galaxy = double((*SignalLog.data)[goodPixels]) ; signal for non-masked points
noise = double((*SignalLog.err)[goodPixels])   ; error for non-masked points
status = 0

; Fill the columns of the design matrix of the least-squares problem
c = rebin(double(mpoly.poly), npix, n_cmp) * tmp  ; Fill with models*mulpol

x = 2d*dindgen(npix)/npix - 1d
c1[0,j] *= legendre(x,j)
endfor
c = [[c1],[c]]
endif

a = c[goodpixels,*]                            ; Used to solve the system
for j=0,(adegree+1)+n_cmp-1 do a[*,j] /= noise ; Weight with errors

;-----------------------------------------------------------------------------
; Experiment with Tikhonov regularization (2011/02/05)
; We add to the system above some linear constraints , weighted with "lambda"
; See the chapters 18.4 and 18.5 of Numerical Recipes
lambda = 0 ; regularization weight : 0 => no regularization
;lambda = 1d-13 * double(npix)/( n_cmp-1)
;
; For the moment, we experimented with 1st order regularization, i.e.
; the weights are equal.
; The solution that we tried should still be finalized: The regularization
; weight should be scaled more intelligently. lambda should be insensitive
; to: (i) the number of wavelength bins (or VELSCALE), (ii) the flux scale
; in the data and (iii) the flux scale in the models. 5This should not be
; an issue to implement this
; We shall also explore these other aspects:
;   - Other types of regularization, like 2nd order (linear dependence)
;     or regularization matrix passed by the user (we may provide a different
;     routine to help the creation of this matrix)
;   - We may want to regularize on the l_weight rather than on the weights.
; We have to find a compromize between the flexibility of the program and
; its ease of use. The regularization is a complex business: How can the
; user choose the proper value of lambda, and the type of regularization?
; If we just give a single regularization (say 1st order) and a default
; lambda ... the user may be happy to have something functional at low
; expense of brain. But it is not fully correct.
if lambda gt 0 then begin;
ng = n_elements(goodpixels)
nrel = ng + n_cmp - 1
aa = dblarr(nrel, n_cmp)
for k = 0, n_cmp-1 do aa[0,k] = a[*, k]
a = temporary(aa)
cr = dblarr((size(c, /DIM)))
for k=ng, nrel-1 do begin
a[k,*] = shift(cr, k-ng)
endfor
goodpixels = [goodpixels, lindgen(n_cmp-1)+npix]
galaxy = [galaxy,fltarr(n_cmp-1)]
noise = [noise,fltarr(n_cmp-1)+1]
endif

status_math = check_math(MASK=32)    ; galaxy/noise may cause an underflow: handle
cmp.weight = uly_fit_fractsolve(a, galaxy/noise, npoly, BOUNDS=cmp.lim_weig, STATE=bvls_state)
bestfit = c # cmp.weight
endif else begin
wght = uly_fit_fractsolve(a, galaxy/noise, npoly, BOUNDS=cmp.lim_weig, STATE=bvls_state)
bestfit = c # wght
endelse
if status_math eq 0 then status_math = check_math(MASK=32)

nan = where(finite(cmp.weight) eq 0 or abs(cmp.weight) gt 1d36, cnt)
if cnt gt 0 then begin
message, /INFO, 'Could not determine the weight of cmp ' + \$
strjoin(strtrim(nan,2),',')+' ... assume 0'
cmp[nan].weight = 0
endif
zeros = where(abs(cmp.weight) lt 1d-36, cnt)
if cnt gt 0 then cmp[zeros].weight = 0d
if abs(total(cmp.weight)) lt 1d-36 then begin
if not keyword_set(quiet) then \$
message, /INFO, 'The total weight is null'
status = 1
common MPFIT_ERROR, error_code
error_code = -2 ; tell MPFIT to stop, the best fit is a null spectrum
return, 0
endif

if keyword_set(lum_weight) and n_elements(cmp.weight) gt 0 then begin
;   we call this only when we want to evaluate the errors on the
;   weights and the luminosity weights
w_notnull = where(cmp.weight ne 0, nw_notnull)
cmp.e_weight = dblarr(n_cmp)
cmp.l_weight = dblarr(n_cmp)
if nw_notnull ne 0 then begin
inv = 0
if adegree gt -1 then \$
else \$
w_nn =  w_notnull
ww = mregress((c[goodPixels,*])[*,w_nn], galaxy, \$
MEASURE_ERRORS=noise, INVERSE=inv, STATUS=status)
cmp.l_weight = (total(tmp, 1) * cmp.weight) / npix
endif
endif else begin
cmp.e_weight = 0
cmp.l_weight = 1
endelse

return, bestfit

end
;------------------------------------------------------------------------------
PAR_LOSVD=par_losvd,        \$
CMP=cmp,                    \$
SPECTRUM=SignalLog,         \$
GOODPIXELS=goodPixels,      \$
VOFF=voff,                  \$
MPOLY=mpoly,                \$
LUM_WEIGHT=lum_weight,      \$
POLPEN=polpen,              \$
CACHE=cache,                \$
STATE=bvls_state,           \$
MODECVG=modecvg,            \$
QUIET=quiet,                \$
STATUS=status

compile_opt idl2

status = 0

; Clear the accumulated stack of math error, in order to diagnos errors
; from the current routine.
if !EXCEPT gt 0 then status_math = check_math(/PRINT) else \$
if check_math() then print, \$
'An arithmetic error occured before calling uly_fit_lin'

npix = (size(*SignalLog.data, /DIM)) ; number of wavelength bins
n_cmp = n_elements(cmp)                 ; number of components
npar_losvd = n_elements(par_losvd)      ; number of free params for LOSVD

galaxy = double((*SignalLog.data)[goodPixels]) ; signal for non-masked points
noise = double((*SignalLog.err)[goodPixels])   ; error for non-masked points

if n_elements(mpoly) eq 0 then begin ; default poly initialization
mpoly = {lmdegree:0, mpolcoefs:dblarr(1), poly:dblarr(npix)}
mpoly.poly = 1d
endif

; check that noise is strictly positive, at least on the place of goodpixels
if min(noise) le 0 then message,  'NOISE vector must be strictly positive'

;-----------------------------------------------------------------------------
; evaluate all the components
if size(cache,/TYPE) ne 8 then begin
cache = {para:ptrarr(n_elements(cmp), /ALLOCATE_HEAP), spectra:dblarr(npix, n_cmp)}
for i=0, n_elements(cmp)-1 do *(cache.para)[i] = -1D35
endif

;init. models spectra
models = dblarr(npix,n_cmp,/NOZERO)

for i=0, n_cmp-1 do begin
cached = 0
if n_elements(*cmp[i].para) gt 0 then begin
if min((*cmp[i].para).value eq *(cache.para)[i]) eq 1 then begin
models[0,i] = cache.spectra[*,i]
cached = 1
endif
endif

if cached eq 0 then begin
;       We treat separately SSPs and stars because it is (a bit) faster
;       than doing everything with call_function
case cmp[i].eval_fun of
'SSP' : mdl = \$
uly_ssp_interp(*cmp[i].eval_data, (*cmp[i].para).value)
'STAR' : mdl = *cmp[i].eval_data
else : mdl = \$
call_function(cmp[i].eval_fun, *cmp[i].eval_data, \$
(*cmp[i].para).value)
endcase
nmd = min([npix, (size(mdl,/DIM))]) - 1
models[0:nmd,i] = mdl[0:nmd]
if n_elements(*cmp[i].para) gt 0 then \$
*(cache.para)[i] = (*cmp[i].para).value
cache.spectra[*,i] = models[*,i]
endif

endfor

if n_elements(par_losvd) gt 0 then begin
;-----------------------------------------------------------------------------
; generate the LOSVD kernel and convolve
; Detailed tests have shown that this kernel is correct as long
; as sigma > 0.8 (it is still acceptable down to sig=0.6)
; For low values of sigma:
;   - The total flux of the kernel is too large (need a normalization)
;   - The function is too much picked (for sigma around 0.3)
;     when the kernel is centered and too flat when the kernel is
;     shifted by half a pixel
; This low-sigma bias of the kernel can be modelled as:
;     sigmac = sqrt(sigma^2 + 0.125*exp(-(sigma-0.28)^2/(0.04)))
; Where sigma is the value used in the program and sigmac the de-biased
; value corresponding to the real broadening. For the lack of simplicity,
; this correction is not applied in the program, but the user should
; be aware of the bias. A healthy solution consist in rebinning the
; spectrum to a scale where sigma>0.6

maxvel = abs(par_losvd)
maxsig = par_losvd
dx = ceil(abs(voff)+abs(maxvel)+5d*maxsig) ; Sample the LOSVD at least to vel+5*sigma
dx = min([dx, (npix-1.)/2.])
n = 2*dx + 1
x = dx - findgen(n)
losvd = dblarr(n,/NOZERO)

vel = voff + par_losvd
sigma_pix = (par_losvd gt 0.1) ? par_losvd : 0.1d
w = (x - vel)/sigma_pix
w2 = w^2
;   Replace the values larger then 5*sigma with zeros
wlarge = where(abs(w) gt 5., cnt, COMPLEMENT=wnorm)
if cnt gt 0 then losvd[wlarge] = 0.
losvd[wnorm] = exp(-0.5d*w2[wnorm])/(sqrt(2d*!dpi)*sigma_pix)

; Hermite polynomials normalized as in Appendix A of van der Marel &
; Franx (1993). They are given e.g. in Appendix C of Cappellari et al. (2002)
nherm = (size(par_losvd,/DIM)) - 2
if nherm gt 0 then begin
poly = 1d + par_losvd/sqrt(3d)*(w*(2d*w2-3d)) ; h3
if nherm gt 1 then \$
poly += par_losvd/sqrt(24d)*(w2*(4d*w2-12d)+3d) ; h4
if nherm gt 2 then \$
poly += par_losvd/sqrt(60d)*(w*(w2*(4d*w2-20d)+15d)) ; h5
if nherm gt 3 then \$    ; h6
poly += par_losvd/sqrt(720d)*(w2*(w2*(8d*w2-60d)+90d)-15d)

losvd *= poly
endif
losvd /= total(losvd) ; Normalized total(Gaussian)=1 (need for sigma<0.6)

;   convolution with the LOSVD
status_math = check_math(MASK=32) ; convolution may cause an underflow: handle
tmp = convol(temporary(models),losvd,/EDGE_TRUNCATE)
if status_math eq 0 then status_math = check_math(MASK=32)

endif else tmp = temporary(models)

; Iterative loop to get the weight and mulcont
;   when n_cmp>1 and md>0 the determination of the weights
;   and of mulpol is not a linear process, and at present we
;   need iteration to converge... and unfortunately, our
;   current approach is slow (we tried to accelarete it but
;   we failed for the moment to obtain a stable solution).
;   However, in most cases, the LM iteration will converge
;   even if uly_fit_lin does not converge at each step.
; So, we are currently offering three options:
;   MODECVG=0 (the default one), this is the fastest, but
;     it happens that the solution is missed
;   MODECVG=1 (the convergence is achieved only once for each
;     LM iteration, it is not reached for the derivatives
;   MODECVG=2, uly_fit_lin always converge, but this is slow
; We are actively working at solving the problem and make
; something fast and reliable.
itermax = 0
if n_elements(modecvg) eq 0 then modecvg=0
if modecvg eq 1 or modecvg eq 2 then itermax=499

fstop = 0
niter = 0
while fstop eq 0 do begin
;  Determination of the weight of each component
bestfitw = uly_fit_lin_weight(SignalLog, goodpixels, tmp, mpoly, \$
LUM_WEIGHT=lum_weight,             \$
CMP=cmp,                           \$
STATE=bvls_state,                  \$
QUIET=quiet,                       \$
STATUS=status)
if status ne 0 then return, bestfitw

posw = where(cmp.weight gt 0, cposw)
if cposw eq 0 then begin
if not keyword_set(quiet) then \$
message, /INFO, 'All the cmp have a null weight... cannot continue'
common MPFIT_ERROR, error_code
error_code = -2          ; tell MPFIT to stop
status = 1
return, 0
endif

; Determination of the multiplicative polynomials
bestfit = uly_fit_lin_mulpol(bestfitw, CMP=cmp, MPOLY=mpoly,    \$
SPECTRUM=SignalLog,                \$
GOODPIXELS=goodPixels,             \$
QUIET=quiet,                       \$
STATUS=status)
if status ne 0 then return, float(bestfit)

mpoly.poly /= mpoly.mpolcoefs
cmp.weight *= mpoly.mpolcoefs
mpoly.mpolcoefs /= mpoly.mpolcoefs
if cposw eq 1 then break ; no need to iterate when only on cmp is active

niter += 1
if abs(1d0-total(bestfit,/NAN)/total(bestfitw,/NAN)) lt 5d-9 or niter gt itermax then fstop = 1
;   if abs(mean(1d0-bestfit/bestfitw)) lt 5d-9 or niter gt itermax then fstop = 1
endwhile

bestfit = float(bestfit)

; The test below is useful to identify that an error occured during
; the execution of uly_fit_lin.
; Note that we cannot use a CATCH handler, because it does not catch
; math errors !!
; Note that to find the actual statement causing an error, the
; variable !EXCEPT can be set to 3 in the IDL session.
if check_math(/PRINT) ne 0 then begin
message, /INFO, 'an arithmetic error occured'
common MPFIT_ERROR, error_code
error_code = -2 ; tell MPFIT to stop
status = 1
return, 0
endif

return, bestfit
end
```
Part of ULySS package