kernel-adaptive-filtering/lorenz.py at master · pin3da/kernel-adaptive-filtering · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
import time
import contextlib

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint

from filters import (APA1, APA2, APA3, APA4, CKLMS, KAPA1, KAPA2, KAPA3, KAPA4,
                     KLMS, KRLS, LMS, QKLMS, get_training_error)


def gen_lorenz_data():
    rho = 28.0
    sigma = 10.0
    beta = 8.0 / 3.0

    def f(state, t):
        x, y, z = state  # unpack the state vector
        return sigma * (y - x), x * (rho - z) - y, x * y - beta * z  # derivatives

    state0 = [1.0, 1.0, 1.0]
    t = np.arange(0.0, 40.0, 0.01)

    states = odeint(f, state0, t)
    return states


def split_training_and_testing(states):
    # time embedding
    TD = 5

    # dataset size
    N_tr = 250
    N_te = 10

    # train data
    X = np.zeros((N_tr, TD))
    for k in range(N_tr):
        X[k, :] = states[k: k + TD]

    # test data
    X_te = np.zeros((N_te, TD))
    for k in range(N_te):
        X_te[k, :] = states[k + N_tr: k + TD + N_tr]

    # Desired signal
    T = np.zeros((N_tr, 1))
    for i in range(N_tr):
        T[i, 0] = states[i]

    # Desire signal testing
    T_te = np.zeros((N_te, 1))
    for i in range(N_te):
        T_te[i, 0] = states[i + N_tr]

    return X, X_te, T, T_te, TD


def plot_lorenz(states):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot(states[:, 0], states[:, 1], states[:, 2])
    plt.show()


@contextlib.contextmanager
def time_it(name):
    start = time.monotonic()
    print(name, end=' ')
    yield
    print('took %f seconds' % (time.monotonic() - start))


if __name__ == '__main__':
    states = gen_lorenz_data()
    # plot_lorenz(states)

    states = np.array([
        x for x, y, z in states
    ])

    X, X_te, T, T_te, TD = split_training_and_testing(states)
    filters = [
        LMS(TD, 0.0001),
        APA1(X[0], T[0], 10, 0.0001),
        APA2(X[0], T[0], 10, 0.02),
        APA3(X[0], T[0], 10, 0.0001, 0.01),
        APA4(X[0], T[0], 10, 0.005),
        # KLMS(TD, X[0], T[0], 0.034, 20),
        # CKLMS(X[0], T[0], 0.03, 10, 10),
        # QKLMS(TD, X[0], T[0], 0.03, 1, 10),
        # KAPA1(X[0], T[0], 10, 0.3, 5),
        # KAPA2(X[0], T[0], 0.2, 2.25, 10),
        # KAPA3(X[0], T[0], 0.1, 2.25, 10),
        # KAPA4(X[0], T[0], 0.2, 2.25, 10),
        # KRLS(X[0], T[0], 1, 1),
    ]
    for fi in filters:
        with time_it(fi.name()):
            err = get_training_error(fi, X, X_te, T, T_te, TD)
        plt.plot(err[50:], label=fi.name())

    plt.legend()
    plt.ylabel('MSE')
    plt.xlabel('iteration')
    plt.title('Lorenz oscillator')
    # plt.savefig('./compare-lorenz-2.png')
    plt.show()