# -*- coding: utf-8 -*-
#
# Isoliner - грид и изолинии (QGIS).
# © 2026 ООО «Информ++» (www.informpp.ru).
# SPDX-License-Identifier: GPL-2.0-or-later
#
# Это свободная программа: вы можете распространять её и/или изменять на
# условиях Стандартной общественной лицензии GNU (GNU GPL), опубликованной
# Фондом свободного ПО (FSF), - либо версии 2 Лицензии, либо (на ваше
# усмотрение) любой более поздней версии.
#
# Программа распространяется в надежде на полезность, но БЕЗ КАКИХ-ЛИБО
# ГАРАНТИЙ, в том числе без подразумеваемой гарантии ТОВАРНОГО СОСТОЯНИЯ или
# ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЁННОЙ ЦЕЛИ. Подробнее см. GNU GPL.
#
# Полный текст лицензии - в файле LICENSE (на английском, юридически значим).
"""
KB2D engine - pure Python/NumPy, no QGIS dependency (headless-testable).

Faithful port of the ArcGIS Isoliner add-in 2D kriging core
(GlobalOperation.vb: cova2_2D / Kriging2DInPoint), itself a port of the
GSLIB kb2d routine.

GSLIB conventions are preserved:
  it (model)  : 1 spherical, 2 exponential, 3 gaussian, 4 power
  ktype       : 0 = simple kriging (uses skmean), 1 = ordinary kriging
  ang         : azimuth, clockwise degrees from Y (north)
  anis        : minor-axis range / major-axis range  (<=1)
"""
import math
import numpy as np

DTOR = math.pi / 180.0
EPS = 1e-11

# UI model codes (0-based) -> GSLIB it (1-based)
MODEL_SPHERICAL = 0
MODEL_EXPONENTIAL = 1
MODEL_GAUSSIAN = 2
MODEL_POWER = 3


class Variogram:
    """Nugget + one or more nested structures (cova2_2D port)."""

    def __init__(self, nugget, structures):
        # structures: list of dict(it=1..4, cc, aa, ang, anis)
        self.c0 = float(nugget)
        self.it = [int(s["it"]) for s in structures]
        self.cc = [float(s["cc"]) for s in structures]
        self.aa = [max(float(s["aa"]), EPS) for s in structures]
        self.ang = [float(s["ang"]) for s in structures]
        self.anis = [max(float(s["anis"]), EPS) for s in structures]
        self.nst = len(structures)
        self.pmx = 9999.0
        self.rotmat = []
        self.maxcov = self.c0
        for i in range(self.nst):
            az = (90.0 - self.ang[i]) * DTOR
            ca, sa = math.cos(az), math.sin(az)
            self.rotmat.append((ca, sa, -sa, ca))
            self.maxcov += self.pmx if self.it[i] == 4 else self.cc[i]

    def cova2(self, dx, dy):
        """Covariance for a separation vector (dx, dy)."""
        if (dx * dx + dy * dy) < EPS:
            return self.maxcov
        cov = 0.0
        for i in range(self.nst):
            r0, r1, r2, r3 = self.rotmat[i]
            dx1 = dx * r0 + dy * r1
            dy1 = (dx * r2 + dy * r3) / self.anis[i]
            h = math.sqrt(max(dx1 * dx1 + dy1 * dy1, 0.0))
            t = self.it[i]
            if t == 1:                       # spherical
                hr = h / self.aa[i]
                if hr < 1.0:
                    cov += self.cc[i] * (1.0 - hr * (1.5 - 0.5 * hr * hr))
            elif t == 2:                     # exponential
                cov += self.cc[i] * math.exp(-h / self.aa[i])
            elif t == 3:                     # gaussian
                cov += self.cc[i] * math.exp(-(h * h) / (self.aa[i] * self.aa[i]))
            else:                            # power
                cov += self.pmx - self.cc[i] * (h ** self.aa[i])
        return cov


def clip_outliers(values, vmin=None, vmax=None, pct=0.0, cap=False):
    """Отсев/срезка ураганных проб. Возвращает (values_out, keep_mask, lo, hi).

    Границы: явные vmin/vmax имеют приоритет; иначе - перцентильные [pct;
    100-pct] (если pct>0); иначе границ нет. Режимы:
      cap=False - отсеять: keep_mask=False у проб вне [lo; hi];
      cap=True  - срезать: значения вне [lo; hi] прижимаются к границе,
                  keep_mask остаётся True для всех.
    Универсальное применение: xs, ys, v = xs[keep], ys[keep], out[keep].
    """
    v = np.asarray(values, float)
    if vmin is None and vmax is None and (pct is None or pct <= 0):
        return v, np.ones(len(v), bool), float("-inf"), float("inf")
    lo, hi = float("-inf"), float("inf")
    if pct and pct > 0 and len(v):
        p = min(max(float(pct), 0.0), 49.0)
        lo = float(np.percentile(v, p))
        hi = float(np.percentile(v, 100.0 - p))
    if vmin is not None:
        lo = float(vmin)
    if vmax is not None:
        hi = float(vmax)
    if lo > hi:
        lo, hi = hi, lo
    if cap:
        return np.clip(v, lo, hi), np.ones(len(v), bool), lo, hi
    return v, (v >= lo) & (v <= hi), lo, hi


def _solve_point(xloc, yloc, xd, yd, vrd, vg, ktype, skmean,
                 ndmin, ndmax, rad2, nodata):
    """Точечный кригинг в одной точке. Возвращает (оценка, дисперсия).

    Дисперсия - это дисперсия ошибки кригинга:
        SK:  σ² = C(0) − Σ λ_i C(x_i, x0)
        OK:  σ² = C(0) − Σ λ_i C(x_i, x0) − μ   (μ - множитель Лагранжа)
    При совпадении узла с пробой σ² = 0. Если оценка не получена (вырожденная
    система) - дисперсия принимается равной априорной (силл). Когда соседей
    меньше ndmin, и оценка, и дисперсия = nodata."""
    dx = xd - xloc
    dy = yd - yloc
    h2 = dx * dx + dy * dy
    order = np.argsort(h2)
    sel = order[h2[order] <= rad2][:ndmax]
    na = len(sel)
    if na < ndmin:
        return nodata, nodata

    xa = xd[sel]
    ya = yd[sel]
    vra = vrd[sel]
    h2sel = h2[sel]
    cbb = vg.maxcov                          # point support: C(0)

    # совпадение узла с пробой -> возвращаем значение пробы, дисперсия 0
    if h2sel[0] < EPS:
        return float(vra[0]), 0.0

    if na == 1:
        cb = vg.cova2(xa[0] - xloc, ya[0] - yloc)
        if ktype == 0:                       # simple
            s = cb / cbb
            return s * vra[0] + (1.0 - s) * skmean, max(cbb - s * cb, 0.0)
        # ordinary, единственная проба: λ=1, μ=cb−C(0) -> σ²=2(C(0)−cb)
        return float(vra[0]), max(2.0 * (cbb - cb), 0.0)

    # допустимый разброс оценки (защита от «разлёта» весов)
    vmin = float(vra.min()); vmax = float(vra.max())
    span = (vmax - vmin) or (abs(vmax) + 1.0)
    lo, hi = vmin - 3.0 * span, vmax + 3.0 * span
    jitter = cbb * 1e-9                       # микро-регуляризация диагонали

    neq = na + ktype                         # +1 row for OK unbiasedness
    A = np.empty((neq, neq))
    r = np.empty(neq)
    for i in range(na):
        A[i, i] = cbb + jitter
        for j in range(i + 1, na):
            c = vg.cova2(xa[i] - xa[j], ya[i] - ya[j])
            A[i, j] = c
            A[j, i] = c
        r[i] = vg.cova2(xa[i] - xloc, ya[i] - yloc)
    if ktype == 1:                           # ordinary kriging
        A[na, :na] = vg.maxcov
        A[:na, na] = vg.maxcov
        A[na, na] = 0.0
        r[na] = vg.maxcov

    est = None
    var = None
    try:
        s = np.linalg.solve(A, r)
        if np.all(np.isfinite(s)):
            w = s[:na]
            e = float(np.dot(w, vra))
            if ktype == 0:                   # simple kriging mean term
                e += (1.0 - w.sum()) * skmean
            v = cbb - float(np.dot(w, r[:na]))
            if ktype == 1:                   # минус множитель Лагранжа
                v -= vg.maxcov * float(s[na])   # (масштабирован на C(0))
            if math.isfinite(e) and lo <= e <= hi:
                est = e                      # принимаем только разумную оценку
                var = max(v, 0.0)
    except np.linalg.LinAlgError:
        pass

    if est is None:                          # запас: обратные расстояния соседей
        wts = 1.0 / np.maximum(h2sel, EPS)
        est = float(np.dot(wts, vra) / wts.sum())
        var = cbb                            # априорная (максимальная) дисперсия
    return est, var


def krige_point(xloc, yloc, xd, yd, vrd, vg, ktype, skmean,
                ndmin, ndmax, rad2, nodata, return_var=False):
    """Точечный кригинг (обёртка). По умолчанию возвращает только оценку;
    при return_var=True - кортеж (оценка, дисперсия)."""
    est, var = _solve_point(xloc, yloc, xd, yd, vrd, vg, ktype, skmean,
                            ndmin, ndmax, rad2, nodata)
    return (est, var) if return_var else est


def cross_validate(xd, yd, vrd, vg, ktype, skmean, ndmin, ndmax,
                   rad2, nodata, progress=None):
    """Скользящий контроль (leave-one-out). Для каждой точки оценка строится
    по всем ОСТАЛЬНЫМ точкам, затем сравнивается с фактическим значением.
    Возвращает (est, var) той же длины, что и входные данные (nodata там, где
    оценка не получена). Это основа подбора вариограммы по ошибке."""
    n = len(xd)
    est = np.full(n, nodata, float)
    var = np.full(n, nodata, float)
    keep = np.ones(n, bool)
    for i in range(n):
        keep[i] = False                      # исключаем саму точку
        e, v = _solve_point(xd[i], yd[i], xd[keep], yd[keep], vrd[keep],
                            vg, ktype, skmean, ndmin, ndmax, rad2, nodata)
        keep[i] = True
        est[i] = e
        var[i] = v
        if progress is not None and (i & 255) == 0:
            progress(i, n)
    return est, var


def build_grid(xd, yd, vrd, vg, ktype, skmean, ndmin, ndmax,
               rad2, nodata, xmn, ymn, cell, nx, ny, progress=None,
               with_variance=False):
    """Sweep the grid, GSLIB order (north row first). Returns float32 (ny,nx).

    При with_variance=True возвращает кортеж (оценка, стд.ошибка), где второй
    грид - стандартная ошибка кригинга = sqrt(дисперсия) (nodata там же, где
    nodata у оценки)."""
    grid = np.full((ny, nx), nodata, dtype=np.float32)
    sgrid = np.full((ny, nx), nodata, dtype=np.float32) if with_variance else None
    total = nx * ny
    done = 0
    for row in range(ny):                    # row 0 = north
        iy = ny - row                        # 1..ny counted from south
        yloc = ymn + (iy - 1) * cell
        for ix in range(nx):
            xloc = xmn + ix * cell
            e, v = _solve_point(
                xloc, yloc, xd, yd, vrd, vg, ktype, skmean,
                ndmin, ndmax, rad2, nodata)
            grid[row, ix] = e
            if with_variance and e != nodata:
                sgrid[row, ix] = math.sqrt(max(v, 0.0))
            done += 1
        if progress is not None:
            progress(done, total)
    return (grid, sgrid) if with_variance else grid


# ===========================================================================
#  Экспериментальная вариограмма и подбор модели (чистый NumPy, без QGIS).
#  Параметризация моделей совпадает с cova2 выше: радиус a здесь - это тот же
#  параметр aa структуры, поэтому подобранный a подставляется в кригинг 1:1.
# ===========================================================================
def variogram_shape(model, h, a):
    """Форма вариограммы γ/C (от 0 до 1) для одной структуры, без наггета.

    model: 0 сферическая, 1 экспоненциальная, 2 гауссова (как MODEL_*).
    Совпадает с cova2: для сферической при h>=a даёт 1; для экспоненциальной
    1-exp(-h/a); для гауссовой 1-exp(-(h/a)^2). Параметр a = aa структуры.
    """
    h = np.asarray(h, float)
    a = max(float(a), EPS)
    hr = h / a
    if model == MODEL_SPHERICAL:
        s = 1.5 * hr - 0.5 * hr ** 3
        return np.where(hr < 1.0, s, 1.0)
    if model == MODEL_EXPONENTIAL:
        return 1.0 - np.exp(-hr)
    if model == MODEL_GAUSSIAN:
        return 1.0 - np.exp(-(hr * hr))
    raise ValueError("fit поддерживает только модели 0/1/2 (sph/exp/gauss)")


def experimental_variogram(xs, ys, vs, n_lags=15, maxlag=None, robust=False,
                           max_pairs=6_000_000, cloud_max=20000, seed=0):
    """Омнинаправленная экспериментальная полувариограмма (оценка Матерона).

    Делит расстояния между всеми парами на n_lags интервалов до maxlag и
    усредняет полудисперсию 0.5*(z_i - z_j)^2 по каждому интервалу.
      maxlag=None/<=0 -> половина диагонали охвата (классический ориентир);
      robust=True     -> устойчивая оценка Кресси-Хокинса (гасит выбросы);
      max_pairs       -> при превышении точки прореживаются (пар ~ n^2/2),
                         чтобы расчёт оставался быстрым;
      cloud_max       -> сколько пар вернуть для облака рассеяния (предпросмотр).
    Возвращает dict: lag, gamma, npairs (по непустым лагам), maxlag, width,
    n_used, subsampled, cloud_h, cloud_g (полудисперсия пары).
    """
    xs = np.asarray(xs, float)
    ys = np.asarray(ys, float)
    vs = np.asarray(vs, float)
    n = len(xs)
    rng = np.random.default_rng(seed)
    subsampled = False
    if n > 2 and n * (n - 1) // 2 > max_pairs:
        m = int((1.0 + math.sqrt(1.0 + 8.0 * max_pairs)) / 2.0)
        m = max(2, min(n, m))
        idx = rng.choice(n, m, replace=False)
        xs, ys, vs = xs[idx], ys[idx], vs[idx]
        n = m
        subsampled = True

    if maxlag is None or maxlag <= 0:
        dx = float(xs.max() - xs.min())
        dy = float(ys.max() - ys.min())
        maxlag = 0.5 * math.hypot(dx, dy)
    if maxlag <= 0:
        maxlag = 1.0
    n_lags = max(int(n_lags), 1)
    width = maxlag / n_lags

    cnt = np.zeros(n_lags)
    s_h = np.zeros(n_lags)
    s_d2 = np.zeros(n_lags)
    s_sqrt = np.zeros(n_lags)
    cloud_h, cloud_g = [], []
    budget = int(cloud_max)

    for i in range(n - 1):
        hx = xs[i + 1:] - xs[i]
        hy = ys[i + 1:] - ys[i]
        h = np.sqrt(hx * hx + hy * hy)
        d = vs[i + 1:] - vs[i]
        sel = (h > 0.0) & (h <= maxlag)
        if not sel.any():
            continue
        h = h[sel]
        d = d[sel]
        b = np.minimum((h / width).astype(np.intp), n_lags - 1)
        np.add.at(cnt, b, 1.0)
        np.add.at(s_h, b, h)
        np.add.at(s_d2, b, d * d)
        if robust:
            np.add.at(s_sqrt, b, np.sqrt(np.abs(d)))
        if budget > 0:
            k = len(h)
            take = min(k, budget)
            if take < k:
                pick = rng.choice(k, take, replace=False)
                cloud_h.append(h[pick]); cloud_g.append(0.5 * d[pick] * d[pick])
            else:
                cloud_h.append(h); cloud_g.append(0.5 * d * d)
            budget -= take

    safe = np.maximum(cnt, 1.0)
    lag = np.where(cnt > 0, s_h / safe, (np.arange(n_lags) + 0.5) * width)
    if robust:
        mean_sqrt = np.where(cnt > 0, s_sqrt / safe, 0.0)
        denom = 0.457 + 0.494 / safe + 0.045 / (safe * safe)
        gamma = 0.5 * (mean_sqrt ** 4) / denom
    else:
        gamma = 0.5 * np.where(cnt > 0, s_d2 / safe, np.nan)

    valid = cnt > 0
    return {
        "lag": lag[valid],
        "gamma": gamma[valid],
        "npairs": cnt[valid].astype(int),
        "maxlag": float(maxlag),
        "width": float(width),
        "n_used": int(n),
        "subsampled": bool(subsampled),
        "cloud_h": (np.concatenate(cloud_h) if cloud_h else np.array([])),
        "cloud_g": (np.concatenate(cloud_g) if cloud_g else np.array([])),
    }


def _wls_c0_c(s, y, w):
    """Взвешенный МНК для γ = c0 + C*s при фикс. форме s, c0>=0, C>=0.
    Возвращает (c0, C, sse)."""
    sw = float(np.sum(w))
    ss = float(np.sum(w * s))
    sss = float(np.sum(w * s * s))
    sy = float(np.sum(w * y))
    ssy = float(np.sum(w * s * y))
    det = sw * sss - ss * ss
    c0 = c = 0.0
    if abs(det) > 1e-30:
        c0 = (sy * sss - ss * ssy) / det
        c = (sw * ssy - ss * sy) / det
    if c0 < 0.0 or c < 0.0 or abs(det) <= 1e-30:
        # граница допустимой области: перебираем варианты с занулением
        cands = []
        # c0=0 -> C по МНК
        cc = ssy / sss if sss > 1e-30 else 0.0
        cands.append((0.0, max(cc, 0.0)))
        # C=0 -> c0 = взвешенное среднее
        cands.append((sy / sw if sw > 0 else 0.0, 0.0))
        # оба >=0 (если вышло)
        if c0 >= 0.0 and c >= 0.0:
            cands.append((c0, c))
        best = None
        for a0, a1 in cands:
            r = y - a0 - a1 * s
            e = float(np.sum(w * r * r))
            if best is None or e < best[2]:
                best = (a0, a1, e)
        return best
    r = y - c0 - c * s
    return c0, c, float(np.sum(w * r * r))


def fit_variogram(lag, gamma, npairs, model="auto", n_a=80, sill_cap=None):
    """Подбор модели вариограммы по экспериментальным точкам (рекомендация).

    γ(h) = c0 + C*shape(h; a). При фиксированном a задача линейна по (c0, C),
    поэтому a сканируется по сетке (лог-шаг) с локальным уточнением, а (c0, C)
    на каждом шаге решаются взвешенным МНК (веса = число пар в лаге; так точнее
    оценённые лаги весомее). Без scipy/skgstat.

    Плато модели c0+C ограничивается сверху (sill_cap), чтобы модель не
    «убегала» в несуществующий силл за счёт длинного радиуса и экстраполяции
    выше наблюдённой вариограммы. По умолчанию cap = 1.15*max(γ). Благодаря
    этому auto не выбирает гауссову там, где она выигрывает лишь экстраполяцией.

    model: 'auto' (выбрать sph/exp/gauss по лучшему R^2) либо 0/1/2.
    Возвращает dict: model, nugget, sill, range, r2, npts. None при нехватке точек.
    """
    lag = np.asarray(lag, float)
    gamma = np.asarray(gamma, float)
    w = np.asarray(npairs, float)
    ok = np.isfinite(lag) & np.isfinite(gamma) & (lag > 0)
    lag, gamma, w = lag[ok], gamma[ok], np.maximum(w[ok], 1.0)
    if len(lag) < 3:
        return None
    if sill_cap is None:
        sill_cap = 1.15 * float(np.max(gamma))

    L = float(lag.max())
    a_lo = max(0.05 * L, float(lag.min()) * 0.5, EPS)
    a_hi = 1.5 * L
    grid = np.geomspace(a_lo, a_hi, max(int(n_a), 8))
    models = [MODEL_SPHERICAL, MODEL_EXPONENTIAL, MODEL_GAUSSIAN] \
        if model == "auto" else [int(model)]

    sw = float(np.sum(w))
    ybar = float(np.sum(w * gamma) / sw)
    sst = float(np.sum(w * (gamma - ybar) ** 2)) or 1.0

    def _eval(m, a):
        s = variogram_shape(m, lag, a)
        c0, c, sse = _wls_c0_c(s, gamma, w)
        total = c0 + c
        if sill_cap and total > sill_cap and total > 0:
            f = sill_cap / total                 # прижать плато к cap
            c0, c = c0 * f, c * f
            r = gamma - c0 - c * s
            sse = float(np.sum(w * r * r))
        return c0, c, sse

    best = None
    for m in models:
        coarse = [(a, _eval(m, a)) for a in grid]
        a_star, (c0, c, sse) = min(coarse, key=lambda t: t[1][2])
        j = [i for i, (a, _) in enumerate(coarse) if a == a_star][0]
        alo = coarse[max(j - 1, 0)][0]
        ahi = coarse[min(j + 1, len(coarse) - 1)][0]
        for _ in range(40):
            a1 = alo * (ahi / alo) ** (1.0 / 3.0)
            a2 = alo * (ahi / alo) ** (2.0 / 3.0)
            if _eval(m, a1)[2] < _eval(m, a2)[2]:
                ahi = a2
            else:
                alo = a1
            if ahi / alo < 1.0001:
                break
        a_ref = math.sqrt(alo * ahi)
        c0, c, sse = _eval(m, a_ref)
        r2 = 1.0 - sse / sst
        cand = {"model": m, "nugget": float(max(c0, 0.0)),
                "sill": float(max(c, 0.0)), "range": float(a_ref),
                "r2": float(r2), "npts": int(len(lag))}
        if best is None or cand["r2"] > best["r2"]:
            best = cand
    return best


def model_curve(vg, hmax, ndir=None, npts=120):
    """Кривая γ(h) заданной модели для наложения на экспериментальную.

    Возвращает (h, gamma_major[, gamma_minor]). По умолчанию вдоль главной оси
    первой структуры (азимут структуры 1); если анизотропия != 1, второй
    кривой добавляется малая ось - так видно обе ветви. Изотропная модель даёт
    одну кривую (ветви совпадают, малая не возвращается).
    """
    h = np.linspace(0.0, float(hmax), int(npts))
    az = vg.ang[0] if vg.nst else 0.0
    ar = az * DTOR
    ux, uy = math.sin(ar), math.cos(ar)          # главная ось (азимут от севера)
    g_major = np.array([vg.maxcov - vg.cova2(hi * ux, hi * uy) for hi in h])
    anis_min = min(vg.anis) if vg.nst else 1.0
    if anis_min < 0.999:
        vx, vy = math.cos(ar), -math.sin(ar)     # перпендикуляр (малая ось)
        g_minor = np.array([vg.maxcov - vg.cova2(hi * vx, hi * vy) for hi in h])
        return h, g_major, g_minor
    return h, g_major
