# -*- coding: utf-8 -*-
from __future__ import annotations


"""
HSDC.py
=======

Optimized, source-consistent, and API-compatible implementation of
Dynamic Complexity.

This module replaces an earlier faulty DC.py implementation while preserving
the original call signature and return format:

    C, F, D = dynamic_complexity(_ts_full=None, win=7, range_=None)

For extended diagnostics (turning points, z-transformed series, parameters),
use:

    result = dynamic_complexity_full(_ts_full=None, win=7, range_=None)

For visualization of already computed Dynamic Complexity values, use:

    fig, axes = plot_dynamic_complexity(ts, C)

For a convenience call that calculates and plots in one step, use:

    C, F, D, fig, axes = compute_and_plot_dynamic_complexity(
        _ts_full=ts,
        win=7,
        range_=6,
    )

Purpose
-------
- Methodologically corrected implementation based on a validated VBA reference
- Full API compatibility with DC.py
- Improved runtime performance through optimized distribution calculation
- Optional extended diagnostics without overhead in the standard call
- Optional plotting helpers separated from the numerical core

Design Principles
-----------------
The numerical core and the visualization layer are deliberately separated.

The function dynamic_complexity(...) performs only the calculation of the
Dynamic Complexity series and its two components:

    C = F * D

where

    F : fluctuation intensity
    D : distribution measure
    C : dynamic complexity

The function plot_dynamic_complexity(...) performs no calculation. It only plots
a previously computed C series together with the original time series. This
avoids unnecessary recalculation when C, F, and D have already been computed.

The function compute_and_plot_dynamic_complexity(...) is only a convenience
wrapper. It calculates C, F, and D exactly once and then forwards the result to
plot_dynamic_complexity(...).

Scale Convention
----------------
The parameter range_ is interpreted as the numerical span of the scale, not as
the number of scale levels.

Examples:

    For a 7-point Likert scale with values 1 ... 7:

        range_ = 6

    For a 0-based scale with values 0 ... 6:

        range_ = 6

    For an observed series where no external scale is specified:

        range_ = None

In the latter case, the observed data range is used.

This convention differs from some older implementations where "range" was used
ambiguously to refer either to the number of possible scale values or to the
numerical value span.

Main Functions
--------------
dynamic_complexity(...)
    Standard API-compatible calculation function.
    Returns C, F, D as Python lists.

dynamic_complexity_full(...)
    Extended diagnostic function.
    Returns an HSDCResult object containing C, F, D, turning-point information,
    z-transformed series, and scale parameters.

plot_dynamic_complexity(...)
    Plot-only helper.
    Requires an already computed C series and performs no recalculation.

compute_and_plot_dynamic_complexity(...)
    Convenience wrapper.
    Computes C, F, D once and then plots C.

Return Formats
--------------
dynamic_complexity(...):

    C, F, D

where each object is a list[float].

dynamic_complexity_full(...):

    HSDCResult(
        C,
        F,
        D,
        turning_points,
        z_data,
        z_F,
        z_D,
        z_C,
        scale_min,
        scale_max,
        sr,
        win,
    )

plot_dynamic_complexity(...):

    fig, axes

compute_and_plot_dynamic_complexity(...):

    C, F, D, fig, axes

Notes
-----
- The standard function computes only C, F, and D for performance reasons.
- Z-transformed outputs are available only via dynamic_complexity_full(...).
- Plotting functions import matplotlib locally, so matplotlib is not required
  for purely numerical use of this module.
- Exact numerical reproduction requires identical window width and scale
  parametrization as in the validated reference logic.
- The implementation does not silently shift zero-based data. A scale from
  0 ... 6 is treated as valid data when the corresponding scale range is used.

Disclaimer
----------
This software is provided "as is", without warranty of any kind, express or
implied, including but not limited to the warranties of merchantability,
fitness for a particular purpose, and noninfringement.

The authors assume no liability for any damages or consequences arising from
the use of this software.

This implementation is intended for research and scientific purposes only.
Users are responsible for validating results before applying them in clinical,
diagnostic, or decision-critical contexts.

This implementation does not constitute medical advice or diagnostic guidance.

References
----------
Schiepek, G., & Strunk, G. (2010).
The identification of critical fluctuations and phase transitions in
short-term and coarse-grained time series.
Biological Cybernetics, 102(3), 197–207.
https://doi.org/10.1007/s00422-009-0362-1

Haken, H., & Schiepek, G. (2010).
Synergetik in der Psychologie: Selbstorganisation verstehen und gestalten.
2nd ed. Göttingen: Hogrefe.
"""

__author__ = "H. Joh. Schöller"
__version__ = "2.0.0"
__year__ = "2018/2026"


from dataclasses import dataclass
from functools import lru_cache
from typing import Iterable, Sequence

import numpy as np

ArrayLike1D = Sequence[float] | np.ndarray | Iterable[float]


@dataclass(frozen=True)
class HSDCResult:
    """
    Extended result object returned by dynamic_complexity_full(...).
    
    Attributes
    ----------
    C, F, D : np.ndarray
        Dynamic Complexity and its components.
    turning_points : np.ndarray
        Detected turning-point / trend-length information.
    z_data : np.ndarray
        Z-transformed input data.
    z_F, z_D, z_C : np.ndarray
        Z-transformed component series.
    scale_min, scale_max : float
        Effective scale limits used internally.
    sr : int
        Scale resolution (number of discrete levels).
    win : int
        Window width used for calculation.
    """
    C: np.ndarray
    F: np.ndarray
    D: np.ndarray
    turning_points: np.ndarray
    z_data: np.ndarray
    z_F: np.ndarray
    z_D: np.ndarray
    z_C: np.ndarray
    scale_min: float
    scale_max: float
    sr: int
    win: int


def _as_1d_float_array(data: ArrayLike1D) -> np.ndarray:
    arr = np.asarray(list(data) if not isinstance(data, np.ndarray) else data, dtype=float)
    if arr.ndim != 1:
        raise ValueError("_ts_full must be one-dimensional.")
    if arr.size == 0:
        raise ValueError("_ts_full must not be empty.")
    return arr


def _zscore(x: np.ndarray) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    if x.size == 0:
        return x.copy()
    if x.size == 1:
        return np.zeros_like(x)
    mean = np.mean(x)
    std = np.std(x, ddof=1)
    if not np.isfinite(std) or std == 0:
        return np.zeros_like(x)
    return (x - mean) / std


def _resolve_scale(arr: np.ndarray, range_: float | None) -> tuple[float, float, int]:
    data_min = float(np.min(arr))
    data_max = float(np.max(arr))

    if range_ is None:
        scale_min = data_min
        scale_max = data_max
        span = scale_max - scale_min
    else:
        span = float(range_)
        if span < 0:
            raise ValueError("range_ must be >= 0.")
        scale_min = data_min
        scale_max = scale_min + span
        if data_max - data_min > span + 1e-12:
            raise ValueError("range_ is smaller than the observed value range of the time series.")

    rounded_span = round(span)
    if abs(span - rounded_span) < 1e-12 and rounded_span >= 1:
        sr = int(rounded_span) + 1
    else:
        sr = 100

    return scale_min, scale_max, sr


def _align_valid_to_full(valid_vals: np.ndarray, n: int, win: int) -> np.ndarray:
    out = np.full(n, np.nan, dtype=float)
    if valid_vals.size == 0:
        return out
    start = win // 2
    out[start:start + valid_vals.size] = valid_vals
    return out


@lru_cache(maxsize=None)
def _dist_index_cache(win: int):
    aa_list = []
    bb_list = []
    glv_list = []

    for cc in range(0, win - 1):
        for DD in range(cc + 1, win):
            for aa in range(cc, DD):
                for bb in range(aa + 1, DD + 1):
                    aa_list.append(aa)
                    bb_list.append(bb)
                    glv_list.append(bb - aa)

    aa_idx = np.asarray(aa_list, dtype=np.int32)
    bb_idx = np.asarray(bb_list, dtype=np.int32)
    glv_unit = np.asarray(glv_list, dtype=float)
    ideal_unit_sum = float(np.sum(glv_unit))
    return aa_idx, bb_idx, glv_unit, ideal_unit_sum


def _map_window_to_scale(window: np.ndarray, min_val: float, max_val: float, sr: int) -> np.ndarray:
    gain = sr / (max_val - min_val + 1.0) if (max_val - min_val + 1.0) != 0 else 0.0
    mapped = np.rint(min_val + (window - min_val) * gain)
    mapped = np.minimum(mapped, sr)
    return mapped.astype(float)


def _dist_single_window_vba_fast(window: np.ndarray, min_val: float, max_val: float, sr: int) -> float:
    win = window.size
    sorted_scaled = np.sort(_map_window_to_scale(window, min_val, max_val, sr))
    aa_idx, bb_idx, glv_unit, ideal_unit_sum = _dist_index_cache(win)

    ss = (sr - 1.0) / (win - 1.0)
    glv = glv_unit * ss
    real = sorted_scaled[bb_idx] - sorted_scaled[aa_idx]
    dev = np.where(real <= glv, glv - real, 0.0)

    gl_sum = ideal_unit_sum * ss
    if gl_sum == 0:
        return 0.0
    return float(1.0 - np.sum(dev) / gl_sum)


def _compute_f_valid(arr: np.ndarray, win: int, scale_min: float, scale_max: float) -> tuple[np.ndarray, np.ndarray]:
    n = arr.size
    out_len = n - win + 1
    fvals = np.zeros(out_len, dtype=float)
    turning_points = np.zeros(n, dtype=float)

    strt = int(np.floor(win / 2))
    stp = n - (win - strt)
    maxind = -1

    for ii in range(strt, stp + 1):
        cumval = 0.0
        start_idx = ii - strt
        start_val = arr[start_idx]
        stop_val = start_val
        vz = 0
        sigvz = 0
        seg_n = 0
        end_idx = ii - strt + win - 1

        for jj in range(start_idx + 1, end_idx + 1):
            diff = arr[jj] - arr[jj - 1]
            if diff > 0:
                vz = 1
            elif diff < 0:
                vz = -1
            else:
                vz = sigvz

            if sigvz == 0:
                sigvz = vz

            if vz == sigvz and jj < end_idx:
                stop_val = arr[jj]
                seg_n += 1

            if vz == -sigvz and vz != 0 and jj < end_idx:
                if seg_n > 0:
                    cumval += abs(stop_val - start_val) / seg_n
                if jj > maxind and jj < n:
                    turning_points[jj] = seg_n
                    maxind = jj
                start_val = arr[jj - 1]
                stop_val = arr[jj]
                seg_n = 1
                sigvz = vz

            if jj >= end_idx:
                seg_n += 1
                if vz == sigvz:
                    stop_val = arr[jj]
                else:
                    stop_val = arr[jj - 1]
                    if seg_n == 1:
                        seg_n = 2
                    cumval += abs(stop_val - start_val) / (seg_n - 1)
                    if jj > maxind and jj < n:
                        turning_points[jj] = seg_n - 1
                        maxind = jj
                    start_val = stop_val
                    stop_val = arr[jj]
                    seg_n = 1
                cumval += abs(stop_val - start_val) / seg_n

        kk = ii - strt
        if (scale_max - scale_min) > 0:
            fvals[kk] = cumval / (scale_max - scale_min) / (win - 1)
        else:
            fvals[kk] = 0.0

    return fvals, turning_points


def _compute_d_valid(arr: np.ndarray, win: int, scale_min: float, scale_max: float, sr: int) -> np.ndarray:
    n = arr.size
    out_len = n - win + 1
    dvals = np.zeros(out_len, dtype=float)
    for i in range(out_len):
        dvals[i] = _dist_single_window_vba_fast(arr[i:i + win], scale_min, scale_max, sr)
    return dvals


def dynamic_complexity(
    _ts_full: ArrayLike1D = None,
    win: int = 7,
    range_: float | None = None,
    plt_flag: bool = False,
):
    """
    Compute Dynamic Complexity (API-compatible with DC.py).
    
    Parameters
    ----------
    _ts_full : 1D array-like
        Input time series (list or NumPy array).
    win : int, default=7
        Window width for moving-window computation.
    range_ : float or None, default=None
        Scale range (span of values). For a 7-point Likert scale,
        typically use range_=6. If None, the observed data range is used.
    plt_flag : bool, default=False
        Backward-compatible plotting flag. If True, the already computed
        Dynamic Complexity series is plotted after calculation. This parameter
        is accepted for compatibility with older PTDA/DC call sites.
    
    Returns
    -------
    C : list[float]
        Dynamic Complexity values.
    F : list[float]
        Fluctuation Intensity values.
    D : list[float]
        Distribution measure values.
    
    Notes
    -----
    - This function preserves the original DC.py interface and accepts
      plt_flag for backward compatibility with older call sites.
    - Only C, F, D are computed for efficiency.
    - If plt_flag=True, no second calculation is performed; the already
      computed C series is forwarded to plot_dynamic_complexity(...).
    - Use dynamic_complexity_full(...) for extended diagnostics.
    """
    
    if _ts_full is None:
        raise ValueError("_ts_full must not be None.")

    arr = _as_1d_float_array(_ts_full)
    if win < 2:
        raise ValueError("win must be >= 2.")
    if arr.size < win:
        return [], [], []

    scale_min, scale_max, sr = _resolve_scale(arr, range_)

    f_valid, _turning_points = _compute_f_valid(arr, win, scale_min, scale_max)
    d_valid = _compute_d_valid(arr, win, scale_min, scale_max, sr)
    c_valid = f_valid * d_valid

    if plt_flag:
        plot_dynamic_complexity(arr, c_valid)

    return c_valid.tolist(), f_valid.tolist(), d_valid.tolist()


def dynamic_complexity_full(_ts_full: ArrayLike1D = None, win: int = 7, range_: float | None = None) -> HSDCResult:

    """
    Compute Dynamic Complexity with extended diagnostics.

    Parameters
    ----------
    _ts_full : 1D array-like
        Input time series.
    win : int, default=7
        Window width.
    range_ : float or None, default=None
        Scale range. If None, the observed data range is used.

    Returns
    -------
    HSDCResult
        Object containing:
        - C, F, D
        - turning_points
        - z-transformed series
        - scale parameters

    Notes
    -----
    - Includes z-transformed outputs (computed with ddof=1).
    - Constant series are mapped to zero z-scores.
    """

    if _ts_full is None:
        raise ValueError("_ts_full must not be None.")

    arr = _as_1d_float_array(_ts_full)
    if win < 2:
        raise ValueError("win must be >= 2.")
    if arr.size < win:
        raise ValueError("The time series is shorter than the window width.")

    scale_min, scale_max, sr = _resolve_scale(arr, range_)

    f_valid, turning_points = _compute_f_valid(arr, win, scale_min, scale_max)
    d_valid = _compute_d_valid(arr, win, scale_min, scale_max, sr)
    c_valid = f_valid * d_valid

    z_data = _zscore(arr)
    z_F = _zscore(_align_valid_to_full(f_valid, arr.size, win))
    z_D = _zscore(_align_valid_to_full(d_valid, arr.size, win))
    z_C = _zscore(_align_valid_to_full(c_valid, arr.size, win))

    return HSDCResult(
        C=c_valid,
        F=f_valid,
        D=d_valid,
        turning_points=turning_points,
        z_data=z_data,
        z_F=z_F,
        z_D=z_D,
        z_C=z_C,
        scale_min=scale_min,
        scale_max=scale_max,
        sr=sr,
        win=win,
    )

# ---------------------------------------------------------------------------
# Plotting helpers
# ---------------------------------------------------------------------------
# These functions are intentionally separated from the numerical core.
# plot_dynamic_complexity(...) does not recalculate C, F, or D.
# compute_and_plot_dynamic_complexity(...) is a convenience wrapper and
# calculates C, F, and D exactly once.

def plot_dynamic_complexity(
    _ts_full: ArrayLike1D,
    C: ArrayLike1D,
    *,
    ax=None,
    title_original: str = "Original Time Series",
    title_dc: str = "Dynamic Complexity",
    show: bool = True,
    pause: float | None = 0.05,
):
    """
    Plot an already computed Dynamic Complexity series.

    This function performs no Dynamic Complexity calculation. It only visualizes
    the original time series and a previously computed C series.

    Parameters
    ----------
    _ts_full : 1D array-like
        Original input time series.
    C : 1D array-like
        Previously computed Dynamic Complexity values.
    ax : None or sequence of matplotlib axes, optional
        If None, a new figure with two subplots is created.
        If provided, it must contain two axes: ax[0] for the original time series
        and ax[1] for Dynamic Complexity.
    title_original : str, default="Original Time Series"
        Title of the upper plot.
    title_dc : str, default="Dynamic Complexity"
        Title of the lower plot.
    show : bool, default=True
        Whether to show/update the plot.
    pause : float or None, default=0.05
        Pause duration for interactive environments. If None, no pause is used.

    Returns
    -------
    fig, axes
        Matplotlib figure and axes objects.

    Notes
    -----
    - No calculation is performed in this function.
    - Use this function when C has already been calculated with
      dynamic_complexity(...) or dynamic_complexity_full(...).
    """

    import matplotlib.pyplot as plt

    ts = _as_1d_float_array(_ts_full)
    c_arr = _as_1d_float_array(C)

    if ax is None:
        fig, axes = plt.subplots(2, 1)
    else:
        axes = ax
        if len(axes) != 2:
            raise ValueError("ax must contain exactly two matplotlib axes.")
        fig = axes[0].figure

    axes[0].plot(ts)
    axes[0].set_title(title_original)

    axes[1].plot(c_arr)
    axes[1].set_title(title_dc)

    fig.tight_layout()

    if show:
        plt.show(block=False)
        if pause is not None:
            plt.pause(pause)

    return fig, axes


def compute_and_plot_dynamic_complexity(
    _ts_full: ArrayLike1D = None,
    win: int = 7,
    range_: float | None = None,
    *,
    ax=None,
    title_original: str = "Original Time Series",
    title_dc: str = "Dynamic Complexity",
    show: bool = True,
    pause: float | None = 0.05,
):
    """
    Compute Dynamic Complexity and plot the result.

    This is a convenience wrapper around dynamic_complexity(...) and
    plot_dynamic_complexity(...).

    Parameters
    ----------
    _ts_full : 1D array-like
        Input time series.
    win : int, default=7
        Window width for moving-window computation.
    range_ : float or None, default=None
        Scale range passed to dynamic_complexity(...).
    ax : None or sequence of matplotlib axes, optional
        If None, a new figure with two subplots is created.
        If provided, it must contain two axes.
    title_original : str, default="Original Time Series"
        Title of the upper plot.
    title_dc : str, default="Dynamic Complexity"
        Title of the lower plot.
    show : bool, default=True
        Whether to show/update the plot.
    pause : float or None, default=0.05
        Pause duration for interactive environments. If None, no pause is used.

    Returns
    -------
    C, F, D, fig, axes
        Dynamic Complexity values, component values, and matplotlib plot objects.

    Notes
    -----
    - This function calculates C, F, D exactly once.
    - If C, F, D have already been calculated, use plot_dynamic_complexity(...)
      instead to avoid recalculation.
    """

    C, F, D = dynamic_complexity(_ts_full=_ts_full, win=win, range_=range_)

    fig, axes = plot_dynamic_complexity(
        _ts_full,
        C,
        ax=ax,
        title_original=title_original,
        title_dc=title_dc,
        show=show,
        pause=pause,
    )

    return C, F, D, fig, axes


__all__ = [
    "dynamic_complexity",
    "dynamic_complexity_full",
    "HSDCResult",
    "plot_dynamic_complexity",
    "compute_and_plot_dynamic_complexity",
]
