# Authors:
#           Wilson Rocha Lacerda Junior <wilsonrljr@outlook.com>
#           Luan Pascoal da Costa Andrade <luan_pascoal13@hotmail.com>
#           Samuel Carlos Pessoa Oliveira <samuelcpoliveira@gmail.com>
#           Samir Angelo Milani Martins <martins@ufsj.edu.br>
# License: BSD 3 clause


import numpy as np
import matplotlib.pyplot as plt
from ..utils.deprecation import deprecated


def compute_residues_autocorrelation(y, yhat):
    e = calculate_residues(y, yhat)
    unnormalized_e_acf = get_unnormalized_e_acf(e)
    half_of_symmetry_autocorr = int(np.floor(unnormalized_e_acf.size / 2))

    e_acf = (
        unnormalized_e_acf[half_of_symmetry_autocorr:]
        / unnormalized_e_acf[half_of_symmetry_autocorr]
    )

    half_of_symmetry_autocorr = int(np.floor(unnormalized_e_acf.size / 2))
    n = len(unnormalized_e_acf) - half_of_symmetry_autocorr
    upper_bound = 1.96 / np.sqrt(len(unnormalized_e_acf))
    lower_bound = upper_bound * (-1)
    return e_acf, upper_bound, lower_bound


def calculate_residues(y, yhat):
    return (y - yhat).flatten()


def get_unnormalized_e_acf(e):
    return np.correlate(e, e, mode="full")


def compute_cross_correlation(y, yhat, arr):
    e = calculate_residues(y, yhat)
    n = len(e) * 2 - 1
    ccf, upper_bound, lower_bound = _input_ccf(e, arr, n)
    return ccf, upper_bound, lower_bound


def _input_ccf(e, a, n):
    ccf = _normalized_correlation(a, e)
    upper_bound = 1.96 / np.sqrt(n)
    lower_bound = upper_bound * (-1)
    return ccf, upper_bound, lower_bound


def _normalized_correlation(a, b):
    """Compute the normalized correlation between two signals.

    Parameters
    ----------
    a : array-like of shape = n_samples.
    b : array-like of shape = n_samples.

    Returns
    -------
    ruy : ndarray of floats:
        The normalized cross correlation between the two signals.

    """
    y = (a - np.mean(a)).flatten()
    u = (b - np.mean(b)).flatten()
    t = int(np.floor(len(a) / 2))
    ruy = np.array(np.zeros(t))
    ruy[0] = np.sum(y * u) / (np.sqrt(np.sum(y ** 2)) * np.sqrt(np.sum(u ** 2)))

    for i in range(1, t):
        y = (a - np.mean(a[:i])).flatten()
        u = (b - np.mean(b[i:])).flatten()
        ruy[i] = np.sum(y[:-i] * u[i:]) / (
            np.sqrt(np.sum(y[:-i] ** 2)) * np.sqrt(np.sum(u[i:] ** 2))
        )

    return ruy


class ResiduesAnalysis:
    """Residues analysis for Polynomial NARX model."""

    @deprecated(
        version="v0.1.7",
        future_version="v0.2.0",
        alternative="from sysidentpy.residues_correlation import compute_cross_correlation, compute_residues_autocorrelation",
    )
    def residuals(self, X, y, yhat):
        """Performs the residual analysis of output to validate model.

        Parameters
        ----------
        y : array-like of shape = n_samples
            The target data used in the identification process.
        yhat : array-like of shape = n_samples
            The prediction values of the identification process.
        X : ndarray of floats
            The input data.

        Returns
        -------
        output_autocorr : ndarray of floats:
            1st column - Residuals normalized autocorrelation.
            2nd/3rd columns - Superior and inferior limits of a
            95% confidence interval.
        output_crosscorr : ndarray of floats:
            1st column - Correlation between residuals and input.
            2nd/3rd columns - Superior and inferior limits of a
            95% confidence interval.

        Examples
        --------
        >>> y = [3, -0.5, 2, 7]
        >>> autocorr(y)
        [62.25 11.5   2.5  21.  ]

        """
        e, e_acf, unnormalized_e_acf = self._residuals_acf(y, yhat)

        xe_ccf = self._input_ccf(X[:, 0], e, len(unnormalized_e_acf))
        ex = e * X[:, 0]
        ee = e ** 2
        x2line = (X[:, 0] ** 2) - (X[:, 0] ** 2).mean()
        eeline = (e ** 2) - (e ** 2).mean()
        e_ex = self._input_ccf(e, ex, len(unnormalized_e_acf))
        x2line_e = self._input_ccf(x2line, e, len(unnormalized_e_acf))
        x2line_ee = self._input_ccf(x2line, ee, len(unnormalized_e_acf))
        ye = y * e.reshape(-1, 1)
        yeline = (y * e.reshape(-1, 1)) - (y * e.reshape(-1, 1)).mean()

        yeline_x2line = self._input_ccf(yeline, x2line, len(unnormalized_e_acf))
        lam = np.sqrt(sum((ee - ee.mean()) ** 2) / sum((ye - ye.mean()) ** 2))
        yeline_eeline = self._input_ccf(yeline, eeline, len(unnormalized_e_acf))
        return (
            e_acf,
            xe_ccf,
            [yeline_x2line, yeline_eeline, e_ex, x2line_e, x2line_ee],
            lam,
        )

    def _input_ccf(self, X, e, len_confidence):
        xe_ccf = np.zeros((int(np.floor(len(X) / 2)), 3), dtype=float)

        xe_ccf[:, 0] = self._normalized_correlation(X, e)

        xe_ccf[:, 1] = np.ones(len(xe_ccf[:, 0])) * (1.96 / np.sqrt(len_confidence))

        xe_ccf[:, 2] = xe_ccf[:, 1] * (-1)
        return xe_ccf

    def _residuals_acf(self, y, yhat):
        e = (y - yhat).flatten()
        unnormalized_e_acf = np.correlate(e, e, mode="full")
        half_of_simmetry_autocorr = int(np.floor(unnormalized_e_acf.size / 2))

        e_acf = np.zeros(
            (len(unnormalized_e_acf) - half_of_simmetry_autocorr, 3), dtype=float
        )

        e_acf[:, 0] = (
            unnormalized_e_acf[half_of_simmetry_autocorr:]
            / unnormalized_e_acf[half_of_simmetry_autocorr]
        )

        e_acf[:, 1] = np.ones(len(e_acf[:, 0])) * (
            1.96 / np.sqrt(len(unnormalized_e_acf))
        )

        e_acf[:, 2] = e_acf[:, 1] * (-1)
        return e, e_acf, unnormalized_e_acf

    def _normalized_correlation(self, signal1, signal2):
        """Compute the normalized correlation between two signals.

        Parameters
        ----------
        signal1 : array-like of shape = n_samples.
        signal2 : array-like of shape = n_samples.

        Returns
        -------
        ruy : ndarray of floats:
            The normalized cross correlation between the two signals.

        """
        y = (signal1 - np.mean(signal1)).flatten()
        u = (signal2 - np.mean(signal2)).flatten()
        t = int(np.floor(len(signal1) / 2))
        ruy = np.array(np.zeros(t))
        ruy[0] = np.sum(y * u) / (np.sqrt(np.sum(y ** 2)) * np.sqrt(np.sum(u ** 2)))

        for i in range(1, t):
            y = (signal1 - np.mean(signal1[:i])).flatten()
            u = (signal2 - np.mean(signal2[i:])).flatten()
            ruy[i] = np.sum(y[:-i] * u[i:]) / (
                np.sqrt(np.sum(y[:-i] ** 2)) * np.sqrt(np.sum(u[i:] ** 2))
            )

        return ruy

    def plot_result(self, y, yhat, e_acf, xe_ccf, figsize=(10, 8), n=100):
        """Plot the free run simulation and residues analysis.

        Parameters
        ----------
        y : array-like of shape = n_samples
            The target data used in the identification process.
        yhat : array-like of shape = n_samples
            The prediction values of the identification process.
        e_acf : ndarray of floats:
            1st column - Residuals normalized autocorrelation.
            2nd/3rd columns - Superior and inferior limits of a
            95% confidence interval.
        xe_ccf : ndarray of floats:
            1st column - Correlation between residuals and input.
            2nd/3rd columns - Superior and inferior limits of a
            95% confidence interval.

        """
        plt.style.use("seaborn-white")
        plt.rcParams["axes.facecolor"] = "white"

        fig, axes = plt.subplots(nrows=2, ncols=2, figsize=figsize, facecolor="white")
        fig.subplots_adjust(hspace=0.7)
        for ax, feature in zip(axes.flatten()[2:], [e_acf, xe_ccf]):
            ax.plot(feature[:, 0], color="#1f77b4")
            ax.axhspan(feature[0, 1], feature[0, 2], color="#ccd9ff", alpha=0.5, lw=0)
            ax.set_xlabel("Lag", fontsize=12)
            ax.set_ylabel("Cross Correlation: ee, ex", fontsize=12)
            # ax = plt.gca()
            ax.set_ylim([-1, 1])
            # ax.grid(color="grey", linestyle="-.", alpha=0.1)
            ax.spines["right"].set_visible(False)
            ax.spines["top"].set_visible(False)
            ax.xaxis.set_ticks_position("bottom")
            ax.yaxis.set_ticks_position("left")
            fig.tight_layout()

        ax = plt.subplot(211)
        ax.plot(
            y[self.max_lag : n],
            c="#1f77b4",
            alpha=1,
            marker="o",
            label="Data",
            linewidth=1.5,
        )
        ax.plot(
            yhat[self.max_lag : n],
            c="#ff7f0e",
            marker="*",
            # linestyle='dashed',
            label="Model",
            linewidth=1.5,
        )
        ax.set_title("Free run simulation", fontsize=18)
        ax.legend()
        ax.tick_params(labelsize=14)
        ax.set_xlabel("Samples", fontsize=14)
        ax.set_ylabel("y, yhat", fontsize=14)
        # ax.grid(color="grey", linestyle="-.", alpha=0.1)
        ax.spines["right"].set_visible(False)
        ax.spines["top"].set_visible(False)
        ax.xaxis.set_ticks_position("bottom")
        ax.yaxis.set_ticks_position("left")
        fig.tight_layout()
        plt.show()
