"""Подбор конкурирующих моделей отклонения и эквивалентного колебательного контура."""  # Численный модуль использует стандартную библиотеку и не зависит от Qt.
from bisect import bisect_left  # Поиск соседних отсчётов сохраняет реальные, в том числе неравномерные временные сетки.
import math  # Экспоненты, тригонометрия и логарифмы задают динамические модели и критерий BIC.
from statistics import median  # Робастные начальные оценки уменьшают влияние отдельных выбросов.

MAX_FIT_POINTS = 2048  # Ограничение относится только к подбору моделей; исходные метрики и поиск событий сохраняют все точки.


def checkpoint(cancel):  # Подбор нелинейных параметров обязан прерываться при смене записи или закрытии окна.
    if cancel():  # Callback передаётся из существующего рабочего потока сравнения.
        raise InterruptedError('Анализ моделей отменён')  # Не публикуем незавершённый подбор как достоверный результат.


def bounded_minimum(operation, left, right, iterations=24):  # Ограниченный золотой поиск уточняет одну нелинейную координату без внешнего оптимизатора.
    ratio = (math.sqrt(5.) - 1.) / 2.  # Золотое сечение переиспользует одну оценку функции на каждом шаге.
    a, b = right - ratio * (right - left), left + ratio * (right - left)  # Две внутренние точки определяют следующую область поиска.
    fa, fb = operation(a), operation(b)  # Операция возвращает модель с обучающей ошибкой или отсутствие допустимого решения.
    error = lambda model: model['training_sse'] if model else math.inf  # Некорректный кандидат не может выигрывать ограниченную оптимизацию.
    for _ in range(iterations):  # Фиксированное число итераций ограничивает время и задаёт воспроизводимую точность.
        if error(fa) < error(fb):  # Минимум остаётся в левой части текущей области.
            right, b, fb = b, a, fa  # Правая внутренняя точка переиспользует уже вычисленное значение.
            a = right - ratio * (right - left)  # Новая точка лежит внутри уменьшенного интервала.
            fa = operation(a)  # Пересчитываем только одну модель за итерацию.
        else:  # Лучшее направление лежит справа или обе оценки одинаковы.
            left, a, fa = a, b, fb  # Левая точка также переиспользуется без повторного подбора.
            b = left + ratio * (right - left)  # Вычисляем новую внутреннюю правую точку.
            fb = operation(b)  # Контрольные наблюдения по-прежнему не участвуют в целевой функции.
    return fa if error(fa) < error(fb) else fb  # Возвращаем лучший допустимый набор параметров найденного интервала.


def interpolate(times, values, stamp):  # Интерполяция для моделей не продолжает кривую за границы записи.
    index = bisect_left(times, stamp)  # Находим первую опору не левее запрошенного времени.
    if index < len(times) and times[index] == stamp:  # Точное попадание не зависит от соседнего пропуска.
        return values[index] if math.isfinite(values[index]) else None  # Неконечная опора остаётся отсутствующим наблюдением.
    if index == 0 or index == len(times):  # За границами данных экстраполяция запрещена.
        return None  # Модель с отсутствующей опорой не получает искусственного крайнего значения.
    left, right = values[index - 1], values[index]  # Две соседние амплитуды определяют отрезок интерполяции.
    if not math.isfinite(left) or not math.isfinite(right):  # Через явный разрыв NaN модель не строится.
        return None  # Невалидный участок будет исключён до сравнения моделей.
    return left + (right - left) * (stamp - times[index - 1]) / (times[index] - times[index - 1])  # Вес вычисляется по фактическому времени, а не номеру точки.


def least_squares(columns, values):  # Небольшой QR-решатель избегает неустойчивых нормальных уравнений при разном масштабе параметров.
    count, size = len(values), len(columns)  # Число наблюдений и неизвестных определяет размеры ортогональной системы.
    if count <= size or not size:  # Для оценки ошибки нужны наблюдения сверх числа параметров.
        return None  # Недоопределённую модель нельзя использовать для аппаратной гипотезы.
    q, upper = [], [[0.] * size for _ in range(size)]  # Храним только малую треугольную матрицу и ортогональные столбцы.
    scales = [math.sqrt(sum(value * value for value in column)) for column in columns]  # Нормируем столбцы перед проверкой независимости параметров.
    for index, column in enumerate(columns):  # Модифицированный Грамм–Шмидт выполняется отдельно для каждого неизвестного.
        if scales[index] <= 1e-30:  # Нулевой столбец означает, что параметр не наблюдается в выбранных данных.
            return None  # Не выдаём произвольный коэффициент для ненаблюдаемого воздействия.
        vector = [value / scales[index] for value in column]  # Безразмерная нормировка снижает влияние единиц В, А и мкс на численную устойчивость.
        for _ in range(2):  # Повторная ортогонализация уменьшает потерю точности у почти зависимых столбцов.
            for previous, basis in enumerate(q):  # Удаляем уже описанные направления изменения сигнала.
                projection = sum(a * b for a, b in zip(basis, vector))  # Скалярное произведение определяет долю предыдущего базиса.
                upper[previous][index] += projection  # Обе ортогонализации относятся к одной треугольной системе.
                vector = [value - projection * base for value, base in zip(vector, basis)]  # Остаток содержит только новое независимое направление.
        norm = math.sqrt(sum(value * value for value in vector))  # Длина остатка показывает идентифицируемость следующего параметра.
        if norm < 1e-9:  # Относительная независимость меньше 10^-9 считается численно вырожденной.
            return None  # Например, постоянный ток не позволяет раздельно определить U0 и сопротивление.
        upper[index][index] = norm  # Диагональ задаёт масштаб нового ортогонального направления.
        q.append([value / norm for value in vector])  # Единичный столбец сохраняет устойчивую проекцию правой части.
    target = [sum(a * b for a, b in zip(basis, values)) for basis in q]  # Проецируем измерение на полученный ортонормированный базис.
    coefficients = [0.] * size  # Обратный ход начинается с последнего неизвестного.
    for index in reversed(range(size)):  # Верхнетреугольная система решается без дополнительной факторизации.
        coefficients[index] = (target[index] - sum(upper[index][j] * coefficients[j] for j in range(index + 1, size))) / upper[index][index]  # Вычитаем уже найденные параметры справа от диагонали.
    return [value / scale for value, scale in zip(coefficients, scales)]  # Возвращаем параметры в исходных физических единицах столбцов.


def evaluate(columns, coefficients):  # Одна формула вычисляет прогноз как для обучающих, так и для контрольных точек.
    return [sum(column[index] * value for column, value in zip(columns, coefficients)) for index in range(len(columns[0]))]  # Сумма базисных функций не обращается к измеренной амплитуде.


def fit_model(code, formula, columns, values, names, extra_parameters=0, params=None):  # Каждый кандидат оценивается на одних и тех же наблюдениях.
    training = [index for index in range(len(values)) if index % 3 != 1]  # Каждая третья точка оставляется для проверки интерполяционного прогноза.
    coefficients = least_squares([[column[i] for i in training] for column in columns], [values[i] for i in training])  # Контрольные амплитуды не участвуют в оценке линейных коэффициентов.
    if coefficients is None:  # Вырожденная модель не должна выигрывать из-за произвольных параметров.
        return None  # Вызывающая сторона честно сообщает отсутствие идентифицируемой оценки.
    prediction = evaluate(columns, coefficients)  # Один набор параметров применяется также к отложенным точкам.
    parameters = dict(zip(names, coefficients))  # Имена отражают физический смысл оценённых коэффициентов.
    parameters.update(params or {})  # Нелинейные параметры включаются в тот же машиночитаемый отчёт.
    return score_model(code, formula, values, prediction, len(columns) + extra_parameters, parameters)  # BIC учитывает и линейные, и нелинейные степени свободы.


def score_model(code, formula, values, prediction, parameters, params):  # Показатели качества имеют одинаковое определение для всех конкурентов.
    residuals = [value - fitted for value, fitted in zip(values, prediction)]  # Остаток всегда сохраняет знак «измерение минус модель».
    square = sum(value * value for value in residuals)  # SSE измеряется в квадрате единиц анализируемого канала.
    validation = [value for index, value in enumerate(residuals) if index % 3 == 1]  # Считаем ошибку только по точкам, не использованным для линейного подбора.
    training_square = sum(value * value for index, value in enumerate(residuals) if index % 3 != 1)  # Нелинейный поиск оптимизирует только обучающую часть.
    training_count = len(values) - len(validation)  # BIC вычисляется на тех же наблюдениях, на которых оценивались параметры модели.
    mean = sum(values) / len(values)  # Среднее определяет постоянную модель для вычисления R².
    variance = sum((value - mean) ** 2 for value in values)  # Полный разброс нужен только для описания качества, а не для постановки диагноза.
    residual_mean = sum(residuals) / len(residuals)  # Среднее остатка вычисляем один раз, сохраняя линейную сложность оценки кандидата.
    centered = [value - residual_mean for value in residuals]  # Автокорреляция оценивается после удаления среднего остатка.
    denominator = sum(value * value for value in centered)  # Нулевой остаток не имеет определённой автокорреляции.
    floor = max(1e-24, max(abs(value) for value in values) ** 2 * 1e-24)  # Численный пол 10^-12 амплитудного масштаба не является допуском исправности.
    return dict(code=code, formula=formula, parameters=params, parameter_count=parameters, count=len(values), sse=square, training_sse=training_square,  # Сохраняем параметры и число степеней свободы для независимого воспроизведения.
                rmse=math.sqrt(square / len(values)), validation_rmse=math.sqrt(sum(value * value for value in validation) / max(1, len(validation))),  # Контрольная ошибка отделена от ошибки по всей выборке.
                bic=training_count * math.log(max(training_square / training_count, floor)) + parameters * math.log(training_count), bic_count=training_count,  # Для гауссовой ошибки сравниваем BIC обучения, сохраняя контрольные точки вне оценки параметров и критерия.
                r_squared=1. - square / variance if variance > floor else None, residual_lag1=sum(a * b for a, b in zip(centered, centered[1:])) / denominator if denominator > floor else None)  # Корреляция остатка показывает неописанную динамику и ограничение гауссового критерия.


def _samples(result, maximum=MAX_FIT_POINTS, cancel=lambda: False):  # Стратифицированная выборка сохраняет экстремумы ошибки каждого временного блока и поддерживает отмену.
    count = len(result.difference.times)  # Метрики и старые признаки уже рассчитаны по всем исходным отсчётам.
    width = max(1, math.ceil(count / max(1, maximum // 4)))  # В каждом блоке оставляем начало, конец и два экстремума остатка.
    selected = set()  # Совпадающие индексы не должны менять вес отдельного наблюдения.
    for start in range(0, count, width):  # Число выбранных точек ограничено независимо от размера исходного CSV.
        checkpoint(cancel)  # Подготовка многомиллионной записи прерывается между небольшими временными блоками.
        stop = min(count, start + width)  # Последний неполный блок также участвует в подборе.
        valid = [i for i in range(start, stop) if math.isfinite(result.expected[i]) and math.isfinite(result.difference.values[i])]  # Пропуски остаются исключёнными из подбора.
        if valid:  # Пустой блок не создаёт фиктивного измерения.
            selected.update((valid[0], valid[-1], min(valid, key=lambda i: result.difference.values[i]), max(valid, key=lambda i: result.difference.values[i])))  # Узкие импульсы сохраняются даже при длинной записи.
    return [(result.difference.times[i] * 1000., result.expected[i], result.expected[i] + result.difference.values[i]) for i in sorted(selected)]  # Модели используют микросекунды и исходные единицы амплитуды.


def fit_reference_models(result, events, cancel=lambda: False):  # Конкурируют нормальный эталон, калибровка, дрейф, задержка, растяжение и разрывы.
    samples = _samples(result, cancel=cancel)  # Подбор ограничен по объёму и не меняет полную статистику сравнения.
    if len(samples) < 24:  # Трёхчастное разбиение требует достаточно обучающих и контрольных точек.
        return dict(models=[], best=None, sampled_points=len(samples), sampling='стратифицированные экстремумы остатка')  # Не называем короткую запись пригодной для идентификации модели.
    duration = samples[-1][0] - samples[0][0]  # Длительность выражена в микросекундах.
    margin = .1 * duration if events else 0.  # Общий интервал опор исключает экстраполяцию при поиске задержки до 10% длительности.
    reference_times = result.reference.times  # Повторно используем полную миллисекундную сетку без многомиллионной дополнительной копии.
    samples = [point for point in samples if reference_times[0] * 1000. + margin <= point[0] <= reference_times[-1] * 1000. - margin]  # Все глобальные модели сравниваются на строго одинаковой области поддержки в мкс.
    if len(samples) < 24:  # Слишком короткое перекрытие не позволяет честно сравнивать кандидатов.
        return dict(models=[], best=None, sampled_points=len(samples), sampling='нет общего интервала моделей')  # Не подменяем отсутствующие точки экстраполяцией.
    times, reference, measured = map(list, zip(*samples))  # Параллельные массивы имеют одинаковую длину и порядок.
    ones = [1.] * len(times)  # Свободный член описывает смещение нуля в единицах канала.
    models = [score_model('reference', 'y(t)=r(t)', measured, reference, 0, {})]  # Исходный эталон является обязательной нулевой гипотезой.
    offset = fit_model('offset', 'y(t)=r(t)+b', [ones], [value - base for value, base in zip(measured, reference)], ['offset'])  # Смещение подбирается по остатку при фиксированном единичном усилении.
    if offset:  # Ошибки остатка совпадают с ошибками исходной аддитивной модели.
        models.append(offset)  # Не вводим лишний коэффициент усиления для постоянного смещения.
    for code, formula, columns, names in [('gain', 'y(t)=a*r(t)+b', [reference, ones], ['gain', 'offset']),  # Аффинная модель отделяет масштаб измерения от аппаратной динамики.
                                           ('drift', 'y(t)=a*r(t)+b+c*(t-t0)', [reference, ones, [stamp - times[0] for stamp in times]], ['gain', 'offset', 'drift_per_us']),  # Линейный дрейф описывает изменение питания или нуля по времени.
                                           ('constant', 'y(t)=b', [ones], ['level'])]:  # Постоянная модель проверяет отсутствие ожидаемого переключения.
        candidate = fit_model(code, formula, columns, measured, names)  # Линейные коэффициенты получают независимую контрольную ошибку.
        if candidate:  # Вырожденные варианты не входят в рейтинг.
            models.append(candidate)  # Все допустимые кандидаты сохраняются для полного отчёта.
    centers = [event[0] * 1000. for event in events]  # Центры эталонных переходов задают локальное растяжение без изменения паузы между импульсами.
    if centers and len(centers) <= 32 and margin > 0. and result.metrics['outside_tolerance']:  # Не подбираем сложную динамику для допустимого сигнала или более 32 неоднозначных переходов.
        def temporal(lag, stretch, code):  # Нелинейные параметры подбираются отдельно от линейных усиления и нуля.
            checkpoint(cancel)  # Каждый кандидат проверяет отмену до большого прохода интерполяции.
            warped = []  # Прогноз эталона строится на той же сетке всех глобальных моделей.
            for stamp in times:  # Сдвиг и растяжение не меняют метки фактического измерения.
                center = min(centers, key=lambda value: abs(stamp - lag - value))  # Ближайший фронт определяет локальную систему времени.
                query = center + (stamp - lag - center) / stretch  # Положительный lag задерживает измерение, stretch>1 растягивает фронт.
                value = interpolate(reference_times, result.reference.values, query / 1000.)  # Перевод мкс в мс сохраняет исходную полную сетку без копирования.
                if value is None:  # Отсутствующая опора делает именно этого кандидата недопустимым.
                    return None  # Нельзя улучшать SSE уменьшением числа контрольных точек.
                warped.append(value)  # Временная модель сохраняет все наблюдения общей области.
            return fit_model(code, 'y(t)=a*r(te+(t-delay-te)/stretch)+b', [warped, ones], measured, ['gain', 'offset'], 1 if code == 'timing' else 2, dict(delay_us=lag, stretch=stretch))  # BIC штрафует оба нелинейных параметра растяжения.
        best_delay = None  # Чистая задержка конкурирует с изменением формы фронта.
        for index in range(21):  # Грубая сетка из 21 точки охватывает ±10% общей длительности.
            candidate = temporal(-margin + 2. * margin * index / 20., 1., 'timing')  # Растяжение фиксировано единицей для отдельной модели задержки.
            if candidate and (best_delay is None or candidate['training_sse'] < best_delay['training_sse']):  # Нелинейный поиск использует только обучающие остатки.
                best_delay = candidate  # Контрольные точки не выбирают оптимальную задержку.
        if best_delay:  # При наличии перекрытия уточняем задержку последовательным уменьшением шага.
            step = margin / 10.  # Начальный шаг равен расстоянию между грубыми кандидатами.
            for _ in range(9):  # Девять уточнений дают разрешение существенно меньше исходного шага фронта.
                for lag in (best_delay['parameters']['delay_us'] - step / 2., best_delay['parameters']['delay_us'] + step / 2.):  # Симметричный поиск допускает опережение и запаздывание.
                    candidate = temporal(lag, 1., 'timing')  # Линейные параметры заново оптимизируются для каждой задержки.
                    if candidate and candidate['training_sse'] < best_delay['training_sse']:  # Контрольная ошибка остаётся независимой от решения.
                        best_delay = candidate  # Принимаем только реальное улучшение обучающего описания.
                step /= 2.  # Локальный шаг уменьшается геометрически.
            models.append(best_delay)  # Чистая задержка получает самостоятельную строку отчёта.
            best_stretch = best_delay  # Единичное растяжение служит допустимой исходной точкой.
            for stretch in (.4, .6, .8, 1., 1.3, 1.7, 2.3, 3., 4.):  # Диапазон 0,4–4 охватывает ускоренные и замедленные фронты без неограниченного переобучения.
                candidate = temporal(best_delay['parameters']['delay_us'], stretch, 'slow_switching')  # Одновременное уточнение задержки выполняется следующим координатным проходом.
                if candidate and candidate['training_sse'] < best_stretch['training_sse']:  # Выбираем старт по обучающим точкам.
                    best_stretch = candidate  # Лучший локальный масштаб определяет следующую область поиска.
            lag_step, scale_step = margin / 10., .5  # Начальные интервалы координатного уточнения заданы в мкс и безразмерном масштабе.
            for _ in range(9):  # Ограниченное число шагов сохраняет предсказуемое время работы.
                lag, stretch = best_stretch['parameters']['delay_us'], best_stretch['parameters']['stretch']  # Обе координаты уточняются относительно текущего решения.
                for next_lag, next_stretch in ((lag - lag_step, stretch), (lag + lag_step, stretch), (lag, max(.25, stretch - scale_step)), (lag, min(5., stretch + scale_step))):  # Ограничения 0,25–5 исключают вырожденное сжатие времени.
                    candidate = temporal(next_lag, next_stretch, 'slow_switching')  # Для каждой комбинации коэффициенты a и b оцениваются QR-методом.
                    if candidate and candidate['training_sse'] < best_stretch['training_sse']:  # Контрольные точки не участвуют в поиске.
                        best_stretch = candidate  # Сохраняем лучшее согласование динамики фронтов.
                lag_step /= 2.; scale_step /= 2.  # Обе сетки последовательно уточняются.
            if best_stretch['code'] == 'slow_switching':  # Совпадение с чистой задержкой не создаёт дублирующую модель.
                models.append(best_stretch)  # Штраф за дополнительное растяжение учитывается при окончательном выборе BIC.
    span = result.metrics['reference_max'] - result.metrics['reference_min']  # Масштаб состояния берётся из полного сравниваемого эталона.
    if span > 0.:  # Для постоянной полки невозможно определить противоположное состояние контакта.
        active, intervals = None, []  # Кусочная модель оценивает границы каждого дополнительного переключения.
        training = [index for index in range(len(times)) if index % 3 != 1]  # Моменты разрывов ищутся только на обучающих точках.
        for index in training:  # Принадлежность состоянию определяется относительно двух эталонных полок.
            target = result.metrics['reference_max'] if reference[index] < (result.metrics['reference_min'] + result.metrics['reference_max']) / 2. else result.metrics['reference_min']  # Противоположная полка соответствует прерыванию нормального состояния.
            contact = abs(measured[index] - target) < .3 * span and abs(measured[index] - reference[index]) > .5 * span  # Пространство состояний ограничено наблюдаемым полным размахом.
            if contact and active is None:  # Новый интервал начинается с первого подтверждённого противоположного состояния.
                active = index  # Сохраняем индекс, а не округлённое время.
            elif not contact and active is not None:  # Возврат завершает кусок неправильного состояния.
                if times[index] > times[active] and len(intervals) < 12:  # Максимум 12 сегментов ограничивает сложность и вычислительную стоимость модели.
                    intervals.append((times[active], times[index]))  # Длительность измеряется по исходной временной сетке.
                active = None  # Следующий разрыв подбирается отдельно.
        if intervals:  # Без завершённых прерываний кусочная модель отсутствует.
            boxes = [[float(left <= stamp < right) for stamp in times] for left, right in intervals]  # Границы, найденные по обучению, применяются также к контрольным временам.
            candidate = fit_model('contact_bounce', 'y(t)=a*r(t)+b+sum(cj*I[tj0,tj1])', [reference, ones, *boxes], measured, ['gain', 'offset', *[f'event_{i}_amplitude' for i in range(len(boxes))]], 2 * len(boxes), dict(events=[dict(start_us=left, stop_us=right, duration_us=right - left) for left, right in intervals]))  # BIC включает амплитуду и обе границы каждого прерывания.
            if candidate:  # Линейно зависимые сегменты не должны попадать в рейтинг.
                models.append(candidate)  # Совместный дефект может оставить звон в остатке этой кусочной модели.
    baseline = models[0]  # Исходный эталон задаёт общую точку отсчёта улучшения.
    floor = max(1e-24, max(abs(value) for value in measured) ** 2 * 1e-24)  # Все глобальные модели используют один численный пол дисперсии даже при подгонке аддитивного остатка.
    for model in models:  # Единый масштаб предотвращает искусственное преимущество модели при почти нулевой численной ошибке.
        count = model['bic_count']  # Число обучающих точек одинаково для всех глобальных конкурентов.
        model['bic'] = count * math.log(max(model['training_sse'] / count, floor)) + model['parameter_count'] * math.log(count)  # Штраф и дисперсия оцениваются на одной и той же обучающей выборке.
    for model in models:  # Все показатели сравнения моделей вычисляются на одном наборе времён.
        model['delta_bic'] = baseline['bic'] - model['bic']  # Положительное значение означает преимущество относительно неизменённого эталона.
        model['explained_error_pct'] = 100. * (1. - model['sse'] / baseline['sse']) if baseline['sse'] > 0. else 0.  # Показываем, какую долю исходной ошибки описала модель.
        model['supported'] = model['delta_bic'] >= 10. and model['validation_rmse'] < .8 * baseline['validation_rmse'] and model['explained_error_pct'] > 50.  # Решение требует штрафованного улучшения и снижения контрольной ошибки минимум на 20%.
    models.sort(key=lambda model: model['bic'])  # Рейтинг показывает конкурентов, а не единственный заранее выбранный диагноз.
    best = models[0] if models[0]['supported'] else baseline  # При недостатке доказательств сохраняем исходную гипотезу без аппаратного вывода.
    return dict(models=models, best=best['code'], sampled_points=len(times), interval_us=[times[0], times[-1]], sampling='до 2048 стратифицированных точек, включая экстремумы; контроль — каждая третья точка')  # Отчёт честно отделяет полную статистику от ограниченного подбора.


def fit_damped_ring(times_us, values, cancel=lambda: False):  # Свободный отклик второго порядка описывается экспонентой и синусом с неизвестной фазой.
    if len(values) < 24 or times_us[-1] <= times_us[0]:  # Короткая или нулевая временная область не позволяет идентифицировать колебание.
        return None  # Отсутствующая оценка не заменяется произвольной частотой.
    full_level = median(values[-max(8, len(values) // 4):])  # Начальную частоту оцениваем до прореживания, чтобы не принять алиас за реальный резонанс.
    full_amplitude = max(abs(value - full_level) for value in values)  # Полоса переходов задаётся полным наблюдаемым откликом.
    full_crossings, previous = [], 0  # Сначала используем полную временную сетку для оценки доступной частоты.
    for index, (stamp, value) in enumerate(zip(times_us, values)):  # Исходные пересечения уровня не зависят от ограниченной выборки нелинейного подбора.
        if index % 16384 == 0:  # Подготовку большого окна также можно прервать.
            checkpoint(cancel)  # Не ждём окончания длинного прохода после запроса отмены.
        sign = 1 if value - full_level > .05 * full_amplitude else -1 if value - full_level < -.05 * full_amplitude else 0  # Гистерезис подавляет поздний шум относительно начального колебания.
        if sign and previous and sign != previous:  # Противоположные значимые полуволны определяют половину периода.
            full_crossings.append(stamp)  # Время сохраняется до любой децимации.
        if sign:  # Промежуточная зона не меняет последнюю подтверждённую полуволну.
            previous = sign  # Дальнейшее пересечение сравнивается с наблюдаемым знаком.
    if len(full_crossings) < 4:  # Короткий звон или одиночный выброс не дают устойчивой начальной частоты.
        return None  # Не создаём резонансную модель без нескольких наблюдаемых полупериодов.
    initial_frequency = .5 / median(b - a for a, b in zip(full_crossings, full_crossings[1:]))  # Частота определяется до возможного понижения частоты дискретизации.
    full_step = median(b - a for a, b in zip(times_us, times_us[1:]))  # Фактический шаг исходных наблюдений задаёт допустимую децимацию.
    if full_step <= 0.:  # Невозрастающая сетка не имеет физически определённой частоты выборки.
        return None  # Не исправляем повреждённые метки времени молча.
    maximum_stride = max(1, int(.2 / (initial_frequency * full_step)))  # Сохраняем минимум пять отсчётов на исходно наблюдаемый период.
    stride = min(maximum_stride, max(1, math.ceil(len(values) / 512)))  # Ограничение памяти не должно создавать ложную низкочастотную модель из высокочастотного звона.
    times = [stamp - times_us[0] for stamp in times_us[::stride][:512]]  # Если безопасная децимация недостаточна, сокращаем окно, сохраняя разрешение периода.
    measured = list(values[::stride][:512])  # Подбор получает не более 512 реально наблюдаемых амплитуд без алиасного понижения частоты.
    duration = times[-1] - times[0]  # Длительность участка выражена в мкс.
    step = median(b - a for a, b in zip(times, times[1:]))  # Медианный шаг устойчив к отдельным неравномерным промежуткам.
    if max(b - a for a, b in zip(times, times[1:])) > 5. * step:  # Через большой ненаблюдавшийся промежуток колебание не подбирается.
        return None  # Спектральная интерпретация разорванного участка была бы неоднозначной.
    level = median(measured[-max(8, len(measured) // 4):])  # Хвост задаёт начальную оценку установившегося уровня.
    amplitude = max(abs(value - level) for value in measured)  # Оценка амплитуды нужна только для устойчивого поиска начальной частоты.
    previous, crossings = 0, []  # Переходы через небольшой гистерезис задают начальную частоту, но не сам диагноз.
    for stamp, value in zip(times, measured):  # Анализируем обе полуволны колебания относительно хвостовой полки.
        sign = 1 if value - level > .05 * amplitude else -1 if value - level < -.05 * amplitude else 0  # Полоса 5% начальной амплитуды подавляет шум после затухания.
        if sign and previous and sign != previous:  # Между противоположными полуволнами проходит приблизительно половина периода.
            crossings.append(stamp)  # Сохраняем времена для робастной начальной оценки периода.
        if sign:  # Нулевая полоса не создаёт дополнительных смен знака.
            previous = sign  # Последняя значимая полуволна остаётся опорной.
    if len(crossings) < 4:  # Менее двух наблюдаемых периодов недостаточно для уверенной идентификации затухания.
        return None  # Одиночный выброс не объявляется резонансом.
    initial = .5 / median(b - a for a, b in zip(crossings, crossings[1:]))  # Частота в циклах/мкс соответствует МГц.
    ones, best = [1.] * len(times), None  # Постоянный уровень подбирается вместе с двумя фазовыми компонентами.
    def candidate(frequency, alpha):  # Метод переменной проекции отделяет нелинейные частоту/затухание от линейных коэффициентов.
        checkpoint(cancel)  # Длительный перебор можно кооперативно отменить.
        if not 1.5 / duration <= frequency <= .2 / step or alpha < 0.:  # Требуем наблюдаемые периоды и минимум пять отсчётов на период, сохраняя запас до границы Найквиста.
            return None  # Неоднозначные или нарастающие решения этой пассивной модели отклоняются.
        cosine = [math.exp(-alpha * stamp) * math.cos(math.tau * frequency * stamp) for stamp in times]  # Первая фазовая компонента затухающего отклика.
        sine = [math.exp(-alpha * stamp) * math.sin(math.tau * frequency * stamp) for stamp in times]  # Вторая компонента допускает произвольную начальную фазу.
        return fit_model('damped_rlc', 'v(t)=Vinf+exp(-alpha*t)*(A*cos(2*pi*f*t)+B*sin(2*pi*f*t))', [ones, cosine, sine], measured, ['level', 'cosine', 'sine'], 2, dict(frequency_hz=frequency * 1e6, alpha_per_us=alpha))  # BIC учитывает пять оцениваемых параметров.
    for fi in range(17):  # Начальная частота уточняется в диапазоне ±35% оценки по переходам через уровень.
        frequency = initial * (.65 + .7 * fi / 16.)  # Ограниченный диапазон снижает риск выбора гармоники вместо основного колебания.
        for alpha in [0.] + [10. ** (-1.5 + index * 3. / 10.) / duration for index in range(11)]:  # Логарифмическая сетка охватывает слабое и быстрое затухание относительно окна наблюдения.
            model = candidate(frequency, alpha)  # Линейные параметры подбираются QR-методом для каждой пары частоты и затухания.
            if model and (best is None or model['training_sse'] < best['training_sse']):  # Контрольные точки не выбирают параметры нелинейного поиска.
                best = model  # Сохраняем лучшее объяснение обучающих наблюдений.
    if best is None:  # Частота могла оказаться выше разрешения исходной записи.
        return None  # Недостаточная дискретизация не компенсируется выдуманными параметрами.
    frequency_step, alpha_step = .05 * initial, max(best['parameters']['alpha_per_us'] * .5, 1. / duration)  # Локальное уточнение начинается с масштаба грубой сетки.
    for _ in range(10):  # Десять координатных уточнений дают воспроизводимое ограниченное время подбора.
        frequency, alpha = best['parameters']['frequency_hz'] / 1e6, best['parameters']['alpha_per_us']  # Возвращаем текущие параметры к внутренним мкс.
        for next_frequency, next_alpha in ((frequency - frequency_step, alpha), (frequency + frequency_step, alpha), (frequency, max(0., alpha - alpha_step)), (frequency, alpha + alpha_step)):  # Частота и затухание уточняются независимо с неотрицательным ограничением потерь.
            model = candidate(next_frequency, next_alpha)  # Линейные коэффициенты пересчитываются при каждом изменении динамики.
            if model and model['training_sse'] < best['training_sse']:  # Принимаем только улучшение на обучающей части.
                best = model  # Контрольный остаток остаётся показателем проверки модели.
        frequency_step /= 2.; alpha_step /= 2.  # Постепенно уменьшаем обе области поиска.
    for _ in range(4):  # Повторное точное уточнение снимает взаимное влияние ошибки частоты и оценённого затухания.
        frequency, alpha = best['parameters']['frequency_hz'] / 1e6, best['parameters']['alpha_per_us']  # Исходная точка уже найдена грубой сеткой и локальными шагами.
        refined = bounded_minimum(lambda value: candidate(value, alpha), frequency * .97, frequency * 1.03)  # Частота уточняется внутри локального диапазона ±3%.
        if refined and refined['training_sse'] < best['training_sse']:  # Золотой поиск не должен ухудшать уже полученное решение.
            best = refined  # Новая частота уменьшает систематическую ошибку фазы по всему окну.
        frequency = best['parameters']['frequency_hz'] / 1e6  # Следующее уточнение затухания использует уже исправленную частоту.
        refined = bounded_minimum(lambda value: candidate(frequency, value), 0., max(3. * alpha, 5. / duration))  # Неотрицательное затухание уточняется независимо от положения грубой логарифмической сетки.
        if refined and refined['training_sse'] < best['training_sse']:  # Сохраняем только реальное улучшение обучающего описания.
            best = refined  # Время затухания больше не наследует грубую ошибку предыдущей сетки частоты.
    params = best['parameters']  # Дополнительные физические параметры вычисляются из найденных полюсов второго порядка.
    omega = math.tau * params['frequency_hz']  # Угловая частота затухающего колебания выражается в рад/с.
    alpha = params['alpha_per_us'] * 1e6  # Переводим скорость затухания из 1/мкс в 1/с.
    natural = math.hypot(omega, alpha)  # Для недодемпфированного отклика omega0²=omega_d²+alpha².
    params.update(amplitude=math.hypot(params['cosine'], params['sine']), decay_time_us=1e6 / alpha if alpha > 0. else None, damping_ratio=alpha / natural,  # Амплитуда относится к началу выбранного участка; отсутствие потерь не создаёт Infinity в JSON.
                  natural_frequency_hz=natural / math.tau, quality_factor=natural / (2. * alpha) if alpha > 0. else None, lc_product=1. / natural ** 2, r_over_l=2. * alpha)  # LC и R/L идентифицируются без задания отдельного R, L или C для последовательного эквивалента.
    constant = fit_model('constant', 'v(t)=Vinf', [ones], measured, ['level'])  # Проверяем, что сложный отклик лучше простой постоянной полки.
    best.update(delta_bic=constant['bic'] - best['bic'], interval_us=[times_us[0], times_us[0] + times[-1]], frequency_resolution_hz=1e6 / duration, input_points=len(values), decimation_stride=stride)  # Отчёт сохраняет фактическое укороченное окно и безопасный шаг подбора.
    best['supported'] = best['delta_bic'] >= 10. and best['r_squared'] is not None and best['r_squared'] >= .8 and best['validation_rmse'] < .5 * constant['validation_rmse']  # Пригодность отклика проверяется штрафом сложности, R² и контрольной ошибкой.
    return best  # Возвращаем также неподтверждённый подбор, чтобы отчёт показывал причину отказа от гипотезы.


def equivalent_rlc(model, circuit=None):  # Параметры относятся только к явно выбранному последовательному эквиваленту второго порядка.
    if not model or not model.get('supported'):  # По плохо описанному колебанию нельзя вычислять электрические элементы.
        return dict(identifiable=False, reason='Затухающий отклик не подтверждён данными.')  # Честно сохраняем отсутствие физической оценки.
    supplied = {key: float(value) for key, value in (circuit or {}).items() if key in ('inductance_h', 'capacitance_f', 'resistance_ohm') and math.isfinite(value) and value > 0.}  # Нулевые значения означают неизвестные элементы, единицы — СИ.
    product, ratio = model['parameters']['lc_product'], model['parameters']['r_over_l']  # Найденные полюса определяют только LC и R/L без дополнительного ограничения.
    if 'inductance_h' in supplied:  # Известная индуктивность однозначно задаёт оставшиеся параметры эквивалента.
        inductance = supplied['inductance_h']  # Используем явно введённое значение в генри.
    elif 'capacitance_f' in supplied:  # Известная ёмкость позволяет восстановить L из произведения LC.
        inductance = product / supplied['capacitance_f']  # Произведение LC имеет единицу с².
    elif 'resistance_ohm' in supplied and ratio > 0.:  # При ненулевом затухании известное сопротивление определяет индуктивность.
        inductance = supplied['resistance_ohm'] / ratio  # Деление R на R/L возвращает генри.
    else:  # Одна осциллограмма без дополнительного элемента не определяет все три неизвестных.
        return dict(identifiable=False, lc_product=product, r_over_l=ratio, supplied=supplied, reason='Нужен хотя бы один известный R, L или C; по отклику определяются только LC и R/L.')  # Не выдаём произвольные значения отдельных элементов.
    estimated = dict(inductance_h=inductance, capacitance_f=product / inductance, resistance_ohm=ratio * inductance)  # Это эквивалентные параметры полюсов, а не автоматически номиналы RC-снаббера.
    mismatch = {key: (estimated[key] / value - 1.) * 100. for key, value in supplied.items()}  # Несколько известных элементов позволяют проверить согласованность предположенной топологии.
    return dict(identifiable=True, topology='последовательный эквивалент RLC', estimated=estimated, supplied=supplied, deviation_pct=mismatch)  # Отделяем введённые значения от восстановленных и не скрываем противоречия.


def fit_ringing_pair(result, events, circuit=None, cancel=lambda: False, delay_us=0.):  # Компенсируем оценённую задержку перед сравнением полюсов, чтобы сдвиг фронта не выглядел отказом снаббера.
    rising = [event for event in events if event[2] == 1][-4:]  # Не более четырёх последних выключений ограничивают стоимость нелинейного подбора.
    pairs = []  # Каждый участок имеет собственные границы, качество и параметры.
    low, high = result.metrics['reference_min'], result.metrics['reference_max']  # Огибающая взята только из реально сравниваемого эталона.
    for center, width, _ in rising:  # Время и длительность события выражены во внутренних миллисекундах.
        checkpoint(cancel)  # Смена записи может отменить подбор между отдельными выключениями.
        start = bisect_left(result.difference.times, center + width / 2.)  # Начинаем поиск установившейся полки после пересечения верхнего порога фронта.
        stop_time = min(result.difference.times[-1], result.measured.times[-1] - delay_us / 1000., center + max(40. * width, .002))  # Ограничиваем обе стороны реальными опорами после компенсации задержки без экстраполяции.
        stop = bisect_left(result.difference.times, stop_time)  # Конечный индекс не выходит за реально наблюдаемое перекрытие.
        for index in range(max(1, start), min(stop - 1, len(result.expected) - 1)):  # Ищем первый максимум эталона после завершения основного роста напряжения.
            if result.expected[index] > low + .85 * (high - low) and result.expected[index - 1] <= result.expected[index] > result.expected[index + 1]:  # Первый пик на высокой полке отделяет колебательный отклик от незавершённого фронта.
                start = index  # Неизвестная фаза синуса допускает начало подбора с любого пика.
                break  # Последующие колебания не должны сдвигать начало до почти нулевой амплитуды.
        indices = range(start, stop)  # Обе кривые подбираются на строго одинаковом временном интервале.
        if stop - start < 24 or any(not math.isfinite(result.expected[i]) or not math.isfinite(result.difference.values[i]) for i in indices):  # Разрыв внутри окна не маскируется удалением отдельных отсчётов колебания.
            continue  # Такому выключению не приписываем достоверный резонансный контур.
        times = [result.difference.times[i] * 1000. for i in indices]  # Подбор использует микросекунды без потери дробной части.
        reference = [result.expected[i] for i in indices]  # Эталон уже интерполирован на фактическую сетку измерения.
        measured = [interpolate(result.measured.times, result.measured.values, result.difference.times[i] + delay_us / 1000.) for i in indices]  # Положительная задержка сдвигает окно измерения к соответствующему выключению.
        if any(value is None for value in measured):  # Через пропуск или отсутствие опоры после сдвига отклик не подбирается.
            continue  # Окно не уменьшается отдельно для одной из сравниваемых моделей.
        fitted_reference = fit_damped_ring(times, reference, cancel)  # Затухание исправного отклика определяется по данным, а не фиксированной константой.
        fitted_measured = fit_damped_ring(times, measured, cancel)  # Частота, затухание и фаза измерения подбираются независимо.
        if fitted_measured:  # Сохраняем также неудачное по качеству согласование для объяснения неопределённости.
            pair = dict(measured=fitted_measured, reference=fitted_reference, equivalent=equivalent_rlc(fitted_measured, circuit), compensated_delay_us=delay_us)  # Эквивалент и компенсация времени сохраняются рядом с качеством подбора.
            if fitted_reference and fitted_reference['supported']:  # Без пригодного эталонного колебания нельзя оценить отношение затуханий.
                a, b = fitted_measured['parameters'], fitted_reference['parameters']  # Параметры имеют одинаковое начало выбранного участка и одинаковые единицы.
                pair.update(amplitude_ratio=a['amplitude'] / b['amplitude'] if b['amplitude'] > 1e-12 else None,  # Отношение амплитуд не зависит от единицы напряжения.
                            decay_ratio=a['decay_time_us'] / b['decay_time_us'] if a['decay_time_us'] is not None and b['decay_time_us'] else None, frequency_change_pct=(a['frequency_hz'] / b['frequency_hz'] - 1.) * 100.)  # Показываем рост длительности звона и изменение резонансной частоты относительно эталона.
            pairs.append(pair)  # Отдельные выключения остаются доступны в JSON для проверки устойчивости оценки.
    return pairs  # Отсутствие пригодных окон является явным отсутствием модели, а не нормальным состоянием снаббера.


def identify_fault_models(result, events, features, physical, role, circuit=None, channels=None, cancel=lambda: False):  # Формируем гипотезы из подтверждённых моделей, а не только из фиксированных порогов формы.
    from set_devices.signal_joint_models import fit_joint_channels  # Ленивый импорт устраняет цикл между совместными моделями и общими численными методами.
    checkpoint(cancel)  # Не начинаем дорогой подбор после отмены исходного сравнения.
    report = fit_reference_models(result, events, cancel)  # Все глобальные конкуренты оцениваются на одинаковых временах.
    best = next((model for model in report['models'] if model['code'] == report['best']), None)  # Временная модель нужна также для корректного выбора окна выключения.
    delay = best['parameters'].get('delay_us', 0.) if best and best['supported'] else 0.  # Компенсируем только подтверждённый по BIC и контролю временной сдвиг.
    rings = fit_ringing_pair(result, events, circuit, cancel, delay) if physical and role == 'uce' else []  # Электрическая интерпретация колебания допустима только для силового напряжения.
    report.update(ringing=rings, assumptions=['BIC — сравнительный критерий при предположении гауссовой ошибки, не вероятность неисправности.', 'Контрольные точки чередуются во времени; это проверка интерполяционного прогноза, не независимый эксперимент.', 'Подвыборка сохраняет экстремумы: статистические показатели условны относительно выбранных точек.', 'RLC — последовательный эквивалент полюсов; номиналы реального RC-снаббера автоматически не определяются.'])  # Ограничения метода сохраняются рядом с численными результатами.
    report['criteria'] = dict(delta_bic_min=10., validation_improvement_min_pct=20., explained_error_min_pct=50., ring_r_squared_min=.8, ring_amplitude_or_decay_ratio=1.5, ring_frequency_change_pct=15., stretch_min=1.15)  # Начальные критерии принятия модели явно доступны пользователю и не выдаются за паспортные нормы транзистора.
    report['circuit_supplied_si'] = dict(circuit or {})  # Известные элементы сохраняются даже при отсутствии пригодного колебательного участка.
    contact_events = best['parameters'].get('events', []) if best and best['code'] == 'contact_bounce' else []  # Сегменты кусочной модели проверяются независимыми физическими каналами.
    joint, joint_findings = fit_joint_channels(result, channels, contact_events, cancel) if physical else (dict(available=False, reason='Физическая роль канала не определена.'), [])  # Не выполняем электрические расчёты по произвольному неизвестному каналу.
    report['joint'] = joint  # Совместная идентификация проводимости, утечки и провалов управления экспортируется вместе с глобальным рейтингом.
    findings = list(joint_findings)  # Несколько моделей могут описывать разные участки совместного дефекта.
    if result.metrics['outside_tolerance'] == 0:  # Даже хорошее улучшение на малом шуме не отменяет пользовательский критерий допустимости.
        return report, findings  # Допустимый выбранный канал не отменяет отдельно подтверждённые совместные гипотезы.
    if best and best['supported']:  # Аппаратная формулировка требует проверки качества, а не только наличия признака.
        code, params = best['code'], best['parameters']  # Параметры выбранной модели будут показаны в подробном отчёте.
        if code == 'contact_bounce' and physical and features['measured_edges'] >= features['reference_edges'] + 2:  # Кусочное объяснение подтверждается полными переходами на всей записи.
            text = f'Модель прерываний: {len(params["events"])} эпизода; возможны плохой контакт (дребезг) или нестабильное управление затвором.'  # Длительности каждого эпизода доступны в параметрах модели.
            findings.append(dict(code=code, text=text, model=code))  # Причина ссылается на конкретный численно проверенный кандидат.
        elif code == 'slow_switching' and physical and params['stretch'] > 1.15:  # Замедление должно превышать 15% при уже доказанном модельном улучшении.
            findings.append(dict(code=code, text=f'Модель растяжения фронтов: ×{params["stretch"]:.3g}; возможны ослабленный драйвер или изменение цепи затвора.', model=code))  # Формулируем гипотезу через фактически оценённый масштаб времени.
        elif code in ('timing', 'slow_switching') and abs(params['delay_us']) > 3. * (features['sample_step_ms'] or 0.) * 1000.:  # Задержка должна превышать три исходных шага, даже если BIC чувствителен к малому шуму.
            findings.append(dict(code='timing', text=f'Модель задержки: {params["delay_us"]:+.4g} мкс; проверьте синхронизацию записей и задержку управления.', model=code))  # Отдельная причина аппаратной задержки по одному каналу не доказывается.
        elif code == 'constant' and physical and features['reference_edges']:  # Постоянная амплитуда выигрывает у динамического эталона с ожидаемыми переходами.
            state = 'высоким' if params['level'] > (result.metrics['reference_min'] + result.metrics['reference_max']) / 2. else 'низким'  # Середина эталонных состояний определяет тип отсутствующего переключения.
            text = f'Модель постоянного уровня: {result.measured.name} остаётся {state}; возможны отказ ключа, обрыв или ошибка управления.'  # Направление неисправности остаётся возможной причиной, а не доказанным компонентом.
            findings.append(dict(code='no_switching', text=text, model=code))  # Сохраняем совместимый код группы отсутствия переключения.
        elif code in ('offset', 'gain', 'drift'):  # Калибровочные и медленные изменения конкурируют с аппаратной динамикой.
            texts = dict(offset=f'Модель смещения: Δ0={params.get("offset", 0.):+.4g}; проверьте ноль датчика, общую точку и режим питания.',  # Постоянное смещение задано в единицах выбранного канала.
                         gain=f'Модель масштаба: коэффициент {params.get("gain", 1.):.5g}, смещение {params.get("offset", 0.):+.4g}; возможна ошибка измерительного тракта или изменение режима.',  # Усиление является безразмерной оценкой соответствия эталону.
                         drift=f'Модель дрейфа: {params.get("drift_per_us", 0.):+.4g} ед./мкс; возможны нестабильность питания или уход нуля измерения.')  # Производная уровня явно содержит единицу времени.
            findings.append(dict(code=code, text=texts[code], model=code))  # Сохраняем идентифицированный вид отклонения отдельно от аппаратной причины.
    gain = abs(best['parameters'].get('gain', 1.)) if best else 1.  # Общий масштаб измерительного тракта не должен ошибочно удваивать амплитуду «неисправного» снаббера.
    for pair in rings:  # Разные выключения могут подтверждать дополнительный резонансный дефект поверх дребезга.
        measured, reference = pair['measured'], pair['reference']  # Каждая сторона имеет собственные показатели качества подбора.
        ratio = pair.get('amplitude_ratio')  # Без пригодного эталонного отклика относительная амплитуда не определена.
        worse = ratio is not None and ratio / max(gain, 1e-9) > 1.5 or pair.get('decay_ratio') is not None and pair['decay_ratio'] > 1.5 or abs(pair.get('frequency_change_pct', 0.)) > 15.  # Параметры должны измениться заметно, а не только улучшить микроскопический остаток.
        if measured['supported'] and reference and reference['supported'] and worse:  # Аппаратная гипотеза требует пригодных моделей обеих осциллограмм.
            params = measured['parameters']  # Оценённые частота и затухание делают вывод проверяемым.
            text = f'Модель звона RLC: f={params["frequency_hz"] / 1e6:.3g} МГц, τ={params["decay_time_us"]:.3g} мкс; изменилось демпфирование — возможны нарушение снаббера или паразитного контура.' if params['decay_time_us'] is not None else f'Модель звона RLC: f={params["frequency_hz"] / 1e6:.3g} МГц, затухание не разрешено; проверьте снаббер и паразитный контур.'  # Нулевое оценённое затухание не превращается в бесконечность JSON.
            findings.append(dict(code='snubber', text=text, model='damped_rlc'))  # Привязываем причину к параметрической модели второго порядка.
            break  # Несколько выключений не дублируют одинаковую строку неисправности.
    if not findings:  # Недостаточно описанная форма должна оставаться неопределённой.
        code = 'noise' if features['residual_reversals'] >= 5 else 'deviation'  # Остаточные колебания обозначают помеху без приписывания конкретного электрического отказа.
        text = 'Остаток не объяснён проверенными моделями: возможны шум, помехи или неучтённая динамика.' if code == 'noise' else 'Ни одна модель не дала достаточного подтверждения причины отклонения; проверьте режим и измерительный тракт.'  # Отсутствие подходящей модели не считается доказательством исправности.
        findings.append(dict(code=code, text=text, model=None))  # Явная неопределённость сохраняется в краткой строке и JSON.
    return report, findings  # Численные результаты и текстовые гипотезы принадлежат одному снимку данных.