"""Causal five-minute price decomposition with a frozen yearly reference.

Dollar differences, never log returns: zero and negative NEM prices are valid.
Calendar coefficients and empirical thresholds only use earlier observations. Online
states use past observations; realized energy includes the observation at t.
Regimes are descriptive percentile bands, not inferred economic market states.
"""
from __future__ import annotations

from datetime import datetime, timedelta
import json
import math
import os
from pathlib import Path
import tempfile

import numpy as np

VERSION = 'seasonal_innovations_v4'
STEP = 300
DAY = 86400
MIN_TRAIN = 28 * 288
COLUMNS = ['time', 'price', 'seasonal', 'level', 'carry', 'ordinary', 'shock',
           'innovation', 'expected_scale', 'total_rms', 'routine_rms', 'shock_rms',
           'volatility_percentile', 'surprise_percentile', 'shock_tier']
COLUMNS.append('historical_shock_risk_percent')
COLUMNS.append('shock_carry')


def features(seconds, annual=True):
    # NEM settlement timestamps are fixed AEST, including in daylight saving.
    days = np.asarray(seconds, dtype=float) / DAY
    daily = 2 * np.pi * days
    weekday = (np.floor(days).astype(np.int64) + 3) % 7  # 1970-01-01 Thursday
    cols = [np.ones(len(days))]
    cols += [(weekday == day).astype(float) for day in range(1, 7)]
    for k in range(1, 9):
        cols += [np.sin(k * daily), np.cos(k * daily)]
    weekend = (weekday >= 5).astype(float)
    for k in range(1, 5):
        cols += [weekend * np.sin(k * daily), weekend * np.cos(k * daily)]
    if annual:
        angle = 2 * np.pi * days / 365.2425
        cols += [np.sin(angle), np.cos(angle), np.sin(2 * angle), np.cos(2 * angle)]
        for k in (1, 2, 3):
            for d in (np.sin(k * daily), np.cos(k * daily)):
                cols += [d * np.sin(angle), d * np.cos(angle)]
    return np.column_stack(cols)


def predict(seconds, model):
    coefficients = np.asarray(model['coefficients'])
    return np.concatenate([features(seconds[i:i + 8192], model['annual']) @ coefficients
                           for i in range(0, len(seconds), 8192)]) if len(seconds) else np.array([])


def robust_calendar(seconds, prices):
    """Tukey IRLS, bounded influence; chunked sufficient statistics bound RAM."""
    annual = seconds[-1] - seconds[0] >= 540 * DAY
    width = features(seconds[:1], annual).shape[1]
    beta = np.zeros(width); beta[0] = np.median(prices)
    scale = max(1.0, 1.4826 * np.median(np.abs(prices - beta[0])))
    for iteration in range(10):
        matrix = np.zeros((width, width)); target = np.zeros(width)
        for left in range(0, len(prices), 8192):
            x = features(seconds[left:left + 8192], annual)
            y = prices[left:left + 8192]
            u = (y - x @ beta) / (4.685 * scale)
            weights = np.maximum(0, 1 - u * u) ** 2
            matrix += x.T @ (weights[:, None] * x)
            target += x.T @ (weights * y)
        penalty = np.eye(width) * max(1e-8, len(prices) * .0001)
        penalty[0, 0] = 1e-8
        candidate = np.linalg.solve(matrix + penalty, target)
        model = {'annual': bool(annual), 'coefficients': candidate.tolist()}
        residual = prices - predict(seconds, model)
        scale = max(.01, 1.4826 * np.median(np.abs(residual - np.median(residual))))
        converged = np.max(np.abs(candidate - beta)) < .001
        beta = candidate
        if iteration >= 2 and converged:
            break
    return {'annual': bool(annual), 'coefficients': beta.tolist(), 'fit_scale': float(scale)}


def rolling_energy(seconds, innovations, shock_mask, window=12):
    """Exactly additive variance energy; require twelve consecutive valid pairs."""
    n = len(seconds)
    total = np.full(n, np.nan); routine = total.copy(); shock = total.copy()
    if n < window:
        return total, routine, shock
    valid = np.isfinite(innovations)
    squares = np.where(valid, innovations, 0) ** 2
    cumulative = np.r_[0., np.cumsum(squares)]
    jumps = np.r_[0., np.cumsum(np.where(shock_mask, squares, 0))]
    counts = np.r_[0, np.cumsum(valid)]
    right = np.arange(window, n + 1); left = right - window
    complete = ((counts[right] - counts[left] == window) &
                (seconds[right - 1] - seconds[left] == (window - 1) * STEP))
    indexes = right[complete] - 1
    variance = np.maximum(0, (cumulative[right[complete]] - cumulative[left[complete]]) / window)
    jump_variance = np.maximum(0, (jumps[right[complete]] - jumps[left[complete]]) / window)
    total[indexes] = np.sqrt(variance)
    shock[indexes] = np.sqrt(jump_variance)
    routine[indexes] = np.sqrt(np.maximum(0, variance - jump_variance))
    return total, routine, shock


def filter_path(seconds, prices, seasonal, model, components=False):
    """Pre-observation local level + bounded AR carry + innovation.

    A jump's influence on the next interval's residual carry is bounded. Sustained shifts can
    move the robust local level, while isolated spikes have bounded influence.
    """
    n = len(prices)
    level = np.zeros(n); carry = np.zeros(n); innovation = np.full(n, np.nan)
    scale = np.full(n, np.nan); shock_carry = np.zeros(n)
    state_level = 0.; previous_clean = 0.; previous_shock = 0.; variance_ratio = 1.
    phi = model['phi']; base_scale = np.asarray(model['slot_scale'])
    level_rate = 1 - 2 ** (-1 / (288 * 7))
    variance_rate = 1 - 2 ** (-1 / 12)
    slots = (seconds // STEP) % 288
    for i in range(n):
        connected = i > 0 and seconds[i] - seconds[i - 1] == STEP
        if not connected:
            previous_clean = 0.; previous_shock = 0.; variance_ratio = 1.
        level[i] = state_level
        carry[i] = phi * previous_clean if connected else 0.
        shock_carry[i] = model.get('jump_phi', 0.) * previous_shock if connected else 0.
        background_sigma = max(.01, float(base_scale[slots[i]]) * math.sqrt(variance_ratio))
        sigma = math.hypot(background_sigma, model.get('jump_scale', 0.) * previous_shock)
        scale[i] = sigma
        residual = float(prices[i] - seasonal[i] - state_level)
        surprise = float(residual - carry[i] - shock_carry[i])
        if abs(surprise) < 1e-8:
            surprise = 0.
        if connected:
            innovation[i] = surprise
        # Only earlier information enters this observation's expectation/scale.
        level_limit = 4 * background_sigma
        state_level += level_rate * min(level_limit, max(-level_limit, residual))
        # Bound the *new innovation*, not the residual level: an ordinary but
        # persistent displacement from the seasonal curve should remain in carry.
        limit = 4 * background_sigma
        previous_clean = float(carry[i]) + min(limit, max(-limit, float(residual - carry[i])))
        previous_shock = residual - previous_clean
        z = surprise / sigma
        variance_ratio = max(.04, min(100., (1 - variance_rate) * variance_ratio +
                                     variance_rate * min(z * z, 9.) * variance_ratio))
    return (level, carry, innovation, scale, shock_carry) if components else (level, carry + shock_carry, innovation, scale)


def rank(values, reference):
    reference = np.asarray(reference)
    result = np.full(len(values), np.nan)
    valid = np.isfinite(values)
    if len(reference):
        # Midrank handles the many exact zero/unchanged settlement prices.
        result[valid] = 50 * (np.searchsorted(reference, values[valid], side='left') +
                              np.searchsorted(reference, values[valid], side='right')) / len(reference)
    return result


def fit_model(seconds, prices, year):
    start = int((datetime(year, 1, 1) - datetime(1970, 1, 1)).total_seconds())
    # Bootstrap from eight weeks of same-year history when earlier years are unavailable.
    start = max(start, int(seconds[0]) + 57 * DAY)
    year_end = int((datetime(year + 1, 1, 1) - datetime(1970, 1, 1)).total_seconds())
    if start >= year_end:
        return None
    history = (seconds < start) & (seconds >= start - 731 * DAY)
    t, y = seconds[history], prices[history]
    if len(y) < MIN_TRAIN or t[-1] - t[0] < 28 * DAY:
        return None
    holdout_start = start - 28 * DAY
    assessment_start = (datetime(1970, 1, 1) + timedelta(seconds=start)).isoformat()
    train = t < holdout_start
    # At least four weeks for fitting plus four untouched weeks for validation.
    if np.count_nonzero(train) < MIN_TRAIN:
        return None
    model = robust_calendar(t[train], y[train])
    model['assessment_start'] = assessment_start
    season = predict(t, model)
    residual = y - season
    connected = np.r_[False, np.diff(t) == STEP]
    # Bound spikes before estimating persistence; estimate on training only.
    centered = residual - np.median(residual[train])
    bounded = np.clip(centered, -4 * model['fit_scale'], 4 * model['fit_scale'])
    pairs = connected & train
    x = np.roll(bounded, 1)[pairs]; target = bounded[pairs]
    model['phi'] = float(np.clip(np.dot(x, target) / max(1e-9, np.dot(x, x)), 0, .995))
    raw_innovations = residual - model['phi'] * np.roll(residual, 1)
    slots = (t // STEP) % 288
    global_scale = max(.01, 1.4826 * np.median(np.abs(raw_innovations[pairs])))
    model['slot_scale'] = [float(max(.01, 1.4826 * np.median(np.abs(raw_innovations[pairs & (slots == slot)]))))
                           if np.count_nonzero(pairs & (slots == slot)) >= 14 else global_scale
                           for slot in range(288)]
    # A second, fast residual channel predicts continuation/reversion after an
    # observed extreme displacement. Never label a repeated plateau as a fresh
    # unpredictable arrival merely because it remains far from the calendar.
    previous = np.roll(centered, 1)
    extreme_level = max(4 * model['fit_scale'], float(np.quantile(np.abs(centered[train]), .995)))
    extreme_pairs = pairs & (np.abs(previous) > extreme_level)
    if np.count_nonzero(extreme_pairs) >= 20:
        ratios = centered[extreme_pairs] / previous[extreme_pairs]
        weights = np.abs(previous[extreme_pairs]); ordering = np.argsort(ratios)
        cumulative = np.cumsum(weights[ordering])
        median_index = int(np.searchsorted(cumulative, cumulative[-1] / 2))
        model['jump_phi'] = float(np.clip(ratios[ordering[median_index]], 0., 1.))
        deviations = np.abs(ratios - model['jump_phi']); ordering = np.argsort(deviations)
        cumulative = np.cumsum(weights[ordering])
        median_index = int(np.searchsorted(cumulative, cumulative[-1] / 2))
        model['jump_scale'] = float(max(.1, 1.4826 * deviations[ordering[median_index]]))
    else:
        model['jump_phi'] = 0.; model['jump_scale'] = 0.
    # Select persistence on a separate block inside the training period. The
    # final four-week holdout below remains untouched by this selection.
    selection_start = holdout_start - 28 * DAY
    context = train & (t >= selection_start - 90 * DAY)
    selection = t[context] >= selection_start
    scores = []
    candidates = sorted(set([0., .5, .8, .9, .95, .98, .995, 1., model['phi']]))
    for phi in candidates:
        candidate = dict(model, phi=phi)
        _, _, candidate_eps, _ = filter_path(t[context], y[context], season[context], candidate)
        valid = selection & np.isfinite(candidate_eps)
        scores.append((float(np.mean(np.abs(candidate_eps[valid]))), phi))
    model['phi'] = min(scores)[1]
    model['selection'] = {'start': (datetime(1970, 1, 1) + timedelta(seconds=selection_start)).isoformat(),
                          'end': (datetime(1970, 1, 1) + timedelta(seconds=holdout_start)).isoformat(),
                          'candidate_mae': [[float(phi), float(score)] for score, phi in scores]}
    _, _, eps, sigma = filter_path(t, y, season, model)
    z = np.abs(eps / sigma)
    # Threshold calibration excludes the temporal holdout as well.
    reference_z = np.sort(z[train & np.isfinite(z)])
    # Empirical rarity is not a Gaussian p-value or a proof of unpredictability.
    threshold = max(4., float(np.quantile(reference_z, .995)))
    severe = max(threshold, float(np.quantile(reference_z, .999)))
    model['shock_threshold'] = threshold; model['severe_threshold'] = severe
    weekday_hour = ((t // DAY + 3) % 7) * 24 + (t % DAY) // 3600
    historical_shocks = z > threshold
    base_rate = float(np.mean(historical_shocks[train & np.isfinite(z)]))
    model['shock_risk'] = [float(100 * (np.count_nonzero(historical_shocks & train & (weekday_hour == slot)) +
                                100 * base_rate) / (np.count_nonzero(train & np.isfinite(z) & (weekday_hour == slot)) + 100))
                           for slot in range(168)]
    model['overall_shock_risk_percent'] = 100 * base_rate
    total, _, _ = rolling_energy(t, eps, z > threshold)
    reference_vol = np.sort(total[train & np.isfinite(total)])
    # Compact monotone empirical quantile grids: percentile resolution ~0.05pp.
    probabilities = np.linspace(0, 1, 2001)
    model['reference_z'] = np.quantile(reference_z, probabilities).tolist()
    model['reference_vol'] = np.quantile(reference_vol, probabilities).tolist()
    model['reference_start'] = (datetime(1970, 1, 1) + timedelta(seconds=int(t[0]))).isoformat()
    model['reference_end'] = (datetime(1970, 1, 1) + timedelta(seconds=int(t[train][-1]))).isoformat()
    model['reference_count'] = int(np.count_nonzero(train))
    model['reference_quantiles'] = {str(p): float(np.quantile(reference_vol, p / 100))
                                    for p in (25, 50, 75, 90, 99)}
    test = ~train & np.isfinite(eps)
    persistence_error = np.r_[np.nan, np.diff(y)]
    usable = test & connected
    mae = float(np.mean(np.abs(eps[usable])))
    naive_mae = float(np.mean(np.abs(persistence_error[usable])))
    model['validation'] = {'start': (datetime(1970, 1, 1) + timedelta(seconds=holdout_start)).isoformat(),
        'end': assessment_start, 'observations': int(usable.sum()),
        'model_mae': mae, 'persistence_mae': naive_mae,
        'skill_percent': 100 * (1 - mae / naive_mae) if naive_mae > 0 else None,
        'shock_rate_percent': float(100 * np.mean(z[usable] > threshold))}
    model['half_life_minutes'] = (float(-5 * math.log(2) / math.log(model['phi']))
                                if 0 < model['phi'] < 1 else 0. if model['phi'] == 0 else None)
    return model


def cache_root():
    return Path(os.environ.get('GREENWINTER_VOLATILITY_CACHE_DIR', '/tmp/greenwinter-volatility-v1'))


def calendar_month_coordinate(seconds):
    dates = np.asarray(seconds, dtype=np.int64).astype('datetime64[s]').astype('datetime64[D]')
    months = dates.astype('datetime64[M]')
    first = months.astype('datetime64[D]').astype(np.int64)
    following = (months + np.timedelta64(1, 'M')).astype('datetime64[D]').astype(np.int64)
    return months.astype(np.int64) % 12 + (dates.astype(np.int64) - first) / (following - first)


def routine_volatility(seconds, innovations, spikes):
    """35-minute sample standard deviation, excluding classified spike samples."""
    result = np.full(len(innovations), np.nan)
    if len(innovations) < 7:
        return result
    windows = np.lib.stride_tricks.sliding_window_view(innovations, 7)
    spike_windows = np.lib.stride_tricks.sliding_window_view(spikes, 7)
    time_windows = np.lib.stride_tricks.sliding_window_view(seconds, 7)
    weights = np.isfinite(windows) & ~spike_windows
    counts = weights.sum(axis=1)
    means = np.where(weights, windows, 0.).sum(axis=1) / np.maximum(1, counts)
    variance = np.where(weights, (windows - means[:, None]) ** 2, 0.).sum(axis=1) / np.maximum(1, counts - 1)
    complete = np.isfinite(windows).all(axis=1) & (time_windows[:, -1] - time_windows[:, 0] == 6 * STEP) & (counts >= 4)
    result[6:] = np.where(complete, np.sqrt(np.maximum(0., variance)), np.nan)
    return result


def signal_reference_for_model(seconds, prices, model):
    """Historical distribution of the same rolling routine volatility on screen."""
    epoch = datetime(1970, 1, 1)
    start = int((datetime.fromisoformat(model['reference_start']) - epoch).total_seconds())
    end = int((datetime.fromisoformat(model['reference_end']) - epoch).total_seconds())
    selected = (seconds >= start) & (seconds <= end)
    t, y = seconds[selected], prices[selected]
    _, _, eps, sigma = filter_path(t, y, predict(t, model), model)
    volatility = routine_volatility(t, eps, np.abs(eps / sigma) > model['shock_threshold'])
    valid = np.isfinite(volatility)
    times, scores = t[valid], volatility[valid]
    phases = calendar_month_coordinate(times) / 12
    slots = ((times // DAY + 3) % 7) * 24 + (times % DAY) // 3600
    probabilities = np.linspace(0, 1, 101)
    grid = np.empty((12, 168, 101))
    pooled = np.quantile(scores, probabilities) if len(scores) else np.zeros(101)
    for month in range(12):
        distance = np.abs(phases - (month + .5) / 12)
        seasonal = np.minimum(distance, 1 - distance) <= 45 / 365.2425
        for slot in range(168):
            matched = (slots == slot) & seasonal
            if np.count_nonzero(matched) < 60:
                matched = slots == slot
            if np.count_nonzero(matched) < 60:
                matched = (slots % 24 == slot % 24) & seasonal
            if np.count_nonzero(matched) < 60:
                matched = slots % 24 == slot % 24
            grid[month, slot] = np.quantile(scores[matched], probabilities) if np.count_nonzero(matched) >= 60 else pooled
    return {'p5': float(np.quantile(scores, .05)) if len(scores) else None,
            'p95': float(np.quantile(scores, .95)) if len(scores) else None,
            'count': len(scores), 'window_minutes': 35, 'measure': 'routine_std_dollars',
            'start': model['reference_start'], 'end': model['reference_end'],
            'calendar_version': 2, 'calendar_grid': np.round(grid[:, :, [5, 25, 50, 75, 95]], 5).tolist(),
            'percentile_grid': np.round(grid, 5).tolist(),
            'season_window_days': 45, 'minimum_group_samples': 60}


def calendar_reference_values(model, times, grid_name='calendar_grid'):
    month_coordinate = calendar_month_coordinate(times) - .5
    month_index = np.floor(month_coordinate).astype(np.int64) % 12
    month_weight = (month_coordinate % 1)[:, None]
    hour_coordinate = ((times // DAY + 3) % 7) * 24 + (times % DAY) / 3600 - .5
    hour_index = np.floor(hour_coordinate).astype(np.int64) % 168
    hour_weight = (hour_coordinate % 1)[:, None]
    grid = np.asarray(model['signal_reference'][grid_name])
    low = grid[month_index, hour_index] * (1 - hour_weight) + grid[month_index, (hour_index + 1) % 168] * hour_weight
    high = grid[(month_index + 1) % 12, hour_index] * (1 - hour_weight) + grid[(month_index + 1) % 12, (hour_index + 1) % 168] * hour_weight
    return low * (1 - month_weight) + high * month_weight


def signal_reference_band(model, first, last):
    epoch = datetime(1970, 1, 1)
    year = datetime.fromisoformat(model['assessment_start']).year
    start = max(first, int((datetime(year, 1, 1) - epoch).total_seconds()))
    end = min(last, int((datetime(year + 1, 1, 1) - epoch).total_seconds()))
    if end <= start:
        return None
    times = np.r_[np.arange(start, end, STEP), end]
    quantiles = calendar_reference_values(model, times)
    return {'year': year, 'quantiles': [5, 25, 50, 75, 95],
            'points': np.column_stack([times - first, np.round(quantiles, 5)]).tolist()}


def model_for_year(seconds, prices, year, region=None, signal_reference=False):
    path = cache_root() / f'{region}-{year}-model.json' if region else None
    model = None
    if path:
        try:
            model = json.loads(path.read_text())
            if model.get('version') == VERSION:
                if not signal_reference or (model.get('signal_reference', {}).get('window_minutes') == 35 and
                                            model.get('signal_reference', {}).get('measure') == 'routine_std_dollars' and
                                            model.get('signal_reference', {}).get('calendar_version') == 2):
                    return model
            else:
                model = None
        except (OSError, ValueError):
            pass
    if model is None:
        model = fit_model(seconds, prices, year)
    if model and signal_reference and (model.get('signal_reference', {}).get('window_minutes') != 35 or
                                      model.get('signal_reference', {}).get('measure') != 'routine_std_dollars' or
                                      model.get('signal_reference', {}).get('calendar_version') != 2):
        model['signal_reference'] = signal_reference_for_model(seconds, prices, model)
    if model and path:
        model['version'] = VERSION
        try:
            path.parent.mkdir(mode=0o700, parents=True, exist_ok=True)
            with tempfile.NamedTemporaryFile(mode='w', dir=path.parent, delete=False) as handle:
                json.dump(model, handle, allow_nan=False)
                temporary = handle.name
            os.replace(temporary, path)
        except OSError:
            pass
    return model


def decompose_year(seconds, prices, year, model):
    start = int((datetime(year, 1, 1) - datetime(1970, 1, 1)).total_seconds())
    end = int((datetime(year + 1, 1, 1) - datetime(1970, 1, 1)).total_seconds())
    # A fixed 90-day warmup ensures identical values in every dashboard view.
    context = (seconds >= start - 90 * DAY) & (seconds < end)
    t, y = seconds[context], prices[context]
    season = predict(t, model)
    level, carry, eps, sigma, shock_carry = filter_path(t, y, season, model, components=True)
    z = np.abs(eps / sigma)
    tiers = np.where(z > model['severe_threshold'], 2, np.where(z > model['shock_threshold'], 1, 0))
    jumps = tiers > 0
    routine = np.where(jumps, 0., eps); shock = np.where(jumps, eps, 0.)
    routine[~np.isfinite(eps)] = np.nan; shock[~np.isfinite(eps)] = np.nan
    total, routine_rms, shock_rms = rolling_energy(t, eps, jumps)
    volatility_rank = rank(total, model['reference_vol']); surprise_rank = rank(z, model['reference_z'])
    weekday_hour = ((t // DAY + 3) % 7) * 24 + (t % DAY) // 3600
    risk = np.asarray(model['shock_risk'])[weekday_hour]
    values = np.column_stack([t, y, season, level, carry, routine, shock, eps, sigma,
                             total, routine_rms, shock_rms, volatility_rank, surprise_rank, tiers, risk, shock_carry])
    assessment_start = int((datetime.fromisoformat(model['assessment_start']) - datetime(1970, 1, 1)).total_seconds())
    values = values[t >= max(start, assessment_start)]
    return values


def cached_year(seconds, prices, year, model, region=None):
    """Reuse canonical full-year signals across views; append freshness by stamp."""
    end = int((datetime(year + 1, 1, 1) - datetime(1970, 1, 1)).total_seconds())
    available = seconds[seconds < end]
    stamp = int(available[-1]) if len(available) else 0
    path = cache_root() / f'{region}-{year}-signals.npz' if region else None
    if path:
        try:
            with np.load(path, allow_pickle=False) as saved:
                if int(saved['stamp']) >= stamp and str(saved['version']) == VERSION:
                    # The filter is causal: a later cached path already contains
                    # the identical prefix needed when navigating backwards.
                    values = saved['values']
                    return values[:np.searchsorted(values[:, 0], stamp, side='right')]
        except (OSError, ValueError, KeyError):
            pass
    values = decompose_year(seconds, prices, year, model)
    if path:
        try:
            path.parent.mkdir(mode=0o700, parents=True, exist_ok=True)
            with tempfile.NamedTemporaryFile(dir=path.parent, delete=False) as handle:
                np.savez_compressed(handle, values=values, stamp=stamp, version=VERSION)
                temporary = handle.name
            os.replace(temporary, path)
        except OSError:
            pass
    return values


def volatility_averages(values, period_start, period_end, view=None):
    """Average original five-minute percentile estimates before display thinning."""
    epoch = datetime(1970, 1, 1)
    first = int((period_start - epoch).total_seconds())
    span = int((period_end - period_start).total_seconds())
    if view is None:
        view = ('hourly' if span <= 3600 else 'weekly' if span <= 7 * DAY else
                'monthly' if span <= 31 * DAY else 'yearly' if span <= 366 * DAY else 'all_time')
    if view == 'all_time':
        edges = [0]
        boundary = period_start.replace(day=1, hour=0, minute=0, second=0, microsecond=0)
        while boundary < period_end:
            boundary = boundary.replace(year=boundary.year + 1, month=1) if boundary.month == 12 else boundary.replace(month=boundary.month + 1)
            edges.append(min(span, int((boundary - period_start).total_seconds())))
        label = 'Monthly averages'
    else:
        interval, label = {'hourly': (STEP, '5-minute averages'),
                           'weekly': (DAY, 'Daily averages'),
                           'monthly': (DAY, 'Daily averages'),
                           'yearly': (7 * DAY, '7-day averages')}[view]
        edges = list(range(0, span, interval)) + [span]
    valid = np.isfinite(values[:, 12]) & (values[:, 0] >= first) & (values[:, 0] < first + span)
    offsets = values[valid, 0] - first
    counts, _ = np.histogram(offsets, bins=edges)
    totals, _ = np.histogram(offsets, bins=edges, weights=values[valid, 12])
    return {'label': label, 'bars': [[int(edges[i]), int(edges[i + 1]), float(total / count), int(count)]
                                   for i, (total, count) in enumerate(zip(totals, counts)) if count]}


def volatility_line_series(seconds, prices, values, first, last):
    """Rolling 60-minute RMS, computed before display thinning; retain gaps/extrema."""
    context = (seconds >= first - 12 * STEP) & (seconds < last)
    t, y = seconds[context], prices[context]
    differences = np.full(len(y), np.nan)
    if len(y) > 1:
        differences[1:] = np.where(np.diff(t) == STEP, np.diff(y), np.nan)
    price_rms, _, _ = rolling_energy(t, differences, np.zeros(len(y), dtype=bool))
    displayed = t >= first
    t, price_rms = t[displayed], price_rms[displayed]
    unexpected = np.full(len(t), np.nan)
    if len(values):
        indexes = np.searchsorted(values[:, 0], t)
        matched = indexes < len(values)
        matched[matched] &= values[indexes[matched], 0] == t[matched]
        unexpected[matched] = values[indexes[matched], 9]
    rows = np.column_stack([t - first, price_rms, unexpected])
    chosen = set()
    width = max(1, int(math.ceil(len(rows) / 240)))
    for left in range(0, len(rows), width):
        block = rows[left:left + width]
        chosen.update((left, left + len(block) - 1))
        for column in (1, 2):
            valid = np.flatnonzero(np.isfinite(block[:, column]))
            if len(valid):
                chosen.update((left + int(valid[np.argmin(block[valid, column])]),
                               left + int(valid[np.argmax(block[valid, column])])))
    if len(rows) > 1:
        transitions = np.flatnonzero(np.any(np.diff(np.isfinite(rows[:, 1:]).astype(int), axis=0), axis=1))
        for index in transitions:
            chosen.update((int(index), int(index + 1)))
    points = [[int(rows[index, 0])] + [round(float(value), 4) if np.isfinite(value) else None
                                      for value in rows[index, 1:]] for index in sorted(chosen)]
    points.extend([[int(t[index] - first + STEP), None, None]
                   for index in np.flatnonzero(np.diff(t) != STEP)])
    points.sort(key=lambda row: row[0])
    return {'columns': ['time', 'price_rms_change', 'unexpected_rms'], 'points': points,
            'window_minutes': 60, 'interval_minutes': 5, 'units': '$/MWh'}


def analyze(rows, period_start, period_end, region=None, point=None, view=None, summary_only=False, include_series=False):
    """Every estimate uses original observations. Rendering keeps extrema only."""
    epoch = datetime(1970, 1, 1)
    valid_rows = ((int((t - epoch).total_seconds()), float(y)) for t, y in rows
                  if y is not None and math.isfinite(float(y)))
    raw = np.fromiter((v for row in valid_rows for v in row), dtype=float).reshape(-1, 2)
    if not len(raw):
        return {'available': False, 'reason': 'No valid five-minute observations.'}
    seconds = raw[:, 0].astype(np.int64)
    prices = raw[:, 1]
    order = np.argsort(seconds, kind='stable')
    seconds, prices = seconds[order], prices[order]
    if np.any(np.diff(seconds) == 0):
        raise ValueError('Five-minute observations must have unique timestamps')
    if np.any(seconds % STEP):
        raise ValueError('Observations must align to the original five-minute grid')
    first = int((period_start - epoch).total_seconds()); last = int((period_end - epoch).total_seconds())
    selected_prices = (seconds >= first) & (seconds < last)
    selected_seconds = seconds[selected_prices]
    observed_prices = prices[selected_prices]
    consecutive_pairs = np.diff(selected_seconds) == STEP
    price_changes = np.diff(observed_prices)[consecutive_pairs]
    measures = {
        'price_rms_change': float(np.sqrt(np.mean(price_changes ** 2))) if len(price_changes) else None,
        'unexpected_rms': None, 'routine_rms': None, 'spike_rms': None,
        'price_observations': int(len(observed_prices)), 'price_pairs': int(len(price_changes)),
        'surprise_observations': 0, 'model_coverage_percent': 0.,
        'units': '$/MWh', 'interval_minutes': 5,
    }
    # Allocate once: all-time requests must not duplicate nine years of arrays
    # when joining the separately cached yearly signals.
    selected_count = int(np.count_nonzero((seconds >= first) & (seconds < last)))
    values = np.empty((selected_count, len(COLUMNS)))
    weekly_routine = np.full(selected_count, np.nan)
    weekly_percentiles = np.full(selected_count, np.nan)
    component_columns = [1, 2, 3, 4, 16, 5, 6]
    weekly_changes = np.full((selected_count, len(component_columns)), np.nan) if view == 'weekly' else None
    previous_year_row = None
    position = 0; metadata = []; reference_bands = []
    first_year = (epoch + timedelta(seconds=point)).year if point is not None else period_start.year
    final_year = first_year if point is not None else min(period_end.year, (epoch + timedelta(seconds=int(seconds[-1]))).year)
    for year in range(first_year, final_year + 1):
        if datetime(year, 1, 1) >= period_end:
            continue
        model = model_for_year(seconds, prices, year, region, signal_reference=view == 'weekly')
        if model is None:
            continue
        if view == 'weekly':
            band = signal_reference_band(model, first, last)
            if band:
                reference_bands.append(band)
        path = cached_year(seconds, prices, year, model, region)
        selection = (path[:, 0] >= first) & (path[:, 0] < last)
        selected = path[selection]
        values[position:position + len(selected)] = selected
        if view == 'weekly':
            changes = np.full((len(path), len(component_columns)), np.nan)
            consecutive = np.diff(path[:, 0]) == STEP
            changes[1:] = np.where(consecutive[:, None], np.diff(path[:, component_columns], axis=0), np.nan)
            if len(path) and previous_year_row is not None and path[0, 0] - previous_year_row[0] == STEP:
                changes[0] = path[0, component_columns] - previous_year_row[component_columns]
            weekly_changes[position:position + len(selected)] = changes[selection]
            if len(path):
                previous_year_row = path[-1]
            routine = routine_volatility(path[:, 0], path[:, 7], path[:, 14] > 0)[selection]
            weekly_routine[position:position + len(selected)] = routine
            references = calendar_reference_values(model, selected[:, 0].astype(np.int64), 'percentile_grid')
            for index in np.flatnonzero(np.isfinite(routine)):
                reference = references[index]
                percentile = 50 * (np.searchsorted(reference, routine[index], side='left') +
                                   np.searchsorted(reference, routine[index], side='right')) / len(reference)
                weekly_percentiles[position + index] = percentile
        position += len(selected)
        metadata.append({key: model[key] for key in ('assessment_start', 'annual', 'phi', 'half_life_minutes', 'shock_threshold',
                        'severe_threshold', 'reference_start', 'reference_end', 'reference_count', 'reference_quantiles',
                        'overall_shock_risk_percent', 'jump_phi', 'jump_scale', 'validation') } |
                        {'year': year, 'signal_reference': {key: value for key, value in model.get('signal_reference', {}).items()
                                                          if key not in ('calendar_grid', 'percentile_grid')} or None})
    values = values[:position]
    weekly_routine = weekly_routine[:position]
    weekly_percentiles = weekly_percentiles[:position]
    if not len(values):
        result = {'available': False,
                'volatility_measures': measures,
                'reason': 'Requires at least eight weeks of earlier five-minute prices: four for fitting and four for temporal validation.'}
        if include_series:
            result['volatility_series'] = volatility_line_series(seconds, prices, values, first, last)
        return result
    usable_surprises = np.isfinite(values[:, 7])
    surprises = values[usable_surprises, 7]
    if len(surprises):
        measures.update({
            'unexpected_rms': float(np.sqrt(np.mean(surprises ** 2))),
            'routine_rms': float(np.sqrt(np.mean(values[usable_surprises, 5] ** 2))),
            'spike_rms': float(np.sqrt(np.mean(values[usable_surprises, 6] ** 2))),
            'surprise_observations': int(len(surprises)),
            'model_coverage_percent': 100 * len(surprises) / max(1, len(observed_prices)),
        })
    if summary_only:
        result = {'available': bool(len(surprises)), 'volatility_measures': measures}
        if include_series:
            result['volatility_series'] = volatility_line_series(seconds, prices, values, first, last)
        return result
    def serial(row):
        return [round(float(v - first) if i == 0 else float(v), 4) if np.isfinite(v) else None for i, v in enumerate(row)]
    if point is not None:
        index = int(np.searchsorted(values[:, 0], point))
        if index == len(values) or values[index, 0] != point:
            return {'available': False, 'reason': 'No valid observation at this five-minute timestamp.'}
        return {'available': True, 'columns': COLUMNS, 'point': serial(values[index]), 'models': metadata}
    average_values = values
    if view == 'weekly':
        average_values = values.copy()
        average_values[:, 12] = weekly_percentiles
    averages = volatility_averages(average_values, period_start, period_end, view)
    # Pixel envelopes preserve extremes across each component; gaps stay broken.
    selected = {0, len(values) - 1}
    width = max(1, int(math.ceil(len(values) / 200)))
    for left in range(0, len(values), width):
        block = values[left:left + width]
        selected.add(left); selected.add(left + len(block) - 1)
        for column in (1, 2, 3, 4, 7, 9, 12):
            indexes = np.flatnonzero(np.isfinite(block[:, column]))
            if len(indexes):
                selected.add(left + int(indexes[np.argmin(block[indexes, column])]))
                selected.add(left + int(indexes[np.argmax(block[indexes, column])]))
    gap_indexes = np.flatnonzero(np.diff(values[:, 0]) != STEP)
    signal_points = None
    price_components = None
    if view == 'weekly':
        # Retain every original observation for additive attribution and hover.
        # Changes use the preceding canonical observation, even outside the view.
        component_points = [serial(row) for row in values]
        component_changes = [[int(row[0] - first)] +
                             [round(float(value), 4) if np.isfinite(value) else None for value in change]
                             for row, change in zip(values, weekly_changes[:position])]
        component_points.extend([[int(values[index, 0] - first + STEP)] + [None] * (len(COLUMNS) - 1)
                                 for index in gap_indexes])
        component_changes.extend([[int(values[index, 0] - first + STEP)] + [None] * len(component_columns)
                                  for index in gap_indexes])
        component_points.sort(key=lambda row: row[0])
        component_changes.sort(key=lambda row: row[0])
        price_components = {'columns': COLUMNS, 'points': component_points,
                            'change_columns': ['time'] + [COLUMNS[column] for column in component_columns],
                            'changes': component_changes, 'units': '$/MWh', 'interval_minutes': 5}
        # Routine standard deviation and instantaneous spike magnitude share
        # dollar units; historical rarity still classifies spike observations.
        signal = np.abs(values[:, 7])
        spikes = values[:, 14] > 0
        routine_signal = weekly_routine
        spike_signal = np.where(spikes, signal, 0.)
        spike_signal[~np.isfinite(signal)] = np.nan
        signal_points = [[int(row[0] - first),
                          round(float(routine), 4) if np.isfinite(routine) else None,
                          round(float(spike), 4) if np.isfinite(spike) else None]
                         for row, routine, spike in zip(values, routine_signal, spike_signal)]
        signal_points.extend([[int(values[index, 0] - first + STEP), None, None] for index in gap_indexes])
        signal_points.sort(key=lambda row: row[0])
    for index in gap_indexes:
        selected.update((int(index), int(index + 1)))
    sampled = [serial(values[index]) for index in sorted(selected)]
    for index in gap_indexes:
        sampled.append([float(values[index, 0] - first + STEP)] + [None] * (len(COLUMNS) - 1))
    sampled.sort(key=lambda row: row[0])
    valid_vol = values[np.isfinite(values[:, 9])]
    eps = values[:, 7]; usable = np.isfinite(eps)
    energy = float(np.sum(eps[usable] ** 2))
    jump_energy = float(np.sum(values[usable, 6] ** 2))
    percentiles = values[np.isfinite(values[:, 12]), 12]
    histogram, edges = np.histogram(percentiles, bins=np.r_[np.arange(0, 100, 5), 99, 100])
    bands = np.where(np.isfinite(values[:, 12]), np.searchsorted([25, 75, 90, 99],
                     np.nan_to_num(values[:, 12]), side='right'), 5)
    boundaries = np.r_[0, np.flatnonzero((np.diff(bands) != 0) | (np.diff(values[:, 0]) != STEP)) + 1, len(values)]
    regimes = []
    for left, right in zip(boundaries[:-1], boundaries[1:]):
        regimes.append([int(values[left, 0] - first), int(values[right - 1, 0] - first + STEP), int(bands[left])])
    observed = np.count_nonzero((seconds >= first) & (seconds < last))
    return {'available': True, 'version': VERSION, 'columns': COLUMNS, 'points': sampled,
        'volatility_measures': measures,
        'volatility_averages': averages,
        'price_components': price_components,
        'volatility_signal': signal_points,
        'volatility_signal_measure': 'routine_std_dollars', 'volatility_window_minutes': 35,
        'routine_volatility_percentiles': [[int(row[0] - first), round(float(value), 4) if np.isfinite(value) else None]
                                          for row, value in zip(values, weekly_percentiles)] if view == 'weekly' else None,
        'signal_reference_bands': reference_bands,
        'latest': serial(values[-1]), 'models': metadata, 'observations': int(len(values)), 'regimes': regimes,
        'coverage_percent': 100 * len(values) / max(1, observed), 'window_minutes': 60,
        'summary': {'median_rms': float(np.median(valid_vol[:, 9])) if len(valid_vol) else None,
            'p95_rms': float(np.quantile(valid_vol[:, 9], .95)) if len(valid_vol) else None,
            'shock_energy_percent': 100 * jump_energy / energy if energy else 0.,
            'shock_count': int(np.count_nonzero(values[:, 14] > 0)),
            'extreme_count': int(np.count_nonzero(values[:, 14] == 2)),
            'high_volatility_percent': float(100 * np.mean(percentiles >= 90)) if len(percentiles) else None},
        'histogram': [[float(edges[i]), float(edges[i + 1]), float(100 * count / max(1, len(percentiles)))]
                      for i, count in enumerate(histogram)]}