# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import metric, motif
plt.rcParams['ps.useafm'] = True
plt.rcParams['pdf.use14corefonts'] = True
plt.rcParams['text.usetex'] = True

# +
# Parameters for synthesized time series
quarter_size = 8
pattern_size = quarter_size * 4
num_samples = 100
std0 = 1
std1 = 0.05
std2 = 2 * std1

# Parameters for motif enumeration methods
window_size = pattern_size
radiuses = np.linspace(0.01, 2.0, 200)

# +
# Define motif-1 (triangular wave)
pattern1 = np.concatenate([
    np.arange(0, quarter_size),
    np.arange(quarter_size, -quarter_size, -1),
    np.arange(-quarter_size, 0)
]).astype('float')

# Define motif-2 (sine wave)
pattern2 = np.sin(np.linspace(0, 2 * np.pi, pattern_size))

# Standardize motifs
pattern1 = metric.standardize(pattern1, axis=0)
pattern2 = metric.standardize(pattern2, axis=0)

# +
# Initialize random seed
np.random.seed(0)

# Generate patterns
X1 = np.random.normal(0, std1, [num_samples, pattern_size])
X1 += pattern1
X2 = np.random.normal(0, std2, [num_samples, pattern_size])
X2 += pattern2

# Standardize patterns
X1std = metric.standardize(X1.copy().T, axis=0).T
X2std = metric.standardize(X2.copy().T, axis=0).T

# Synthesize time seirs
period = 3 * pattern_size
X = np.random.normal(0, std0, [num_samples, period])
begin1 = 2 * quarter_size
end1 = 6 * quarter_size
begin2 = 8 * quarter_size
end2 = 12 * quarter_size
X[:, begin1:end1] = X1
X[:, begin2:end2] = X2
X = X.reshape([-1, 1])
Xtail = np.random.normal(0, std0, [2 * quarter_size, 1])
X = np.concatenate([X, Xtail])
length = len(X)

# Visualize motifs and synthesized time series
fig = plt.figure(figsize=(4, 3))
ax = fig.add_subplot(2, 2, 1)
ax.plot(pattern1, c='tab:blue')
ax.set_xticks(np.arange(0, pattern_size + 1, quarter_size))
ax.set_ylim(-2, 2)
ax.set_yticks([-1, 0, 1])
ax.set_title('Motif 1')
ax = fig.add_subplot(2, 2, 2)
ax.plot(pattern2, c='tab:orange')
ax.set_xticks(np.arange(0, pattern_size + 1, quarter_size))
ax.set_ylim(-2, 2)
ax.set_yticks([-1, 0, 1])
ax.set_title('Motif 2')
ax = fig.add_subplot(2, 1, 2)
ax.plot(X[:26 * quarter_size], c='0.5', lw=1)
ax.plot(np.arange(2 * quarter_size, 6 * quarter_size), X[2 * quarter_size:6 * quarter_size], c='tab:blue')
ax.plot(np.arange(8 * quarter_size, 12 * quarter_size), X[8 * quarter_size:12 * quarter_size], c='tab:orange')
ax.plot(np.arange(14 * quarter_size, 18 * quarter_size), X[14 * quarter_size:18 * quarter_size], c='tab:blue')
ax.plot(np.arange(20 * quarter_size, 24 * quarter_size), X[20 * quarter_size:24 * quarter_size], c='tab:orange')
ax.set_xlabel('Time')
ax.set_title('Synthesized time series')
plt.tight_layout()
plt.savefig('../figures/synthetic_ts.pdf', bbox_inches='tight', pad_inches=0)

# +
# Get ground truth labels
actuals1 = np.full(num_samples, 0)
actuals2 = np.full(num_samples, 1)
actuals = np.concatenate([actuals1, actuals2])

# Get ground truth indices
indices1 = np.arange(2 * quarter_size, length, period)
indices2 = np.arange(8 * quarter_size, length, period)
indices = np.concatenate([indices1, indices2])
# -

# Visualize MDS
mds = motif.MultiDimensionalScaling(markers=['^', 'o'])
mds.fit(np.concatenate([X1std, X2std], axis=0), actuals)
mds.plot(width=4, height=3, filename='../figures/synthetic_mds.pdf')

# +
# ScanMK
mk = motif.ScanMK()
mk.prepare([X], window_size)
mk.scan()

# Enumerate motifs for each radius
mk_purities = {}
for radius in radiuses:
    print('\nradius =', radius)
    mk.find(radius)
    print(mk.motifs)
    if len(mk.motifs) > 0:
        # Visualize MDS
#         predictions = mk.get_nearest_labels(indices)
#         mds.plot(predictions, width=4, height=3,
#                  title='ScanMK (radius = {:g})'.format(radius),
#                  filename='../figures/synthetic_mds_mk_{:.2f}.pdf'.format(radius))
        # Compute purities
        purities = [np.nan]
        for num_motifs in range(1, len(mk.motifs) + 1):
            predictions = mk.get_nearest_labels(indices, num_motifs)
            purity = motif.get_purity(actuals, predictions)
            purities.append(purity)
        mk_purities[radius] = purities
        print(purities)

# +
# SetFinder
sf = motif.SetFinder()
sf.prepare([X], window_size)
sf.count(radiuses)

# Enumerate motifs for each radius
sf_purities = {}
for radius in radiuses:
    print('\nradius =', radius)
    sf.find(radius)
    print(sf.motifs)
    if len(sf.motifs) > 0:
        # Visualize MDS
#         predictions = sf.get_nearest_labels(indices)
#         mds.plot(predictions, width=4, height=3,
#                  title='SetFinder (radius = {:g})'.format(radius),
#                  filename='../figures/synthetic_mds_sf_{:.2f}.pdf'.format(radius))
        # Compute purities
        purities = [np.nan]
        for num_motifs in range(1, len(sf.motifs) + 1):
            predictions = sf.get_nearest_labels(indices, num_motifs)
            purity = motif.get_purity(actuals, predictions)
            purities.append(purity)
        sf_purities[radius] = purities
        print(purities)

# +
# HubFinder
hf = motif.HubFinder()
hf.compute([X], window_size)
print(hf.motifs)

# Compute purities
hf_purities = [np.nan]
for num_motifs in range(1, len(hf.motifs) + 1):
    predictions = hf.get_nearest_labels(indices, num_motifs)
    # Visualize MDS
#     mds.plot(predictions, width=4, height=3,
#              title='HubFinder (number of motifs = {:g})'.format(num_motifs),
#              filename='../figures/synthetic_mds_hf_{:g}.pdf'.format(num_motifs))
    # Compute purities
    purity = motif.get_purity(actuals, predictions)
    hf_purities.append(purity)
print(hf_purities)

# +
# Get radius vs purity for #motifs=2
result_mk_radiuses = []
result_mk_purities = []
for radius in radiuses:
    try:
        result_mk_purities.append(mk_purities[radius][2])
        result_mk_radiuses.append(radius)
    except KeyError:
        continue
    except IndexError:
        result_mk_purities.append(mk_purities[radius][-1])
        result_mk_radiuses.append(radius)
result_sf_radiuses = []
result_sf_purities = []
for radius in radiuses:
    try:
        result_sf_purities.append(sf_purities[radius][2])
        result_sf_radiuses.append(radius)
    except KeyError:
        continue
    except IndexError:
        result_sf_purities.append(sf_purities[radius][-1])
        result_sf_radiuses.append(radius)

# Get the best radius and purity
mk_argmax = np.array(result_mk_purities).argmax()
mk_best_radius = result_mk_radiuses[mk_argmax]
mk_best_purity = result_mk_purities[mk_argmax]
sf_argmax = np.array(result_sf_purities).argmax()
sf_best_radius = result_sf_radiuses[sf_argmax]
sf_best_purity = result_sf_purities[sf_argmax]
hf_purity = hf_purities[2]
print('purity (num_motifs = 2)')
print('ScanMK: {:g} (@ radius = {:g})'.format(mk_best_purity, mk_best_radius))
print('SetFinder: {:g} (@ radius = {:g})'.format(sf_best_purity, sf_best_radius))
print('HubFinder: {:g}'.format(hf_purity))
mk.find(mk_best_radius)
sf.find(sf_best_radius)
# -

# Visualize radius vs purity
fig, ax = plt.subplots(figsize=(4, 3))
ax.plot(result_mk_radiuses, result_mk_purities, '--', c='tab:blue', label='ScanMK')
ax.plot(result_sf_radiuses, result_sf_purities, '-.', c='tab:orange', label='SetFinder')
ax.plot(mk_best_radius, mk_best_purity, 'o', c='tab:blue')
ax.plot(sf_best_radius, sf_best_purity, 'o', c='tab:orange')
ax.axhline(hf_purity, ls='-', c='tab:green', label='HubFinder')
ax.set_xlabel('Radius')
ax.set_ylabel('Purity')
ax.set_ylim(-0.1, 1.1)
ax.legend()
plt.savefig('../figures/synthetic_radius_vs_purity.pdf', bbox_inches='tight', pad_inches=0)

# Visualize #motifs vs purity (ScanMK)
fig, ax = plt.subplots(figsize=(4, 3))
for radius in radiuses:
    try:
        ax.plot(mk_purities[radius], lw=1, c='0.8')
    except KeyError:
        pass
ax.axvline(2, ls=':', c='0.5')
ax.plot(mk_purities[mk_best_radius], '--', c='tab:blue',
        label='ScanMK (radius = {:g})'.format(mk_best_radius))
ax.plot(hf_purities, '-', c='tab:green', label='HubFinder')
ax.plot(2, mk_purities[mk_best_radius][2], 'o', c='tab:blue')
ax.plot(2, hf_purities[2], 'o', c='tab:green')
ax.set_xlabel('Number of motifs')
ax.set_ylabel('Purity')
ax.set_ylim(-0.1, 1.1)
ax.legend()
plt.savefig('../figures/synthetic_motif_vs_purity_mk.pdf', bbox_inches='tight', pad_inches=0)

# Visualize #motifs vs purity (SetFinder)
fig, ax = plt.subplots(figsize=(4, 3))
for radius in radiuses:
    try:
        ax.plot(sf_purities[radius], lw=1, c='0.8')
    except KeyError:
        pass
ax.axvline(2, ls=':', c='0.5')
ax.plot(sf_purities[sf_best_radius], '-.', c='tab:orange',
        label='SetFinder (radius = {:g})'.format(sf_best_radius))
ax.plot(hf_purities, '-', c='tab:green', label='HubFinder')
ax.plot(2, sf_purities[sf_best_radius][2], 'o', c='tab:orange')
ax.plot(2, hf_purities[2], 'o', c='tab:green')
ax.set_xlabel('Number of motifs')
ax.set_ylabel('Purity')
ax.set_ylim(-0.1, 1.1)
ax.legend()
plt.savefig('../figures/synthetic_motif_vs_purity_sf.pdf', bbox_inches='tight', pad_inches=0)

# +
# Visualize motifs
fig = plt.figure(figsize=(4, 3))

# ScanMK
ax = fig.add_subplot(3, 1, 1)
ax.plot(X, c='0.5', lw=1)
for i in np.where(mk.labels[1] == 0)[0]:
    ax.plot(range(i, i + window_size), X[i:i + window_size], c='tab:blue')
for i in np.where(mk.labels[1] == 1)[0]:
    ax.plot(range(i, i + window_size), X[i:i + window_size], c='tab:orange')
ax.set_xlim(20*period, 25*period)
ax.set_ylabel('ScanMK')

# SetFinder
ax = fig.add_subplot(3, 1, 2)
ax.plot(X, c='0.5', lw=1)
for i in np.where(sf.labels[1] == 0)[0]:
    ax.plot(range(i, i + window_size), X[i:i + window_size], c='tab:blue')
for i in np.where(sf.labels[1] == 1)[0]:
    ax.plot(range(i, i + window_size), X[i:i + window_size], c='tab:orange')
ax.set_xlim(20*period, 25*period)
ax.set_ylabel('SetFinder')

# HubFinder
ax = fig.add_subplot(3, 1, 3)
ax.plot(X, c='0.5', lw=1)
for i in np.where(hf.labels[1] == 0)[0]:
    ax.plot(range(i, i + window_size), X[i:i + window_size], c='tab:blue')
for i in np.where(hf.labels[1] == 1)[0]:
    ax.plot(range(i, i + window_size), X[i:i + window_size], c='tab:orange')
ax.set_xlim(20*period, 25*period)
ax.set_ylabel('HubFinder')

# Save
plt.tight_layout()
plt.savefig('../figures/synthetic_motifs.pdf', bbox_inches='tight', pad_inches=0)
# -

# Visualize motif labels
mk.plot_labels(2, xscale=0.0015, yscale=0.5, ylabel='ScanMK', filename='../figures/synthetic_mk_motifs.pdf')
sf.plot_labels(2, xscale=0.0015, yscale=0.5, ylabel='SetFinder', filename='../figures/synthetic_sf_motifs.pdf')
hf.plot_labels(2, xscale=0.0015, yscale=0.5, ylabel='HubFinder', filename='../figures/synthetic_hf_motifs.pdf')



StarcoderdataJupyterScriptsDedupFiltered