(FILECREATED "23-Feb-86 17:27:28" {DSK}<LISPFILES>DINDESMOOTH.;7 17874
changes to: (FNS SPLIT RunningSmooth PairSmooth S3RSSH S4 S2 S5 S3 HANN RUNMED MEDOF4 ENDPTS)
(VARS DINDESMOOTHCOMS)
(MACROS MEDOF3)
previous date: " 3-Feb-86 18:52:44" {DSK}<LISPFILES>DINDESMOOTH.;5)
(* Copyright (c) 1986 by Massachusetts Institute of Technology. All rights reserved.)
(PRETTYCOMPRINT DINDESMOOTHCOMS)
(RPAQQ DINDESMOOTHCOMS ((FNS RSM RunningSmooth PairSmooth S3RSSH S4253H S4 S2 S5 S3 S3R HANN ENDPTS
SPLIT RUNMED MEDOF3)
(MACROS MEDOF3)))
(DEFINEQ
(RSM
[LAMBDA (Y VERSN SMOOTH ROUGH) (* SCP "16-Feb-86 18:26")
(* * NONLINEAR SMOOTHERS. REFERENCE: HOAGLIN AND VELLEMAN, "THE ABC'S OF EDA")
(* * ON ENTRY: -
THE RANK 1 FLONUM CMLARRAY Y IS A DATA SEQUENCE OF N VALUES. VERSN SPECIFIES THE SMOOTHER TO BE USED: EITHER
3RSSH,TWICE OR 4253H,TWICE)
(* * ON EXIT: -
RETURNS THE CONS OF SMOOTH AND ROUGH WHERE CMLARRAYS SMOOTH AND ROUGH CONTAIN THE RESULT FROM THE SMOOTHING
OPERATION. Y (I) = SMOOTH (I) + ROUGH (I) FOR EACH I FROM 0 TO N-1.)
(if (NOT (EQP (ARRAYRANK Y)
1))
then (ERROR "Array not 1-Dimensional:" Y))
(if (NOT (EQP (ARRAYELEMENTTYPE Y)
(QUOTE FLONUM)))
then (ERROR "Array not FLONUM:" Y))
(if (LEQ (ARRAYTOTALSIZE Y)
6)
then (ERROR "The CMLARRAY passed to RSM has fewer than 7 components:" Y))
(if (NULL SMOOTH)
then (SETQ SMOOTH (MAKEARRAY (ARRAYTOTALSIZE Y)
(QUOTE ELEMENTTYPE)
(QUOTE FLONUM)))
else (if (NOT (EQP (ARRAYDIMENSIONS SMOOTH)
(ARRAYDIMENSIONS Y)))
then (ERROR "Array supplied for SMOOTH is not the same shape as data."))
(if (NOT (EQP (ARRAYELEMENTTYPE SMOOTH)
(QUOTE FLONUM)))
then (ERROR "Array not FLONUM:" Y)))
(if (NULL ROUGH)
then (SETQ ROUGH (MAKEARRAY (ARRAYTOTALSIZE Y)
(QUOTE ELEMENTTYPE)
(QUOTE FLONUM)))
else (if (NOT (EQP (ARRAYDIMENSIONS ROUGH)
(ARRAYDIMENSIONS Y)))
then (ERROR "Array supplied for ROUGH is not the same shape as data."))
(if (NOT (EQP (ARRAYELEMENTTYPE ROUGH)
(QUOTE FLONUM)))
then (ERROR "Array not FLONUM:" Y)))
(for I from 0 to (SUB1 (ARRAYTOTALSIZE Y)) do (\FLOATASET (\FLOATAREF Y I)
SMOOTH I))
(SELECTQ VERSN
(3RSSH,TWICE (S3RSSH SMOOTH))
(4253H,TWICE (S4253H SMOOTH))
(ERROR "The choice of smooth is invalid:" VERSN))
(* * COMPUTE ROUGH FROM FIRST SMOOTHING)
(for I from 0 to (SUB1 (ARRAYTOTALSIZE Y)) do (\FLOATASET (DIFFERENCE (\FLOATAREF Y I)
(\FLOATAREF SMOOTH I))
ROUGH I))
(* * REROUGH SMOOTHERS ("TWICING"))
(SELECTQ VERSN
(3RSSH,TWICE (S3RSSH ROUGH))
(4253H,TWICE (S4253H ROUGH))
(ERROR "The choice of smooth is invalid:" VERSN))
(for I from 0 to (SUB1 (ARRAYTOTALSIZE Y))
do (\FLOATASET (PLUS (\FLOATAREF SMOOTH I)
(\FLOATAREF ROUGH I))
SMOOTH I)
(\FLOATASET (DIFFERENCE (\FLOATAREF Y I)
(\FLOATAREF SMOOTH I))
ROUGH I))
(CONS SMOOTH ROUGH])
(RunningSmooth
[LAMBDA (Y Smoother Twicing?) (* SCP "23-Feb-86 17:08")
(* * Nonlinear smoothers. Reference: Hoaglin and Velleman, "the ABCs of EDA")
(* * Y is a rank 1 FLONUM CMLArray to be smoothed.)
(* * Smoother specifies the smoother to be used: currently one of S3, S3R, S4, S5, HANN, and two compound smoothers
S3RSSH and S4253H. Default is S3R.)
(* * Twicing? if present indicates that the smooth should be "reroughed" in the sense that the rough is also
smoothed and the result added back into the original smooth.)
(* * Returns a new CMLArray containing the smoothed values of Y.)
(LET ((N (ARRAYTOTALSIZE Y))
(SmootherName (if (NULL Smoother)
then (QUOTE S3R)
else Smoother))
SmoothedY)
(* * Error checking on the incoming arrays.)
(if (NOT (EQP (ARRAYRANK Y)
1))
then (ERROR "Array not 1-Dimensional:" Y))
(if (NOT (EQP (ARRAYELEMENTTYPE Y)
(QUOTE FLONUM)))
then (ERROR "Array not FLONUM:" Y))
(if (LEQ N 6)
then (ERROR "The CMLARRAY passed to RunningSmooth must have at least 7 components:" Y))
(SETQ SmoothedY (MAKEARRAY N (QUOTE ELEMENTTYPE)
(ARRAYELEMENTTYPE Y)
(QUOTE INITIALCONTENTS)
Y))
(for I from 0 to (SUB1 N) do (\FLOATASET (\FLOATAREF Y I)
SmoothedY I))
(* * Now compute the appropriate smooth.)
(SELECTQ SmootherName
(S3R (S3R SmoothedY))
(S3RSSH (S3RSSH SmoothedY))
(S4253H (S4253H SmoothedY))
(S3 (S3 SmoothedY))
(S4 (S4 SmoothedY))
(S5 (S5 SmoothedY))
(HANN (HANN SmoothedY))
(ERROR "The choice of smooth is invalid:" SmootherName))
(* * Rerough if twicing is requested.)
[if Twicing?
then (LET [(RoughOfY (MAKEARRAY N (QUOTE ELEMENTTYPE)
(QUOTE FLONUM]
(for I from 0 to (SUB1 N) do (\FLOATASET (DIFFERENCE (\FLOATAREF Y I)
(\FLOATAREF SmoothedY I))
RoughOfY I))
(* * Having computed the RoughOfY from the smooth.)
(SELECTQ SmootherName
(S3R (S3R RoughOfY))
(S3RSSH (S3RSSH RoughOfY))
(S4253H (S4253H RoughOfY))
(S3 (S3 RoughOfY))
(S4 (S4 RoughOfY))
(S5 (S5 RoughOfY))
(HANN (HANN RoughOfY))
(ERROR "The choice of smooth is invalid:" SmootherName))
(for I from 0 to (SUB1 N) do (\FLOATASET (PLUS (\FLOATAREF SmoothedY I)
(\FLOATAREF RoughOfY I))
SmoothedY I]
(* * Finally return the smooth)
SmoothedY])
(PairSmooth
[LAMBDA (Y X Smoother Twicing?) (* SCP "23-Feb-86 17:20")
(* Returns a CMLArray, the smooth version of the Y smoothed along the sorted values of X, using the smoother given
by the argument Smoother. Default Smoother is 3RSSH,TWICE. This really makes sense only when X is
(nearly) equally spaced, otherwise a scatterplot smoother should be used!)
(LET ((N (ARRAYTOTALSIZE Y)))
(if (NULL X)
then (RunningSmooth Y Smoother Twicing?)
elseif (NOT (EQP N (ARRAYTOTALSIZE X)))
then (ERROR "CMLArrays X and Y don't have the same dimensions.")
else (LET [(WorkY (MAKEARRAY N (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)))
(WorkX (MAKEARRAY N (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)))
(Index (MAKEARRAY N (QUOTE ELEMENTTYPE)
(QUOTE FLONUM] (* Set up some work space since lower sorting functions
are by side-effect i.e. here they change their
arguments.)
(for I from 0 to (SUB1 N)
do (\FLOATASET (\FLOATAREF X I)
WorkX I)
(\FLOATASET (\FLOATAREF Y I)
WorkY I)
(\FLOATASET (FLOAT I)
Index I))
(UFQuickPairSort WorkX Index) (* Store original order of the X's since PairSort is
destructive to WorkX and Index.)
(for I from 0 to (SUB1 N) do (\FLOATASET (\FLOATAREF X I)
WorkX I))
(* Get the X's back again.)
(UFQuickPairSort WorkX WorkY) (* Now sort the Y's according to X)
(SETQ WorkY (Smooth WorkY Smoother Twicing?))
(* Finally Smooth Y and then sort the result so that it
again matches X.)
(UFQuickPairSort Index WorkY)
WorkY])
(S3RSSH
[LAMBDA (Y) (* SCP "23-Feb-86 16:55")
(* * SMOOTH CMLARRAY Y BY 3RSSH.)
(S3R Y)
(if (SPLIT Y)
then (S3R Y)
(if (SPLIT Y)
then (S3R Y)))
(HANN Y])
(S4253H
[LAMBDA (Y) (* SCP " 3-Dec-85 15:39")
(* * SMOOTH BY 4253H)
(LET ((ENDSAV (S4 Y)))
(S2 Y ENDSAV)
(S5 Y)
(S3 Y)
(ENDPTS Y)
(HANN Y])
(S4
[LAMBDA (Y) (* SCP "23-Feb-86 16:54")
(* * SMOOTH BY RUNNING MEDIANS OF 4)
(* * EVEN LENGTH MEDIAN OFFSET THE OUTPUT SEQUENCE TO THE HIGH END, -
SINCE THEY CANNOT BE SYMMETRIC. ENDSAV IS LEFT HOLDING Y (N-1) SINCE -
THERE IS NO OTHER ROOM FOR IT. Y (0) IS UNCHANGED.)
(LET* [(N (ARRAYTOTALSIZE Y))
(ENDSAV (AREF Y (SUB1 N)))
(ENDM1 (AREF Y (DIFFERENCE N 2]
(RUNMED Y 4)
(ASET (TIMES (PLUS (AREF Y 0)
(AREF Y 1))
.5)
Y 1)
(ASET (TIMES (PLUS ENDM1 ENDSAV)
.5)
Y
(SUB1 N))
ENDSAV])
(S2
[LAMBDA (Y ENDSAV) (* SCP "23-Feb-86 16:54")
(* * SMOOTH BY RUNNING MEDIANS (MEANS) OF 2 -
USED TO RECENTER RESULTS OF RUNNING MEDIANS OF 4 -
ENDSAV HOLDS THE ORIGINAL Y{N-1}.)
(LET ((N (ARRAYTOTALSIZE Y)))
(for I from 1 to (DIFFERENCE N 2) do (\FLOATASET (TIMES (PLUS (\FLOATAREF Y (ADD1 I))
(\FLOATAREF Y I))
.5)
Y I))
(\FLOATASET ENDSAV Y (SUB1 N])
(S5
[LAMBDA (Y) (* SCP "23-Feb-86 16:53")
(* * SMOOTH BY RUNNING MEDIANS OF 5)
(LET ((N (ARRAYTOTALSIZE Y))
YMED1 YMED2)
(SETQ YMED1 (MEDOF3 (AREF Y 0)
(AREF Y 1)
(AREF Y 2)))
[SETQ YMED2 (MEDOF3 (AREF Y (SUB1 N))
(AREF Y (DIFFERENCE N 2))
(AREF Y (DIFFERENCE N 3]
(RUNMED Y 5)
(ASET YMED1 Y 1)
(ASET YMED2 Y (DIFFERENCE N 2])
(S3
[LAMBDA (Y) (* SCP "23-Feb-86 16:52")
(* * COMPUTE RUNNING MEDIAN OF 3 ON CMLARRAY Y. -
RETURNS T IF ANY CHANGE IS MADE.)
(bind ANS Y1 (CHANGED ← NIL)
(Y2 ←(\FLOATAREF Y 0))
(Y3 ←(\FLOATAREF Y 1)) declare (TYPE FLOATP ANS Y1 Y2 Y3) for I from 1
to (DIFFERENCE (ARRAYTOTALSIZE Y)
2)
do (SETQ Y1 Y2)
(SETQ Y2 Y3)
(SETQ Y3 (\FLOATAREF Y (ADD1 I)))
(SETQ ANS (MEDOF3 Y1 Y2 Y3))
(if (NOT (UFEQP ANS Y2))
then (SETQ CHANGED T)
(\FLOATASET ANS Y I))
finally (RETURN CHANGED])
(S3R
[LAMBDA (Y) (* RWO "19-Nov-85 15:31")
(* * COMPUTE REPEATED RUNNING MEDIANS OF 3.0)
(repeatwhile (S3 Y))
(ENDPTS Y])
(HANN
[LAMBDA (Y) (* SCP "23-Feb-86 16:52")
(* * 3-POINT SMOOTH BY MOVING AVERAGES WEIGHTED 1/4, 1/2, 1/4. -
THIS IS CALLED HANNING.)
(bind Y1 (Y2 ←(\FLOATAREF Y 0))
(Y3 ←(\FLOATAREF Y 1)) declare (TYPE FLOATP Y1 Y2 Y3) for I from 1
to (DIFFERENCE (ARRAYTOTALSIZE Y)
2)
do (SETQ Y1 Y2)
(SETQ Y2 Y3)
(SETQ Y3 (\FLOATAREF Y (ADD1 I)))
(\FLOATASET (TIMES (PLUS Y1 Y2 Y2 Y3)
.25)
Y I])
(ENDPTS
[LAMBDA (Y) (* SCP "16-Feb-86 17:42")
(* * ESTIMATE SMOOTHED VALUES FOR BOTH END POINTS OF THE SEQUENCE IN Y USING THE END POINT EXTRAPOLATION RULE.
ALL THE OTHER VALUES IN Y HAVE ALREADY BEEN SMOOTHED.)
(LET ((N (ARRAYTOTALSIZE Y))
Y0)
(DECLARE (TYPE FLOATP Y0))
(* * LEFT END)
[SETQ Y0 (DIFFERENCE (TIMES 3.0 (\FLOATAREF Y 1))
(TIMES 2.0 (\FLOATAREF Y 2]
(\FLOATASET (MEDOF3 Y0 (\FLOATAREF Y 0)
(\FLOATAREF Y 1))
Y 0)
(* * RIGHT END)
[SETQ Y0 (DIFFERENCE (TIMES 3.0 (\FLOATAREF Y (IDIFFERENCE N 2)))
(TIMES 2.0 (\FLOATAREF Y (IDIFFERENCE N 3]
(\FLOATASET (MEDOF3 Y0 (\FLOATAREF Y (SUB1 N))
(\FLOATAREF Y (IDIFFERENCE N 2)))
Y
(SUB1 N])
(SPLIT
[LAMBDA (Y) (* SCP "23-Feb-86 17:26")
(* * FIND 2-FLATS IN Y AND APPLY SPLITTING ALGORITHM.)
(* W IS A WINDOW 6 WIDE WHICH SLIDES ALONG Y.)
(bind W0 (W1 ←(\FLOATAREF Y 2))
(W2 ←(\FLOATAREF Y 0))
(W3 ←(\FLOATAREF Y 1))
(W4 ←(\FLOATAREF Y 2))
(W5 ←(\FLOATAREF Y 3))
(N ←(DIFFERENCE (ARRAYTOTALSIZE Y)
1))
(NM2 ←(DIFFERENCE (ARRAYTOTALSIZE Y)
3))
CHANGE Y1 ANS declare (TYPE FLOATP W0 W1 W2 W3 W4 W5 Y1 ANS) for I from 0
to (SUB1 (DIFFERENCE (ARRAYTOTALSIZE Y)
1))
do
(* * IF Y (0) =Y (1) .NE. Y (2), TREAT FIRST 2 LIKE A 2-FLAT WITH END PT RULE)
[if (UFEQP W2 W3)
then (if (OR (AND (UFLESSP W2 W1)
(UFLESSP W3 W4))
(AND (UFGREATERP W2 W1)
(UFGREATERP W3 W4)))
then
(* * W (2) AND W (3) FORM A 2-FLAT)
[if (GEQ I 2)
then
(* * APPLY RIGHT END PT RULE AT I)
(SETQ Y1 (DIFFERENCE (TIMES 3.0 W1)
(TIMES 2.0 W0)))
(SETQ ANS (MEDOF3 Y1 W2 W1))
(\FLOATASET ANS Y I)
(SETQ CHANGE (NOT (UFEQP ANS W2]
(if (LESSP I NM2)
then
(* * APPLY LEFT END PT RULE AT I+1)
(SETQ Y1 (DIFFERENCE (TIMES 3.0 W4)
(TIMES 2.0 W5)))
(SETQ ANS (MEDOF3 Y1 W3 W4))
(\FLOATASET ANS Y (ADD1 I))
(SETQ CHANGE (NOT (UFEQP ANS W3]
(* * SLIDE WINDOW)
(SETQ W0 W1)
(SETQ W1 W2)
(SETQ W2 W3)
(SETQ W3 W4)
(SETQ W4 W5)
(if (LESSP I (DIFFERENCE N 3))
then (SETQ W5 (\FLOATAREF Y (PLUS I 4)))
else
(* * APPLY RULE TO LAST 2 POINTS IF NEEDED.)
(SETQ W5 W2))
finally (RETURN CHANGE])
(RUNMED
[LAMBDA (Y LEN) (* SCP "21-Feb-86 21:33")
(* * SMOOTH CMLARRAY Y BY RUNNING MEDIANS OF LENGTH LEN. NOTE: USE S3 FOR RUNNING MEDIANS OF 3 INSTEAD OF RUNMED.)
(* SAVE ACTS LIKE A WINDOW OF LENGTH LEN WHICH SLIDES DOWN Y. IT WORKS CIRCULARLY, THE OLDEST ELEMENT IS REPLACED BY
THE NEW ARRIVAL.)
(* THE COMPONENTS OF WORK ARE THE SAME AS THOSE OF SAVE
POSSIBLY IN A DIFFERENT ORDER AS A SIDE EFFECT OF
UFSelectMedian.)
(LET [(SAVE (MAKEARRAY LEN (QUOTE ELEMENTTYPE)
(QUOTE FLONUM)))
(WORK (MAKEARRAY LEN (QUOTE ELEMENTTYPE)
(QUOTE FLONUM]
(for I from 0 to (SUB1 LEN)
do (\FLOATASET (\FLOATAREF Y I)
SAVE I)
(\FLOATASET (\FLOATAREF Y I)
WORK I))
(bind (SAVEPT ← 0)
TEMP declare (TYPE FLOATP TEMP) for I from LEN to (SUB1 (ARRAYTOTALSIZE Y)) as SMOPT
from (QUOTIENT LEN 2)
do (\FLOATASET (SETQ TEMP (UFSelectMedian WORK))
Y SMOPT)
[bind (TEMP ←(\FLOATAREF SAVE SAVEPT)) declare (TYPE FLOATP TEMP) for J from 0
to (SUB1 LEN) until (UFEQP (\FLOATAREF WORK J)
TEMP)
finally ((\FLOATASET (\FLOATAREF Y I)
WORK J)
(\FLOATASET (\FLOATAREF Y I)
SAVE SAVEPT)
(SETQ SAVEPT (IMOD (ADD1 SAVEPT)
LEN]
finally (\FLOATASET (SETQ TEMP (UFSelectMedian WORK))
Y SMOPT])
(MEDOF3
[LAMBDA (X1 X2 X3) (* SCP "16-Feb-86 17:35")
(* * RETURN THE MEDIAN OF X1, X2, X3)
(if (GREATERP X2 X1)
then (if (GREATERP X3 X2)
then X2
elseif (GREATERP X3 X1)
then X3
else X1)
elseif (LESSP X3 X2)
then X2
elseif (LESSP X3 X1)
then X3
else X1])
)
(DECLARE: EVAL@COMPILE
(PUTPROPS MEDOF3 DMACRO ((X1 X2 X3)
(if (UFGREATERP X2 X1)
then
(if (UFGREATERP X3 X2)
then X2 elseif (UFGREATERP X3 X1)
then X3 else X1)
elseif
(UFLESSP X3 X2)
then X2 elseif (UFLESSP X3 X1)
then X3 else X1)))
)
(PUTPROPS DINDESMOOTH COPYRIGHT ("Massachusetts Institute of Technology" 1986))
(DECLARE: DONTCOPY
(FILEMAP (NIL (581 17468 (RSM 591 . 3564) (RunningSmooth 3566 . 6538) (PairSmooth 6540 . 8559) (S3RSSH
8561 . 8850) (S4253H 8852 . 9123) (S4 9125 . 9843) (S2 9845 . 10375) (S5 10377 . 10910) (S3 10912 .
11612) (S3R 11614 . 11829) (HANN 11831 . 12403) (ENDPTS 12405 . 13308) (SPLIT 13310 . 15400) (RUNMED
15402 . 17028) (MEDOF3 17030 . 17466)))))
STOP