aboutsummaryrefslogtreecommitdiff
path: root/pkg/xtools/rmturlach.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/xtools/rmturlach.x')
-rw-r--r--pkg/xtools/rmturlach.x417
1 files changed, 417 insertions, 0 deletions
diff --git a/pkg/xtools/rmturlach.x b/pkg/xtools/rmturlach.x
new file mode 100644
index 00000000..bb23511e
--- /dev/null
+++ b/pkg/xtools/rmturlach.x
@@ -0,0 +1,417 @@
+# Turlach -- Running median library.
+#
+# The algorithm is that described by Haerdle und Steiger (1995) and the
+# implementation is after Turlach. The starting point was the GNU General
+# Pubic Licensed code from the R Foundation (see copyright heritage below).
+# Besides the language recoding the structure has been significantly changed.
+#
+# Copyright (C) 1995 Berwin A. Turlach <berwin@alphasun.anu.edu.au>
+# Copyright (C) 2000-2 Martin Maechler <maechler@stat.math.ethz.ch>
+# Copyright (C) 2003 The R Foundation
+
+include <mach.h>
+
+define RMT_OFFSET 4 # Offset to data
+define RMT_LEN (RMT_OFFSET+5*$1+3) # Structure length
+define RMT_BOX Memi[$1] # Running box size
+define RMT_DATA Memi[$1+1] # Sorted data (ptr)
+define RMT_IN Memi[$1+2] # Mapping to input (ptr)
+define RMT_OUT Memi[$1+3] # Mapping to output (ptr)
+
+define DATA Memr[RMT_DATA($1)+$2]
+define IN Mems[RMT_IN($1)+$2]
+define OUT Mems[RMT_OUT($1)+$2]
+
+
+# RMTURLACH -- Compute running median value using the Turlach algorithm.
+
+real procedure rmturlach (rm, index, data)
+
+pointer rm #I Method pointer
+int index #I Index of new data
+real data #I Input data value
+
+short nrnew, box, outnext, out, leaf, one
+data one/1/
+
+begin
+ nrnew = index - 1
+ box = RMT_BOX(rm)
+ outnext = mod (nrnew, box)
+ out = OUT(rm,outnext)
+ DATA(rm,out) = data
+
+ leaf = out - box
+ if (out > box) {
+ if (data >= DATA(rm,box))
+ call rm_uoui (leaf, box, DATA(rm,1), OUT(rm,1), IN(rm,1))
+ else
+ call rm_uodi (leaf, box, nrnew, outnext, data,
+ DATA(rm,1), OUT(rm,1), IN(rm,1))
+ } else if (out < box) {
+ if (data < DATA(rm,box))
+ call rm_dodi (leaf, box, DATA(rm,1), OUT(rm,1), IN(rm,1))
+ else
+ call rm_doui (leaf, box, nrnew, outnext, data,
+ DATA(rm,1), OUT(rm,1), IN(rm,1))
+ } else if (DATA(rm,box) > DATA(rm,box+1)) {
+ call rm_swap (box, box+one, DATA(rm,1), OUT(rm,1), IN(rm,1));
+ call rm_uptoleaf (one, box, DATA(rm,1), OUT(rm,1), IN(rm,1));
+ } else if (DATA(rm,box) < DATA(rm,box-1)) {
+ call rm_swap (box, box-one, DATA(rm,1), OUT(rm,1), IN(rm,1));
+ call rm_downtoleaf (-one, box, DATA(rm,1), OUT(rm,1), IN(rm,1));
+ }
+
+ return (DATA(rm,box))
+end
+
+
+# RMT_OPEN -- Open Turlach running median algorithm.
+
+pointer procedure rmt_open (box, data)
+
+int box #I Running median box
+real data #I Initial data value
+pointer rm #R Method pointer
+
+short i, halfbox
+#short i, j, k, halfbox, one
+#data one/1/
+
+begin
+ call malloc (rm, RMT_LEN(box), TY_STRUCT)
+
+ RMT_BOX(rm) = box
+ RMT_DATA(rm) = rm + RMT_OFFSET
+ RMT_IN(rm) = P2S(RMT_DATA(rm) + 2 * box + 1)
+ RMT_OUT(rm) = RMT_IN(rm) + 2 * box + 1
+
+ halfbox = (box - 1) / 2
+
+ do i = 1+halfbox, box+halfbox {
+ DATA(rm,i) = data
+ IN(rm,i) = i-halfbox-1
+ OUT(rm,i-halfbox-1) = i
+ }
+
+ do i = 0, halfbox {
+ DATA(rm,i) = -MAX_REAL
+ DATA(rm,i+box+halfbox+1) = MAX_REAL
+ }
+
+ return (rm)
+end
+
+
+# RMT_CLOSE -- Close Turlach running median algorithm.
+
+procedure rmt_close (rm)
+
+pointer rm #I Method pointer
+
+begin
+ call mfree (rm, TY_STRUCT)
+end
+
+
+# RMT_DUMP -- Dump data structure.
+
+procedure rmt_dump (rm, unsorted, sorted, in, out)
+
+pointer rm #I Method pointer
+bool unsorted #I Dump data in unsorted order?
+bool sorted #I Dump data in sorted order?
+bool in #I Dump in list?
+bool out #I Dump out list?
+
+int i, box, halfbox
+
+begin
+ box = RMT_BOX(rm)
+ halfbox = box / 2
+ if (unsorted) {
+ do i = 1+halfbox, halfbox+box {
+ call eprintf (" %3.0f")
+ call pargr (DATA(rm,OUT(rm,i-halfbox-1)))
+ }
+ call eprintf ("\n")
+ }
+ if (sorted) {
+ #do i = 0, 2*box {
+ do i = 1+halfbox, halfbox+box {
+ call eprintf (" %3.0f")
+ call pargr (DATA(rm,i))
+ }
+ call eprintf ("\n")
+ }
+ if (in) {
+ #do i = 0, 2*box {
+ do i = 1+halfbox, halfbox+box {
+ call eprintf (" %3d")
+ call pargs (IN(rm,i))
+ }
+ call eprintf ("\n")
+ }
+ if (out) {
+ do i = 0, box-1 {
+ call eprintf (" %3d")
+ call pargs (OUT(rm,i))
+ }
+ call eprintf ("\n")
+ }
+end
+
+
+# RM_SWAP -- Swap positions `l' and `r'.
+
+procedure rm_swap (l, r, window, outlist, nrlist)
+
+short l #I Index to swap
+short r #I Index to swap
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short nl, nr
+real w
+
+begin
+ w = window[l]; window[l] = window[r]; window[r] = w
+ nl = nrlist[l]; nr = nrlist[r]; nrlist[l] = nr; nrlist[r] = nl
+ outlist[nl] = r; outlist[nr] = l
+end
+
+
+# RM_SIFTUP -- Used only in the initial sorting.
+
+procedure rm_siftup (l, r, window, outlist, nrlist)
+
+short l #I Left index
+short r #I Right index
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short i, j, nrold
+real w
+
+begin
+ i = l
+ j = 2 * i
+ w = window[i]
+ nrold = nrlist[i]
+ while (j <= r) {
+ if (j < r) {
+ if (window[j] < window[j+1])
+ j = j + 1
+ }
+ if (w >= window[j])
+ break
+
+ window[i] = window[j]
+ outlist[nrlist[j]] = i
+ nrlist[i] = nrlist[j]
+ i = j
+ j = 2 * i
+ }
+
+ window[i] = w
+ outlist[nrold] = i
+ nrlist[i] = nrold
+end
+
+
+# RM_UOUI - Upper Out Upper In
+
+procedure rm_uoui (leaf, box, window, outlist, nrlist)
+
+short leaf #I Leaf
+short box #I Box size
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short i, j, k
+
+begin
+ call rm_uptoleaf (leaf, box, window, outlist, nrlist)
+
+ i = leaf
+ j = i + box
+ k = i / 2 + box
+ while (window[j] < window[k]) {
+ call rm_swap (j, k, window, outlist, nrlist)
+ i = (k - box)
+ j = i + box
+ k = i / 2 + box
+ }
+end
+
+
+# RM_DODI - Down Out Down In
+
+procedure rm_dodi (leaf, box, window, outlist, nrlist)
+
+short leaf #I Leaf
+short box #I Box size
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short i, j, k
+
+begin
+ call rm_downtoleaf (leaf, box, window, outlist, nrlist)
+
+ i = leaf
+ j = i + box
+ k = i / 2 + box
+ while (window[j] > window[k]) {
+ call rm_swap (j, k, window, outlist, nrlist)
+ i = (k - box)
+ j = i + box
+ k = i / 2 + box
+ }
+end
+
+
+# RM_UODI -- Upper Out Down In
+
+procedure rm_uodi (leaf, box, nrnew, outnext, in, window, outlist, nrlist)
+
+short leaf #I Leaf
+short box #I Box size
+short nrnew #I nrnew
+short outnext #I outnext
+real in #I Input value
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short one
+data one/1/
+
+begin
+ call rm_toroot (leaf, box, nrnew, outnext, in, window, outlist,
+ nrlist)
+
+ if (window[box] < window[box-1]) {
+ call rm_swap (box, box-one, window, outlist, nrlist)
+ call rm_downtoleaf (-one, box, window, outlist, nrlist)
+ }
+end
+
+
+# RM_DOUI -- Down Out Upper In
+
+procedure rm_doui (leaf, box, nrnew, outnext, in, window, outlist, nrlist)
+
+short leaf #I Leaf
+short box #I Box size
+short nrnew #I nrnew
+short outnext #I outnext
+real in #I Input value
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short one
+data one/1/
+
+begin
+ call rm_toroot (leaf, box, nrnew, outnext, in, window, outlist,
+ nrlist)
+
+ if (window[box] > window[box+1]) {
+ call rm_swap (box, box+one, window, outlist, nrlist)
+ call rm_uptoleaf (one, box, window, outlist, nrlist)
+ }
+end
+
+# RM_TOROOT
+
+procedure rm_toroot (leaf, box, nrnew, outnext, in, window, outlist, nrlist)
+
+short leaf #I Leaf
+short box #I Box size
+short nrnew #I nrnew
+short outnext #I outnext
+real in #I Input value
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short i, j, k
+
+begin
+ i = leaf
+ repeat {
+ j = i + box
+ k = i / 2 + box
+ window[j] = window[k]
+ outlist[nrlist[k]] = j
+ nrlist[j] = nrlist[k]
+ i = k - box
+ } until (i == 0)
+
+ window[box] = in
+ outlist[outnext] = box
+ nrlist[box] = outnext
+end
+
+
+# RM_DOWNTOLEAF
+
+procedure rm_downtoleaf (leaf, box, window, outlist, nrlist)
+
+short leaf #I Leaf
+short box #I Box size
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short i, j, childl, childr
+
+begin
+ i = leaf
+ repeat {
+ j = i + box
+ childl = 2 * i + box
+ childr = childl - 1
+ if (window[childl] < window[childr])
+ childl = childr
+ if (window[j] >= window[childl])
+ break
+ call rm_swap (j, childl, window, outlist, nrlist)
+ i = childl - box
+ }
+end
+
+
+# RM_UPTOLEAF
+
+procedure rm_uptoleaf (leaf, box, window, outlist, nrlist)
+
+short leaf #I Leaf
+short box #I Box size
+real window[ARB] #U Work array
+short outlist[ARB] #U Work array
+short nrlist[ARB] #U Work array
+
+short i, j, childl, childr
+
+begin
+ i = leaf
+ repeat {
+ j = i + box
+ childl = 2 * i + box
+ childr = childl + 1
+ if (window[childl] > window[childr])
+ childl = childr
+ if (window[j] <= window[childl])
+ break
+ call rm_swap (j, childl, window, outlist, nrlist)
+ i = childl - box
+ }
+end
+
generated by cgit (git 2.34.1) at 2025-09-20 08:41:42 -0400