enigma_tools.py

# coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as LA
import scipy.optimize as optim


class UnivFitFunc:

    def __init__(self, err=5, Nmax=1000):
        self.err = err
        self.Nmax = Nmax

    def mean_fun(self, X):

        T = np.shape(X)[0]
        x_mean = np.zeros((T,))
        for d in range(T):
            x_mean[d] = np.mean(X[d, :])
        return x_mean

    def __slopes_routine(self, X, flag):

        x_mean = self.mean_fun(X)

        M = np.shape(X)[1]
        T = np.shape(X)[0]
        slope = np.zeros((M,))
        for b in range(M):
            slope[b] = np.dot(X[:, b].reshape(1, len(X[:, b])), x_mean) \
                       / np.dot(x_mean.reshape(1, len(x_mean)), x_mean)
        slope_sort = np.sort(slope)[::-1]

        y_mean = np.mean(slope_sort)
        y_max = np.max(slope_sort)
        y_min = np.min(slope_sort)

        y_up = y_mean + (y_max - y_mean) / 3
        y_dn = y_mean - np.abs(y_min - y_mean) / 3

        Nmn = 0
        Nup = 0
        Ndn = 0

        for b in range(len(slope)):
            if y_up >= slope[b] >= y_dn:
                Nmn = Nmn + 1
            if y_max >= slope[b] > y_up:
                Nup = Nup + 1
            if y_dn > slope[b] >= y_min:
                Ndn = Ndn + 1

        Ndn = M - Ndn

        if flag == 'test':
            Rt = (Nmn / M) * 100
            y = [y_up, y_mean, y_dn]
            sl = [slope, slope_sort]
            return M, y, sl, Rt
        elif flag == 'fit':
            idx_slope_sort = np.argsort(slope)[::-1]
            X = X[:, idx_slope_sort]
            F_up = (1 / Nup) * np.sum(X[:, 0:Nup + 1], axis=1)
            F_mn = (1 / Nmn) * np.sum(X[:, Nup + 1:Ndn + 1], axis=1)
            F_dn = (1 / Ndn) * np.sum(X[:, Ndn + 1:M + 1], axis=1)
            F = [F_up, F_mn, F_dn]
            return M, T, F

    def test_set(self, X):

        M, y, sl, Rt = self.__slopes_routine(X, flag='test')

        print('Reliability: ' + str(round(Rt, 2)) + '%')

        y_up = np.ones((M,)) * y[0]
        y_mean = np.ones((M,)) * y[1]
        y_dn = np.ones((M,)) * y[2]

        slope = sl[0]
        slope_sort = sl[1]

        m = np.array([i for i in range(M)]) + 1

        plt.figure(figsize=(10, 5), dpi=300)
        plt.plot(m, slope_sort, '-o', label='Sorted slopes')
        plt.plot(m, y_mean, label='Mean realization')
        plt.plot(m, y_up, label='Upper bound')
        plt.plot(m, y_dn, label='Lower bound')
        plt.scatter(m, slope, label='Slopes', color="orange")
        plt.title("Slopes and maximal deviations of the signal")
        plt.ylabel('Signal slopes')
        plt.xlabel('Number of experiments')
        plt.grid()
        plt.legend()
        plt.show()

    def fit(self, X):

        x_mean = self.mean_fun(X)
        M, T, F = self.__slopes_routine(X, flag='fit')
        A = np.dot(F[0], LA.pinv([F[2], F[1]]))
        A1 = A[0]
        A0 = A[1]

        # Find the Kappas
        p = [1, -A1, -A0]
        kappa_roots = np.roots(p)

        # Find T period
        T_min = int(T / 2)  # minimal period
        T_max = 2 * T  # maximal period

        P = [i for i in range(T_min, T_max + 1)]  # massive of possible periods
        S = len(P)

        K = 30  # initial number of modes
        J_fit = np.zeros(np.shape(x_mean))
        s = 1

        flag = 0
        N = 0
        while flag == 0 and N != self.Nmax:
            E = np.zeros((1 + 2 * K, T), dtype=complex)
            T_s = T_min + (s / S) * (T_max - T_min)

            """ 
                According to the formulas (16), (18) and (21) [1]:
                    <y(x)> ∼= F(x; K, T_s) = A_0*E_0(x) + sum_(k=1)^(K>>1){ (Ac_k*Ec_k(x) + As_k*Es_k(x)) },
            
                    E_0(x) = [kappa_1(x)]^(x/T_s) + [|kappa_2(x)|]^(x/T_s) * cos(pi*x/T_s),
                    Ec_k = E_0(x)*cos(2*pi*k*x/T_s), 
                    Es_k = E_0(x)*sin(2*pi*k*x/T_s).
            
                Matrix form:
            
                    F = A*E = [ A_0 Ac^(1) As^(1) Ac^(2) As^(2) ].T * [ E_0 Ec^(1) Es^(1) Ec^(2) Es^(2) ]
            
                Dimensions: 
                    A: 1 x (1 + 2K)
                    E: (1 + 2K) x X
                     
            """
            for x in range(1, T + 1):  # change the timeslots

                E_0 = complex(kappa_roots[0]) ** (x / T_s) \
                      + np.cos(np.pi * x / T_s) * complex(np.abs(kappa_roots[1])) ** (x / T_s)

                Ec = E_0 * np.array([np.cos(2 * np.pi * k * (x / T_s))
                                     for k in range(K)], dtype=complex)
                Es = E_0 * np.array([np.sin(2 * np.pi * k * (x / T_s))
                                     for k in range(K)], dtype=complex)

                E[:, x - 1] = [E_0] + list(Ec) + list(Es)  # concatenation

            A = -optim.lsq_linear(E.T, -x_mean).x
            J_fit = np.dot(A.T, E)

            RelErr = (np.std(x_mean - np.real(J_fit)) / (np.mean(abs(x_mean)))) * 100
            if RelErr > self.err and s < S:
                s = s + 1
                N = N + 1
            elif RelErr > self.err and s == S:
                s = 1
                K = K + 1
                N = N + 1
            else:
                flag = 1
        print('Relative Error: ' + str(round(RelErr, 2)) + '%')
        return J_fit


""" Literature

 1. Nigmatullin, R.R., Zhang, W. and Striccoli, D., 2017. 
 “Universal” Fitting Function for Complex Systems: Case of the Short Samplings. 
 Journal of Applied Nonlinear Dynamics, 6(3), pp.427-443.
 
 2. Nigmatullin, R.R., Zhang, W. and Striccoli, D., 2015. 
 General theory of experiment containing reproducible data: The reduction to an ideal experiment. 
 Communications in Nonlinear Science and Numerical Simulation, 27(1-3), pp.175-192.
 
 3. Nigmatullin, R.R., Maione, G., Lino, P., Saponaro, F. and Zhang, W., 2017. 
 The general theory of the Quasi-reproducible experiments: How to describe the measured data of complex systems?. 
 Communications in Nonlinear Science and Numerical Simulation, 42, pp.324-341.

"""