Viewing contents of file '../idllib/ssw/allpro/amoeba_c.pro'
PRO AMOEBA_C,X,Y,FNAME,PARAM,ACCURACY=ACCURACY,MAX_ITER=MAX_ITER, $
POISSON=POISSON,ERROR=ERR,LAMBDA=LAMBDA0,CHISQR=CHISQR, $
N_ITER=ITER,NOPRINT=NOPRINT,ABSOLUTE=ABSOLUTE,LORENTZ=LORENTZ,$
PRANGE=K_PRANGE
;+
; Project : SOHO - CDS
;
; Name : AMOEBA_C
;
; Purpose : Reiteratively fits an arbitrary function
;
; Explanation : Fits an arbitrary function to a series of data points via a
; least-squares reiterative technique.
;
; The procedure used is taken from Numerical Recipes.
;
; Use : AMOEBA_C, X, Y, FNAME, PARAM
;
; Inputs : X = Positions.
; Y = Data values.
; FNAME = Name of function to be fitted (string variable).
; PARAM = Parameters of fit. Passed as first guess. Returned
; as fitted values.
;
; Opt. Inputs : None.
;
; Outputs : PARAM = Parameters of fit. See note above.
;
; Opt. Outputs: None.
;
; Keywords :
; ACCURACY = Accuracy to cut off at. Defaults to 1E-5.
; MAX_ITER = Maximum number of reiterations. Defaults to 20.
; POISSON = If set, then a Poisson error distribution is assumed, and
; the weights are set accordingly to 1/Y.
; ERROR = Array of errors. The weights are set accordingly to
; 1/ERROR^2 (normal distribution). Overrides POISSON.
; LAMBDA = Initial step sizes for PARAM, or if scalar then fraction of
; PARAM. Defaults to 1E-2. When passed as an array, this
; parameter can be used to hold parameters constant by setting
; LAMBDA(I)=0 for those parameters.
; NOPRINT = If set, then no printout is generated.
; CHISQR = Returned value of chi-squared. Only relevant if ERROR
; passed explicitly.
; N_ITER = Number of iterations used.
; ABSOLUTE = If set, then the sum of the absolute differences is
; minimized instead of the sum of the squares. This is
; equivalent to assuming a double-sided exponential
; distribution.
; LORENTZ = If set, then a Lorentz distribution is used instead of a
; normal distribution. Not truely meaningful unless ERROR is
; passed.
; PRANGE = Range of acceptable parameter values. Must have the
; dimensions (NPAR,2), where PRANGE(*,0) are the minimum
; values and PRANGE(*,1) are the maximum values. Only those
; ranges where PRANGE(*,1) are larger than PRANGE(*,0) are
; implemented--all other parameters are considered to be
; unbounded.
;
; Calls : FORM_CHISQR, FORM_SIGMAS
;
; Common : None
;
; Restrictions: The user defined function is passed by name as a character
; string in the variable FNAME. The function must have the form.
;
; Y = F(X,PARAM)
;
; where X is the independent variable, and PARAM is the vector
; containing the parameters of the fit.
;
; Unless LAMBDA is passed as an array, if the initial guess for
; PARAM contains any zeroes, then those parameters will be kept
; constant at zero.
;
; Side effects: The statistical interpretation of CHISQR is unclear if either
; ABSOLUTE or LORENTZ is set.
;
; Category : Utilities, Curve_Fitting
;
; Prev. Hist. :
; William Thompson, August, 1989.
; William Thompson, June 1991, modified to use keywords.
; William Thompson, December 1992, modified to use other fitting
; strategies besides minimizing the root-mean-square. Also
; removed keyword WEIGHTS as this was incompatible with this
; strategy.
;
; Written : William Thompson, GSFC, August 1989
;
; Modified : Version 1, William Thompson, GSFC, 9 January 1995
; Incorporated into CDS library
; Version 2, William Thompson, GSFC, 6 October 1995
; Fixed typo.
; Version 3, William Thompson, GSFC, 10 December 1997
; Renamed to AMOEBA_C to avoid conflict with IDL/v5
; routine of the same name.
; Version 4, William Thompson, GSFC, 05-Jun-1998
; Fixed bug where LAMBDA was not being properly accepted
; if a vector.
; Version 5, William Thompson, GSFC, 27-Apr-1999
; Added keyword PRANGE.
; Version 6, William Thompson, GSFC, 04-May-1999
; Allow for case where either PARAM or LAMBDA is set to
; zero to keep a parameter constant.
;
; Version : Version 6, 04-May-1999
;-
;
ON_ERROR,2
;
; Check the number of parameters passed.
;
IF N_PARAMS(0) NE 4 THEN BEGIN
PRINT,' This procedure must be called with four parameters:'
PRINT,' X, Y, FNAME, PARAMS'
RETURN
ENDIF
;
; Set up the default parameters.
;
ACCUR = 1E-5
IF N_ELEMENTS(ACCURACY) EQ 1 THEN $
IF ACCURACY GT 0 THEN ACCUR = ACCURACY
;
IF N_ELEMENTS(MAX_ITER) EQ 1 THEN NMAX = MAX_ITER ELSE NMAX = 20
;
LAMBDA = 1E-2
IF N_ELEMENTS(LAMBDA0) EQ N_ELEMENTS(PARAM) THEN LAMBDA = LAMBDA0
IF N_ELEMENTS(LAMBDA0) EQ 1 THEN $
IF LAMBDA0 GT 0 THEN LAMBDA = LAMBDA0
;
; Define NDATA and NPAR from the input arrays. Start ITER at zero.
;
NDATA = N_ELEMENTS(X) < N_ELEMENTS(Y)
NPAR = N_ELEMENTS(PARAM)
IF NPAR EQ 1 THEN PARAM = MAKE_ARRAY(PARAM)
ITER = 0
;
; Calculate the sigmas to use for each data point.
;
FORM_SIGMAS, Y, SIG, POISSON=POISSON, ERROR=ERR
;
; Check the array or scalar LAMBDA.
;
NLAMB = N_ELEMENTS(LAMBDA)
IF (NLAMB NE NPAR) AND (NLAMB NE 1) THEN BEGIN
PRINT,'*** LAMBDA must have 1 or ' + FIX(NPAR) + $
' elements, routine AMOEBA_C.'
RETURN
ENDIF
;
; Calculate the array of parameters to start with.
;
P = PARAM # REPLICATE(1.,NPAR+1)
FOR I = 0,NPAR-1 DO BEGIN
IF NLAMB EQ NPAR THEN P(I,I+1) = PARAM(I) + LAMBDA(I) $
ELSE P(I,I+1) = PARAM(I) * (1. + LAMBDA)
ENDFOR
;
; Determine the allowed parameter range, and make sure that all the parameters
; are within that range.
;
PR_COUNT = 0
SZ = SIZE(K_PRANGE)
IF SZ(0) EQ 2 THEN IF (SZ(1) EQ NPAR) AND (SZ(2) EQ 2) THEN BEGIN
PRANGE=K_PRANGE
W_PR = WHERE(PRANGE(*,1) GT PRANGE(*,0), PR_COUNT)
IF PR_COUNT GT 0 THEN FOR I = 0,PR_COUNT-1 DO BEGIN
J = W_PR(I)
P(J,*) = PRANGE(J,0) > P(J,*) < PRANGE(J,1)
ENDFOR
ENDIF
;
; Initialize the array containing chi-squared.
;
CHI = 0.*[0,PARAM]
FOR I = 0,NPAR DO BEGIN
TEST = EXECUTE('F = ' + FNAME + '(X,P(*,I))')
CHI(I) = FORM_CHISQR( (F - Y)/SIG, ABSOLUTE=ABSOLUTE, $
LORENTZ=LORENTZ)
ENDFOR
;
; Print the header.
;
IF (NMAX GT 0) AND (NOT KEYWORD_SET(NOPRINT)) THEN BEGIN
PRINT,FORMAT="(1H1,3X,'I',10X,'log',20X,'log',16X,'I')"
PRINT,FORMAT="(1X,I4,9X,'ERROR',17X,'CHISQR',7X,I9)",NMAX,NMAX
PRINT,FORMAT="(2X,4(1H-),5X,10(1H-),5X,25(1H-),4X,4(1H-))"
ENDIF
;
; Starting point for all iterations. Find the highest, next highest, and
; lowest values of CHI.
;
ITERATE:
ITER = ITER + 1
ILO = 0
IF CHI(0) GT CHI(1) THEN BEGIN
IHI = 0
INHI = 1
END ELSE BEGIN
IHI = 1
INHI = 0
ENDELSE
FOR I = 0,NPAR DO BEGIN
IF CHI(I) LT CHI(ILO) THEN ILO = I
IF CHI(I) GT CHI(IHI) THEN BEGIN
INHI = IHI
IHI = I
END ELSE IF (CHI(I) GT CHI(INHI)) AND (I NE IHI) THEN INHI = I
ENDFOR
;
; Find the average of P for all points other than IHI.
;
PBAR = 0*PARAM
FOR I = 0,NPAR DO IF I NE IHI THEN PBAR = PBAR + P(*,I)
PBAR = PBAR / NPAR
;
; Reflect from the high point through PBAR.
;
PR = 2*PBAR - P(*,IHI)
IF PR_COUNT GT 0 THEN PR(W_PR) = $
PRANGE(W_PR,0) > PR(W_PR) < PRANGE(W_PR,1)
TEST = EXECUTE('F = ' + FNAME + '(X,PR)')
CR = FORM_CHISQR( (F - Y)/SIG, ABSOLUTE=ABSOLUTE, LORENTZ=LORENTZ)
;
; If an improvement was achieved, try an additional extrapolation.
;
IF CR LE CHI(ILO) THEN BEGIN
PRR = 2*PR - PBAR
IF PR_COUNT GT 0 THEN PRR(W_PR) = $
PRANGE(W_PR,0) > PRR(W_PR) < PRANGE(W_PR,1)
TEST = EXECUTE('F = ' + FNAME + '(X,PRR)')
CRR = FORM_CHISQR( (F - Y)/SIG, ABSOLUTE=ABSOLUTE, $
LORENTZ=LORENTZ)
;
; If an additional improvement was achieved, use the second extrapolation.
;
IF CRR LT CHI(ILO) THEN BEGIN
P(0,IHI) = PRR
CHI(IHI) = CRR
;
; Otherwise use the first extrapolation (reflection).
;
END ELSE BEGIN
P(0,IHI) = PR
CHI(IHI) = CR
ENDELSE
;
; If the reflection yields a value between the highest value and the next
; highest value, use it.
;
END ELSE IF CR GE CHI(INHI) THEN BEGIN
IF CR LT CHI(IHI) THEN BEGIN
P(0,IHI) = PR
CHI(IHI) = CR
ENDIF
;
; Look for an intermediate lower point. If it's an improvement use it.
;
PRR = (P(*,IHI) + PBAR) / 2.
IF PR_COUNT GT 0 THEN PRR(W_PR) = $
PRANGE(W_PR,0) > PRR(W_PR) < PRANGE(W_PR,1)
TEST = EXECUTE('F = ' + FNAME + '(X,PRR)')
CRR = FORM_CHISQR( (F - Y)/SIG, ABSOLUTE=ABSOLUTE, $
LORENTZ=LORENTZ)
IF CRR LT CHI(IHI) THEN BEGIN
P(0,IHI) = PRR
CHI(IHI) = CRR
;
; Nothing seems to help. Contract about the lowest point.
;
END ELSE BEGIN
FOR I = 0,NPAR DO IF I NE ILO THEN BEGIN
PR = (P(*,I) + P(*,ILO)) / 2.
P(0,I) = PR
TEST = EXECUTE('F = ' + FNAME + '(X,PRR)')
CHI(I) = FORM_CHISQR( (F - Y)/SIG, $
ABSOLUTE=ABSOLUTE, LORENTZ=LORENTZ)
ENDIF
ENDELSE
;
; The reflection yields an intermediate value. Use it.
;
END ELSE BEGIN
P(0,IHI) = PR
CHI(IHI) = CR
ENDELSE
;
; Print out the iteration information.
;
DENOM = P(*,IHI)^2 + P(*,ILO)^2
W = WHERE(DENOM NE 0)
ERROR = TOTAL( (P(W,IHI) - P(W,ILO))^2 / DENOM(W) )
BANG_C = !C
CHI_HI = MAX(CHI)
CHI_LO = MIN(CHI)
I_LOW = !C
!C = BANG_C
IF ERROR GT 0 THEN ERRORLG = ALOG10(ERROR)/2 ELSE ERRORLG = -999
IF CHI_HI GT 0 THEN CHILOG = ALOG10(CHI_HI) ELSE CHILOG = -999
IF CHI_LO GT 0 THEN CLOLOG = ALOG10(CHI_LO) ELSE CLOLOG = -999
FORMAT = "(I5,3F15.5,I8)"
IF NOT KEYWORD_SET(NOPRINT) THEN PRINT,FORMAT=FORMAT,ITER,ERRORLG, $
CHILOG,CLOLOG,ITER
IF (ERROR GT ACCUR^2) AND (ITER LT NMAX) THEN GOTO,ITERATE
;
; The program has either converged or reached its maximum number of
; iterations.
;
PARAM = P(*,I_LOW)
CHISQR = CHI(I_LOW)
END