rrcf-demo-run.py

# Copyright (c) 2022 Kohei "Max" MATSUSHITA (ma2shita+git@ma2shita.jp)
# Released under the MIT license
# https://opensource.org/licenses/mit-license.php

import sys
# Params for sin()
hz = 1.0
length_sec = 4.1
sampling_frequency = 10
# Generate
import numpy as np
np.set_printoptions(precision=2, floatmode='fixed', suppress=True)
data = np.sin(2 * np.pi * hz * np.arange(0, length_sec, 1/sampling_frequency))

data[31:33] = -1 # inject anomaly

# Hyper-params for RRCF
shingle_size = 4
num_trees = 100
tree_size = 60 # おまじない

# Compute CoDisp by RRCF
import rrcf
forest = []
for _ in range(num_trees):
    tree = rrcf.RCTree()
    forest.append(tree)
points = rrcf.shingle(data, size=shingle_size)
avg_codisp = {}
for index, point in enumerate(points):
    # For each tree in the forest...
    for tree in forest:
        # If tree is above permitted size...
        if len(tree.leaves) > tree_size:
            # Drop the oldest point (FIFO)
            tree.forget_point(index - tree_size)
        # Insert the new point into the tree
        tree.insert_point(point, index=index)
        # Compute codisp on the new point...
        new_codisp = tree.codisp(index)
        # And take the average over all trees
        if not index in avg_codisp:
            avg_codisp[index] = 0
        avg_codisp[index] += new_codisp / num_trees
print(avg_codisp, file=sys.stderr) # dump

# Rendering plot
import matplotlib.pyplot as plt
_, ax1 = plt.subplots()
ax1.plot(data, color='tab:cyan', marker="")
ax2 = ax1.twinx()
ax2.set_ylim(0, 20) # suppressing axis top
ax2.plot(avg_codisp.values(), color='tab:red')
plt.savefig(sys.stdout.buffer)