Есть код, который физики передают друг другу с 2002 года. Кристиан Мэтцлер из Бернского университета написал набор MATLAB‑функций для расчёта рассеяния Ми, физики, из‑за которой небо голубое, а молоко белое. В 2009-м Стивен Жак обернул его для тканевой оптики, и с тех пор эти полторы сотни строк разошлись по тысячам работ, от атмосферных исследований до лазерной медицины.
Я перенёс этот код на Python. Вместе со всей проверочной обвязкой ушло несколько дней. Сам перевод в этой истории скучный: любая языковая модель делает его за минуту. Интересный вопрос начинается дальше. Откуда известно, что новый код считает то же самое? Не «кривые на графике вроде совпадают», а строго, с числом совпавших знаков и чьей‑то фамилией под ним.
Ниже: как выглядит дифференциальное тестирование на практике, как правило округления NumPy едва не сломало порт молча и почему двадцатилетняя библиотека возвращает NaN в одном конкретном случае и почему этот NaN правильный.
Доверие не переносится вместе с кодом
Схема проверки старая и простая. Прогнать оба кода на одинаковых входах и сравнить выходы. Называется differential testing. Всё сложное живёт в деталях.
Первая деталь — корпус входов. Я собрал 714 тестовых случаев: примеры из документации и демо‑скриптов, граничные ситуации и большой случайный свип с зафиксированным сидом. Диапазон брал с запасом, чтобы потом не оправдываться: размерный параметр x идёт от 0.05 до 150, частицы могут быть поглощающими (комплексный показатель преломления), а один случай сидит в вырожденной точке, где частица оптически неотличима от среды. У каждого файла корпуса в манифесте лежит SHA-256, так что подменить входы задним числом ради красивого результата не выйдет.
Вторая деталь — допуски. У каждой выходной величины свой порог, и у каждого порога письменное обоснование. Не «1e-9, потому что проходит», а «1e-9, потому что вот механизм, который накапливает ошибку». Примеры ниже.
Третья — правило прохождения. Элемент проходит, если укладывается в относительный порог или в абсолютный. Это «или» не от лени. В прозрачной среде поглощение равно ровно нулю, и чисто относительная метрика делит на ноль и улетает в бесконечность именно тогда, когда обе реализации совпали идеально. Нули сравнивают по абсолютной шкале. В отчёте каждый статус говорит, какая ветка его вытащила: PASS(rel), PASS(abs), PASS(both).
Результат, чтобы сразу закрыть вопрос: 714 случаев, ноль провалов, совпадение примерно до девяти значащих цифр в худшей точке. Теперь про места, где эти цифры чуть не потерялись.
Ловушка первая: round округляет не туда, куда вы думаете
Число членов ряда в коде Мэтцлера задаётся формулой прямиком из Бохрена и Хаффмана:
matlab
nmax = round(2 + x + 4*x^(1/3));
Выглядит безобидно. Но MATLAB round толкает половинки от нуля, и round(2.5) даёт 3. NumPy округляет к чётному, и numpy.round(2.5) даёт 2. Пока аргумент не попадает ровно на половинку, ничего не происходит. А когда попадает, две реализации строят ряды разной длины, и сравнивать становится нечего: массивы коэффициентов не совпадают даже по форме.
Поэтому в корпусе есть семейство фикстур edge_round_half, где входы подобраны так, чтобы 2 + 4x^(1/3) + x садилось точно на половинку. В порте вместо numpy.round стоит явное округление от нуля, а совпадение nmax работает жёстким гейтом: прежде чем сравнивать хоть какие‑то числа, харнесс проверяет, что обе стороны договорились о длине ряда. Эти самые edge_round_half оказались худшим случаем в половине таблиц итогового отчёта, что неплохо оправдывает привычку тестировать границы намеренно, а не в надежде.
Ловушка вторая: Бесселю дрейфовать можно, вам нельзя
Коэффициенты Ми считаются через сферические функции Бесселя полуцелого порядка, а у поглощающих частиц аргумент ещё и комплексный. MATLAB и SciPy реализуют их по‑разному, и на порядках около 175 (это x под 150) расхождение доходит примерно до 2e-9 в относительных величинах. Это ничей не баг: две корректные реализации одной спецфункции просто расходятся в последних битах.
Отсюда разные допуски у разных выходов. Коэффициенты an, bn, cn, dn получают 1e-8, с запасом над измеренным дрейфом Бесселя. Эффективности qext и qsca держат 1e-11, потому что собираются суммированием положительных членов без сокращений, и дрейф в хвосте ряда тонет в сумме. Их разность qabs = qext минус qsca ослаблена до 1e-9: вычитание двух почти равных чисел съедает значащие цифры всякий раз, когда поглощение мало на фоне экстинкции. Каждый порог в tolerances.yaml несёт свою строчку‑обоснование, файл открытый, так что спорить можно с любой строкой.
Ловушка третья, любимая: честный NaN лучше уверенного мусора
Задайте частице показатель преломления среды. m = 1, рассеивать нечего. Физически анизотропия рассеяния в этой точке не определена, примерно как средняя температура по пустой палате. А код всё равно обязан что‑то вернуть.
MATLAB доходит до NaN честной дорогой. Суммы коэффициентов схлопываются в точный ноль, qsca равен нулю, деление asy/qsca превращается в 0/0, то есть в NaN. Python до нуля не добирается. Его qsca андерфлоит примерно до 1e-32, и деление возвращает конечное число, уверенное, аккуратное и абсолютно бессмысленное, потому что это шум с плавающей точкой, делённый на другой такой же шум.
Порт по умолчанию воспроизводит NaN: ниже порога qsca в 1e-20 анизотропия объявляется неопределённой, как в оригинале. Конечное значение доступно за флагом fix_nonscattering, и когда флаг включён, затронутые строки отчёта получают статус DEVIATION с именем флага рядом. Это общая политика всего проекта. По умолчанию bug‑for‑bug, любое улучшение за флагом и с пометкой в таблицах, а не в сноске. Тот, кто мигрирует двадцатилетний код, хочет сначала точную копию, а разговор о починках уже потом.
Археология: что ещё было в оригинале
В проверенном диапазоне я не нашёл ни одной ошибки в самих числах, и это тоже результат. Ещё в 2017-м Роб Браун прошёлся по версии 2009 года и убрал то, на чём она спотыкалась: демо‑скрипт падал на необъявленной переменной, а вызов Mie с большой буквы ломался на регистрозависимых файловых системах. До меня от той эпохи доехали только окаменелости в комментариях.
Вещи потоньше нашлись, все они лежат в файле находок.
Библиотечная функция печатает в stdout на каждом вызове, так что свип из пары сотен точек превращает консоль в водопад. Порт по умолчанию молчит (verbose=False), и это помечено как отклонение по вводу‑выводу, а не по числам.
Докстринг обещает, что показатель преломления среды может быть комплексным. Проверка говорит, что не может. Комплексный nmed делает комплексным размерный параметр x, после чего сравнение x > 0 в MATLAB тихо смотрит только на вещественную часть, а round от комплексного выдаёт мусорный nmax. Задокументированная возможность, которая никогда не работала, стоит дороже обычного бага: теперь граница применимости написана явно.
Для отрицательного x в коде нет ни одной ветки. Функция просто не присваивает результат, и MATLAB падает с невнятным сообщением. Порт поднимает ValueError с читаемым. Он воспроизводит сам факт отказа, не воспроизводя невнятность.
Что я из этого вынес
Репозиторий открыт: Python‑порт, весь корпус с хешами, tolerances.yaml с обоснованиями, отчёт и скрипты, которые перегенерируют всё это у вас на машине, включая MATLAB‑сторону, если есть лицензия, и Octave, если нет. Отчёт читается сверху вниз: сводка для тех, кто решает, и таблицы худших случаев для тех, кто проверяет.
Главный вывод у меня получился не про Python и не про MATLAB. Перевод кода стал дешёвым. Доказательство эквивалентности осталось дорогим, и теперь именно оно отделяет «код переписали» от «коду можно верить». Формула nmax с её половинками, андерфлоу вместо нуля и Бессель, гуляющий в девятом знаке, невидимы на глаз и невидимы на графике. Видимыми их делает корпус из 714 случаев с порогами, за которыми стоят причины.
Репозиторий на Гитхабе. Найдёте дыру в методологии, заводите issue. Ради этого он и лежит открытым.