Window length restrictions and changepoint detection

[xrulesin]

Hi,
Firstly amazing work and this is great piece of code. I have a couple of questions/comments. I am moving from matlab and have been using findchangepts till now. What I noticed is that 1. matlab still allows for single point windows while I get notEnoughpoints exception in python if I try that.

Can we have an option that replicates the closest matlab solution ? Reason I am looking for a single point window is that there might be abrupt changes in the signal (anamoly).. and it is better that it stays isolates as compared to clubbed with 2 points on the left and right.

I am new to this code and I may not be doing things correctly, so please guide me to achieve the desired results.

[deepcharles]

Firstly amazing work and this is great piece of code.

Thanks!

Can you share the code that you run (only the part where you use ruptures)? Also, provide the length of the signal (in number of samples).

[xrulesin]

Thanks Charles,
here is the matlab snippet first (I have 8760 points in my signal).
[fs, fr] = findchangepts(np,'MaxNumChanges',1499,'Statistic','mean','MinDistance',1)
and some of the breakpoints are .... 8203,8218,8219....
where you can see that I get consecutive breakpoint between 8218 and 8219..
Now my python implementation (using Bottomup) and using l2 cost function

import ruptures as rpt
import numpy 
from ruptures.base import BaseCost

signal = numpy.loadtxt("np.txt")
signal = numpy.array(signal)
model='l2'
algo = rpt.BottomUp(model=model,min_size=1, jump=1).fit(signal)
result = algo.predict(n_bkups=1499)

and some of the breakpoints are ... 8203, 8216, 8218, 8224, ...

I notice is that minimum window length for BottomUp is 2 points (8216 and 8218)... irrespective of the cost function that I may try to choose

I am trying to get to a method which can capture the single points in a window (as it may be needed for my project). Can you suggest what changes I can do to achieve that.

[oboulant]

Hi @xrulesin ,
Thanks for you interest in ruptures.

Indeed, in most costs, the min_size is set to 2 because to compute an error on a segment, the error() method for your cost object should be compatible with a segment of size min_size (which is not obvious when you allow for single point segments).

Now, considering your problem, you can create you own Custom Cost and setting the min_size to your desired value. Just be aware that it is then your responsibility to make sure that when using this custom cost, it behaves accordingly in cases where segments of size min_size are considered.

In case of L2 cost, it seems okay to allow to set min_size to 1 since the computed variance on a single point segment is then 0.0.

An exemple here :

import ruptures as rpt
import numpy as np
from ruptures.base import BaseCost
from ruptures.costs import NotEnoughPoints


class MyCost(BaseCost):

    r"""
    Least squared deviation.
    """

    model = "l2"

    def __init__(self):
        """Initialize the object."""
        self.signal = None
        self.min_size = 1

    def fit(self, signal) -> "CostL2":
        """Set parameters of the instance.

        Args:
            signal (array): array of shape (n_samples,) or (n_samples, n_features)

        Returns:
            self
        """
        if signal.ndim == 1:
            self.signal = signal.reshape(-1, 1)
        else:
            self.signal = signal

        return self

    def error(self, start, end) -> float:
        """Return the approximation cost on the segment [start:end].

        Args:
            start (int): start of the segment
            end (int): end of the segment

        Returns:
            segment cost

        Raises:
            NotEnoughPoints: when the segment is too short (less than `min_size` samples).
        """
        if end - start < self.min_size:
            raise NotEnoughPoints

        return self.signal[start:end].var(axis=0).sum() * (end - start)

signal = np.random.normal(0, 0.1, 2000)
algo = rpt.BottomUp(custom_cost=MyCost(),min_size=1, jump=1).fit(signal)
result = algo.predict(n_bkps=1499)
res = [True if result[i] == (result[i-1]+1) else False for i in range(1, len(result)) ]
nb_1sample_segment = np.sum(res)
print(f"There are {nb_1sample_segment} segments of size 1.")

Hope it helps !