Proposal: Implementation of the ODAC (Online Divisive-Agglomerative Clustering) algorithm

[Bezum30]

Introduction

Hi everyone,

This is a proposal for an implementation of the ODAC algorithm for hierarchical clustering of time-series data streams. The implementation focuses on the ODAC algorithm explanation in:

• P. P. Rodrigues, J. Gama and J. Pedroso, "Hierarchical Clustering of Time-Series Data Streams," in IEEE Transactions on Knowledge and Data Engineering, vol. 20, no. 5, pp. 615-627, May 2008, doi: 10.1109/TKDE.2007.190727

To use this algorithm, the odac.py file was put into [path to river]/cluster/. Then, into the __init__.py file (in the same folder), was added:

# import ODAC
from .odac import ODAC

# add "ODAC" in __all__ list
__all__ = [..., "ODAC"]

Algorithm implementation

The content of the odac.py file is the following

from __future__ import annotations

from math import sqrt, log
from typing import Union
from river import base
import numpy as np

class Statistics(base.Base):
    def __init__(self, ts_names:np.ndarray):
        self.sum = {}
        self.prd = {}
        self.n = 0
        for i in range(ts_names.size):
            self.sum[ts_names[i]] = 0.
            for j in range(i,ts_names.size):
                self.prd[(ts_names[i],ts_names[j])] = 0.
      
    def update_sum(self, x:dict) -> None:
        for k in self.sum:
            self.sum[k] += float(x[k])

    def update_product(self, x:dict) -> None:
        for k in self.prd:
            self.prd[k] += float(x[k[0]]) * float(x[k[1]])

    def update_n(self, n=1) -> None:
        self.n += n

"""

Each Node has 1 or more time series, and it is a cluster or not.

"""
class Node(base.Base):

    def __init__(self, name:str, parent: Union['Node', None] = None):
        self.cluster = True
        self.name = name
        self.parent = parent
        self.children = None

        self.d1: float = None
        self.d2: float = None
        
        self.list_of_timeseries_names: np.ndarray = np.array([])

    def __str__(self) -> str:
        return self.__repr__()

    def __repr__(self) -> str:
        pre_0 = "    "
        pre_1 = "│   "
        pre_2 = "├── "
        pre_3 = "└── "
        node = self
        prefix = pre_2 if node.parent is not None and id(node) != id(node.parent.children.second) else pre_3
        while node.parent is not None and node.parent.parent is not None:
            if id(node.parent) != id(node.parent.parent.children.second):
                prefix = pre_1 + prefix
            else:
                prefix = pre_0 + prefix
            node = node.parent
        
        if self.parent is None:
            representation = f"{self.name} d1={self.d1 if self.d1 is not None else '<Not calculated>'}{' d2=' + str(self.d2) if self.d2 is not None else ''}"
        else:
            representation = f"{prefix}{self.name} d1={self.d1 if self.d1 is not None else '<Not calculated>'}{' d2=' + str(self.d2) if self.d2 is not None else ''}"
        if self.cluster is True:
            return representation + f" {self.list_of_timeseries_names}\n"
        else:
            representation += " [NOT A CLUSTER]\n"
        if self.children is not None:
            for child in [self.children.first, self.children.second]:
                representation += child.__repr__()
        return representation

    def __call__(self, ts_names:np.ndarray):
        self.list_of_timeseries_names: np.ndarray = np.sort(ts_names)
        self.statistics = Statistics(ts_names=self.list_of_timeseries_names)

        self.e: float = 0.

        self.d0: float = None
        self.avg: float= None

        self.pivot_0: tuple = None
        self.pivot_1: tuple = None
        self.pivot_2: tuple = None

    def update_statistics(self, x:dict) -> None:
        self.statistics.update_sum(x)
        self.statistics.update_n()
        self.statistics.update_product(x)

    def __calculate_corr_dict(self) -> dict:
        corr_dict = {}

        for i in range(self.list_of_timeseries_names.size):
            for j in range(i+1,self.list_of_timeseries_names.size):
                corr_dict[(self.list_of_timeseries_names[i],self.list_of_timeseries_names[j])] = 0.

        for k in corr_dict:
            p = self.statistics.prd[(k[0], k[1])]
            a = self.statistics.sum[k[0]]
            a2 = self.statistics.prd[(k[0], k[0])]
            b = self.statistics.sum[k[1]]
            b2 = self.statistics.prd[(k[1], k[1])]

            term_p = p - ((a*b)/self.statistics.n)
            term_a = sqrt(a2 - ((a*a)/self.statistics.n))
            term_b = sqrt(b2 - ((b*b)/self.statistics.n))

            corr_dict[(k[0], k[1])] = term_p/(term_a*term_b)
        return corr_dict  
        
    def __calculate_rnomc_dict(self) -> dict:
        rnomc_dict = {}

        for i in range(self.list_of_timeseries_names.size):
            for j in range(i+1,self.list_of_timeseries_names.size):
                rnomc_dict[(self.list_of_timeseries_names[i], self.list_of_timeseries_names[j])] = 0.

        corr_dict = self.__calculate_corr_dict()

        for k in rnomc_dict:
            rnomc_dict[k] = sqrt( (1-corr_dict[k]) / 2 )
            
        return rnomc_dict

    def calculate_coefficients(self, confidence_level:float) -> None:
        rnomc_dict = self.__calculate_rnomc_dict()

        if bool(rnomc_dict):
            self.avg = sum(rnomc_dict.values()) / self.statistics.n

            self.pivot_0 = min(rnomc_dict, key=rnomc_dict.get)
            self.d0 = rnomc_dict[min(rnomc_dict, key=rnomc_dict.get)]

            self.pivot_1 = max(rnomc_dict, key=rnomc_dict.get)
            self.d1 = rnomc_dict[max(rnomc_dict, key=rnomc_dict.get)]

            rnorm_copy = rnomc_dict.copy()

            rnorm_copy.pop(self.pivot_1, None)

            if rnorm_copy is not None and bool(rnorm_copy):
                self.pivot_2 = max(rnorm_copy, key=rnorm_copy.get)
                self.d2 = rnorm_copy[max(rnorm_copy, key=rnorm_copy.get)]
            else:
                self.d2 = None
                self.pivot_2 = None

            r_sqrd = 1
            self.e = sqrt(r_sqrd * log(1/confidence_level) / (2 * self.statistics.n))
        elif None not in [self.avg, self.d0, self.pivot_0, self.d1, self.pivot_1, self.d2, self.pivot_2, self.e]:
            self.avg = None
            self.d0 = None
            self.pivot_0 = None
            self.d1 = None
            self.pivot_1 = None
            self.d2 = None
            self.pivot_2 = None
            self.e = None

    def __get_closer_cluster(self, pivot_1:str, pivot_2:str, current:str, rnormc_dict:dict):
        dist_1 = rnormc_dict[(min(pivot_1, current), max(pivot_1, current))]
        dist_2 = rnormc_dict[(min(pivot_2, current), max(pivot_2, current))]
        return 2 if dist_1 >= dist_2 else 1
    
    def __split_this_cluster(self, pivot_1: str, pivot_2: str, rnormc_dict: dict):
        pivot_set = {pivot_1, pivot_2}
        pivot_1_list:np.ndarray = np.array([pivot_1])
        pivot_2_list:np.ndarray = np.array([pivot_2])

        for ts_name in self.list_of_timeseries_names:
            if ts_name not in pivot_set:
                cluster = self.__get_closer_cluster(pivot_1=pivot_1, pivot_2=pivot_2, current=ts_name, rnormc_dict=rnormc_dict)
                if cluster == 1:
                    pivot_1_list = np.append(pivot_1_list, ts_name)
                else:
                    pivot_2_list = np.append(pivot_2_list, ts_name)

        if self.name == "ROOT":
            new_name = "1"
        else:
            new_name = str(int(self.name.split('_')[-1]) + 1)

        cluster_child_1 = Node("CH1_LVL_" + new_name, parent=self)
        cluster_child_1(pivot_1_list)

        cluster_child_2 = Node("CH2_LVL_" + new_name, parent=self)
        cluster_child_2(pivot_2_list)

        self.children = Children(cluster_child_1, cluster_child_2)

        self.cluster = False

        self.avg = None
        self.d0 = None
        self.pivot_0 = None
        self.pivot_1 = None
        self.pivot_2 = None

        del self.statistics

    def __aggregate_this_cluster(self):
        self.statistics = Statistics(ts_names=self.list_of_timeseries_names)
        self.cluster = True
        if self.children is not None:
            self.children.reset_parent()
            del self.children
            self.children = None

    def test_split(self, tau:float, verbose:bool):
        """
        if d1 - d2 > εk or t > εk then
            if (d1 - d0)|(d1 - d_avg) - (d_avg - d0) > εk then -> split
        """
        
        if self.d2 is not None:
            if ( (self.d1 - self.d2) > self.e ) or ( tau > self.e ) :
                if ( (self.d1 - self.d0) * abs( self.d1 + self.d0 - 2*self.avg) ) > self.e:
                    x1 = self.pivot_1[0]
                    x2 = self.pivot_1[1]
                    if verbose is True:
                        print("#################")
                        print("##### SPLIT #####")
                        print("#################")
                        print(f"self={self.name} parent={self.parent.name if self.parent else 'None'} >>> PIVOTS: x1={x1} x2={x2}")
                        print("#################")
                    self.__split_this_cluster(pivot_1 = x1, pivot_2 = x2, rnormc_dict=self.__calculate_rnomc_dict())
                    return True
        return False

    def test_aggregate(self, verbose:bool):
        """
        The Splitting Criteria guarantees that cluster's diameters monotonically decrease.
            - Assume Clusters: cj with descendants ck and cs;
            - If diameter (ck) - diameter (cj) > ε OR diameter (cs) - diameter (cj ) > ε:
                -> Change in the correlation structure!
                -> Merge clusters ck and cs into cj.
        """

        if self.parent is not None and self.d1 is not None:
            if self.d1 - self.parent.d1 > max(self.parent.e, self.e):
                if verbose is True:
                    print("#################")
                    print("##### AGGR #####")
                    print("#################")
                    print(f"self={self.name} parent={self.parent.name} >>> VARIABLES: c_k={self.d1} c_j={self.parent.d1} e_k={self.e} e_j={self.parent.e}")
                    print("#################")
                self.parent.__aggregate_this_cluster()
                return True
        return False

class Children(base.Base):
    def __init__(self, first:Node, second:Node):
        self.first = first
        self.second = second

    def reset_parent(self) -> None:
        self.first.parent = None
        self.second.parent = None

class ODAC(base.Clusterer):
    def __init__(self, confidence_level = 0.9, n_min = 100, tau = 0.1, verbose = False):
        self.__root_node = Node("ROOT")
        self.confidence_level = confidence_level
        self.n_min = n_min
        self.__observations_threshold = n_min
        self.__n_observations = 0
        self.tau = tau
        self.verbose = verbose
    
    @property
    def n_clusters(self):
        return self.__count_leaves(self.__root_node)
    
    @property
    def n_nodes(self):
        return self.__count_nodes(self.__root_node)
    
    @property
    def height(self) -> int:
        return self.__calculate_height(self.__root_node)
    
    @property
    def summary(self):
        summary = {
            "n_nodes": self.n_nodes,
            "n_clusters": self.n_clusters,
            "height": self.height,
            "supervised": self._supervised,
        }
        return summary

    def __calculate_height(self, node:Node):
        if node.children is not None:
            child_heights = [self.__calculate_height(child) for child in [node.children.first, node.children.second]]
            return max(child_heights) + 1
        else:
            return 0

    def __count_nodes(self, node:Node) -> int:
        count = 1
        if node.children is not None:
            for child in [node.children.first, node.children.second]:
                count += self.__count_nodes(child)
        return count

    def __count_leaves(self, node:Node) -> int:
        if node.cluster is True:
            return 1
        elif node.children is not None:
            return sum(self.__count_leaves(child) for child in [node.children.first, node.children.second])
        else:
            return 0

    def __find_all_leaf_nodes(self, node:Node):
        if node.cluster is True:
            yield node
        elif node.children is not None:
            for child in [node.children.first, node.children.second]:
                yield from self.__find_all_leaf_nodes(child)

    def learn_one(self, new_data: dict):
        """Update the model with a set of features `new_data`.

        Parameters
        ----------
        new_data
            A dictionary of features.

        """
        pass_threshold:bool = False

        if self.__n_observations == 0:
            self.__root_node(np.array(list(new_data.keys())))
            if self.verbose is True:
                print("#################")
                print("##### INIT #####")
                print("#################")
                self.draw()
                print(f"{{'n_nodes': {self.n_nodes}, 'n_clusters': {self.n_clusters}, 'height': {self.height}}}")
                print("#################")
        
        for leaf in self.__find_all_leaf_nodes(self.__root_node):

            if leaf.cluster is False: continue

            leaf.update_statistics(new_data)

            if self.__n_observations >= self.__observations_threshold:
                leaf.calculate_coefficients(confidence_level=self.confidence_level)

                if leaf.test_aggregate(verbose=self.verbose) or leaf.test_split(verbose=self.verbose, tau=self.tau):
                    if self.verbose is True:
                        print(f"ODAC at observation #{self.__n_observations}")
                        self.draw()
                        print(f"{{'n_nodes': {self.n_nodes}, 'n_clusters': {self.n_clusters}, 'height': {self.height}}}")
                        print("#################")
                if pass_threshold is False: 
                    pass_threshold = True
        self.__n_observations += 1
        if pass_threshold is True:
            self.__observations_threshold += self.n_min

    def predict_one(self, x:dict) -> str:
        out:str = None
        if self.__n_observations > 0 and bool(x):
            out = "Not Implemented yet"
        return out

    def draw(self) -> None:
        print(self.__root_node, end="")

Testing

This algorithm was evaluated on Artificial Data, presented in the paper mentioned before (in Section 4.2.1). It was created three files: closed_dataset.csv, two_cluster_dataset.csv, and four_cluster_dataset.csv. You can produce the tests with the following code

from river.cluster import ODAC
from river.stream import iter_csv

verbose = True

model = ODAC(verbose=verbose, confidence_level=0.9, n_min=50)

#################
# optional
params = {
    'converters': { 
        'S0': float,
        'S1': float,
        'S2': float,
        'S3': float,
        'S4': float,
        'S5': float,
        'S6': float,
        'S7': float,
        'S8': float,
        'S9': float,
    },
}
#################

dataset = iter_csv("closed_dataset.csv", **params) # change the file path to test for the three different datasets

for X, _ in dataset:
    model.learn_one(X)
    
model.draw()

print("n_nodes = {}".format(model.n_nodes))

print("n_clusters = {}".format(model.n_clusters))

print("height = {}".format(model.height))

print("summary = {}".format(model.summary))

Please review the proposed implementation and share your thoughts.

Thanks!

[smastelini]

Hi @Bezum30! I received your e-mail and added some review points about your code. Please, send my regards to Professor João Gama!

---

For the record, Gonçalo (@Bezum30) and a colleague are working under the supervision of João Gama to bring the ODAC algorithm to River. A proper PR will be created once we align some details about style and functionality.