use MERLIN to get motifs

[NimaSarajpoor]

Hi,

In issue #417, the MERLIN algorithm is mentioned and in [WIP] PR #505, as of now, it is being implemented.

Here, I would like to share my idea on using MERLIN to get top-motif (not discord). This is a draft and, for now, I focus on top-motif. It needs to be improved for returning top-k motif.

Main Idea:
Let us assume we have a time series T and we are looking for the top-motif of length m (so we have k = n - m + 1 subsequences). Can we agree that if, for a given min_dist, the algorithm MERLIN returns k-1 candidates, the one that is left is motif? (because it has at least one NN to which its distance is less than min_dist, so already it has priority compared to the ones that are selected as candidates for discords). And, that's it! That is the core idea. Below, I explain how we can implement it.

Implementation:
We can start with min_dist and find some candidates of discords (recall that these cannot be top-motif. So, we just need to scan their right/left neighbors in the first phase very quickly). We store them in a variable that I call it is_cands_cache. It is a Boolean vector that keep track of already-discovered candidates of discords. Then, we can exclude them and search among the remaining ones with lower min_dist. So, we are trying to exclude the candidates of discords!... please see below:

def _find_motif_cands(T, m, min_dist=None, decay=0.1, include=None):
    T, M_T, Σ_T = core.preprocess(T, m)
    
    excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
    
    if min_dist is None:
         min_dist = np.inf

    n = T.shape[0]
    k = n - m + 1
    
    if include is None:
        include = np.ones(k, dtype=bool)
    include_original = include.copy()
    
    approx_P= _get_approx_P(T, m, M_T, Σ_T, s=int(0.001 * len(T)))
    min_dist = min(min_dist, np.max(approx_P[include]))
    
    max_attempts = int(np.ceil(1.0 / decay)) + 1
    min_distances = min_dist * (1 - np.arange(max_attempts) * decay)
    
    min_distances = np.clip(min_distances, a_min=0.0, a_max=None)

    is_cands_cache_before = np.zeros(k, dtype=bool)
    is_cands_cache_updated = np.zeros(k, dtype=bool)
    
    for min_dist in min_distances[:-1]:
        is_cands_init = include & (approx_P >= min_dist) 
        #Note: I am not worried about having an index of an infinite subsequence in `is_cands_init`. 
        #It will be taken care of in _find_candidates(...)
        
        is_cands = _find_candidates(T, m, M_T, Σ_T, min_dist, is_cands_init, direction = 'right')
        is_cands = _find_candidates(T, m, M_T, Σ_T, min_dist, is_cands, direction = 'left')
        
        is_cands_cache_updated[is_cands] = True
        
        if (np.all(is_cands_cache_updated[include_original]: 
            break
        
        is_cands_cache_before[is_cands] = True 
        include[is_cands_cache_before] = False
        
               
    is_cands_cache_before[~include_original]=True
    is_motifs = ~is_cands_cache_before
    
    return is_motifs, min_dist

And,

def find_top_motif(T, m, M_T, Σ_T, is_motifs, approx_P):
    
    motifs_idx = np.flatnonzero(is_motifs)
    
    n = T.shape[0]
    k = n - m + 1 # number of subsequences
    
    excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
    
    #initialize outputs 
    motif_idx = -1
    motif_dist = np.inf
    motif_nn_idx = -1
    
    for idx in motifs_idx:
        if approx_P[idx] > motif_dist:
            continue #because idx cannot be the top-discord among cands_idx
        
        Q = T[idx:idx+m]
        QT = core.sliding_dot_product(Q, T)
        D = core._mass(Q, T, QT, M_T[idx], Σ_T[idx], M_T, Σ_T)
        
        core.apply_exclusion_zone(D, idx, excl_zone, val=np.inf)

        nn_idx = np.argmin(D) 
        
        if D[nn_idx] < motif_dist:
            motif_idx = idx
            motif_dist = D[nn_idx]
            motif_nn_idx = nn_idx
               
    
    return motif_idx, motif_dist, motif_nn_idx

Let us test it on random data:

np.random.seed(seed = 10)
T = np.random.uniform(low=-100.0, high=100.0, size=10000)
m = 500

T, M_T, Σ_T = core.preprocess(T, m)
is_cands_init = np.ones(T.shape[0] - m + 1, dtype=bool)

approx_P = _get_approx_P(T, m, M_T, Σ_T, s = int(0.001*T.shape[0]))

is_motifs, min_dist = _find_motif_cands(T, m, min_dist=np.inf, decay=0.5, include=None)
is_motifs_refined, _ = _find_motif_cands(T, m, min_dist=min_dist, decay=0.25, include=is_motifs.copy())

find_top_motif(T, m, M_T, Σ_T, is_motifs_refined, approx_P)

They give same answer. In stumpy, it takes 0.15 seconds and in the proposed MURLIN-based approach, the running time is 65 second.

---

Please let me know what you think.

[seanlaw]

@ninimama Can you help me understand why we would want to do it this way? What is the benefit compared with stumpy.stump? What are you trying to acheive?

[NimaSarajpoor]

Well, simple yet effective question :) I thought again and apparently, in MURLIN-based approach we are still calculating all pair-wise distances any way and I did not notice it. (I was trying to get rid of some of pair-wise distance calculation and I now think there is no way around it. I mean we HAVE TO find ALL pair-wise distances). And I thought that "I thought enough"!