Статья про использование Python в научных вычислениях подтолкнула меня написать эту статью. Это история, случившаяся со мной и с коллегами 6 лет назад. На тот момент я уже достаточно подразобрался с Delphi и Python, но только теперь я ощущаю что достаточно поработал с C/C++, чтобы здраво оценить время на «ремонт» сломанного кода и вообще — общее время разработки. Да, это статья про код, который был написан разными людьми на Delphi, Python и C++ для одной и той же задачи, внутри одной команды.
Перед описанием ситуации, хочу оставить вторую преамбулу для тех программистов, которые не сталкивались с программированием в науке, и не имеют базовых представлений об этом. Осторожно, далее будет несколько мини-инфарктов:
— В коде просто обязано иметь много переменных. И многие из них должны иметь названия типа mu_e, mu_m, E, T_f, T и так далее — просто потому, что в формулы эти величины входят именно так, и учёным, которые в этой области занимаются исследованиями, могут потратить много времени вспоминая как расшифровывается название коэффициента в полном виде, зато безошибочно узнают букву, которой он обозначается во всех статьях, и скажет — «Ага! Это наша мю...»
— Код никто не защищает от неправильного ввода и нет никакого UX. Код пишется для того чтобы с ним работали учёные, аспиранты, и специально обученные студенты, и никогда — типичные_пользователи и хакеры.
— Технического задания у вас нет. Вы ведёте исследование, значит вы пишете один код, смотрите на результат, который он выдал, анализируете, и в зависимости от результатов формируется следующая задача, под которую будет дописываться существующий проект.
— Большая часть кода будет использоваться один, или два раза. Один — если аспирант, или научный сотрудник делает некий расчёт. Два — если студент делает некий расчёт, а затем демонстрирует работу кода своему научному руководителю. Соответственно, код почти никогда не абстрактен, и в нём может повстречаться костыль любого вида и размера.
На первое время хватит, начну свой рассказ. Я был слегка моложе, голова моя была чистой, сердце — горячим, а руки были холодными и чесались что-то делать. И так случилось, что одним из первых заданий в новой исследовательской группе (после смены руководителя), для меня стал разбор чужого кода, который по неясным его создателю причинам — не работал. Далее последуют три истории отладки и настройки программы. Одна из них слегка уходит в сторону размышлений, две другие про реально предпринятые действия, и столь же достоверно описывают среднестатистический опыт работы учёного с соответствующим языком программирования.
Задача была одна на всех — решение дифференциального уравнения (у которого не существует аналитического решения) с разными начальными условиями. Уравнение было большим и пёстрым — f(x) и x, описывающие состояние системы в некоем нестационарном процессе, входили в нём и в производные, и в производные их частного, и в экспоненты, и в логарифмы, и в разнообразные степенные функции. Некоторые коэффициенты при степенях неявно также зависели от них, и их необходимо было пересчитывать на каждом шаге. В дипломах и диссертациях такие уравнения встречаются довольно часто, и обычно в таких выражениях вводят замены, чтобы уместить в одну строку, но часто и замены не помогают — тогда уже расписывают на весь лист. Предсказать вид такой функции (он же — поведение системы) просто по виду задающего её уравнения — задача абсолютно неподъёмная для человека, пусть даже и человека с самыми развитыми способностями к анализу функций и самыми благими намерениями. Соответственно задача программы — тут же вывести этот самый вид, начертив график.
Решаются такие задачи просто: есть одна приблизительно известная точка, которая берётся за исходные данные. Затем от неё на каждой итерации цикла делается шаг в сторону, с расчётом новых значений из старых. В случае метода Эйлера величина шага пропорциональна производной в исходной точке. В случае метода Рунге-Кутты зависимость несколько сложнее, чтобы быть слегка точнее. Есть и другие методы, когда точка, строго говоря, вам неизвестна. Но вам известна y координата одной точки, и x координата другой. Тогда применяются слегка более извращённые методы, но сводятся они всегда к перебору всех возможных вариантов и множественному использованию Эйлера, или РК. Ну а найдя точки — надо тут же отправлять их на интерфейс в график. Звучит несложно, приступим.
Начнём с С++
После проверки всего кода, с множеством переменных с одно- и двухбуквенными названиями, вы приходите к выводу, что программа должна правильно рассчитать все коэффициенты, правильно найти значения для следующего шага, правильно перейти на следующий шаг, и правильно вывести данные на график. Следовательно, программа в целом должна отработать правильно, решаете вы, и запускаете.
Компьютер зависает и происходит какая-то очень нехорошая ошибка каким-то непонятным образом связанна с памятью.
Вы остаётесь с задачкой как это всё отладить. Очевидный путь — включить дебаггер и пройтись по потоку вычислений — в данном случае пока не подходит. Ведь вы не знаете сколько раз успел отработать цикл, вдруг миллион, два? Выходит, чтобы начать отладку, вам нужно найти точку вылета. То есть намечается целых четыре новых этапа: переписать код так, чтобы он выводил в лог номер цикла, чтобы заметить последний номер итерации, где происходит вылет, затем переписать код ещё раз, чтобы остановиться на этой итерации и выставить красную точечку в дебаггере, и пройтись по последней итерации дебаггером, надеясь на то что ошибка станет ясной. При этом
а) у вас нет гарантий что лог сохранится, в связи с падением ОС. Возможно придётся проводить работу визуально, разбив первый этап на на два: сначала выводить номер цикла в консоль и пытаться запомнить последнее число примерно, затем переписать код чтобы включить в цикл искусственный замедлитель после этой запомненной итерации — тогда на следующем тесте можно будет увидеть как числа перед падением начинают поступать медленнее, и успеть прочитать последнюю итерацию (скажем, 6524-ая), значит количество падений будет не 3, а 4.
б) у вас нет гарантий, что вы поймёте, или даже найдёте ошибку с первого раза. Скорее всего где-то перед самой ошибкой код начнёт уводить вас вглубь кроличьей норы — по извилистым лабиринтам стандартных библиотек, и хоть вы пройдётесь всего по одной ветке, вам потребуется изучить значительную часть дерева возможностей вокруг этой ветки, чтобы понять что вообще произошло, ведь функции в ней зачастую просто перебрасывают какие-то ошмётки данных от одних проверок к другим, для понимания многих из них вам понадобиться ещё и разбираться в физическом устройстве машинных операций, а понимание значения иных приходит к человеку только после попытки написать свою ОС. Более того, после того как вы решите что нашли в чём проблема, нет никаких гарантий, что вы правы. Начался длительный процесс отладки, и вам предстоит ещё много тестов… Вместо четырёх падений мы получаем n, где n>3.
в) необходимые для тестов падения, вообще говоря, вещь довольно стрёмная. Мало ли что произойдёт? А тем более, когда предстоит сделать это много раз.
Теперь вы сидите и думаете, а стоит ли вообще приступать к реализации такого плана? Он выглядит громоздко, и неизвестность вас слегка пугает — мало ли где зарылась ошибка.
Pascal/Delphi
Проверка кода. Вообще говоря она может занимать больше времени, но в данном конкретном случае вам нравится, что все 66 коэффициентов объявлены и заданы в начале, а не в любом месте кода, например перед использованием. Вы сверили коэффициентики, отбросили табличку и теперь просто сверили все формулы. Опять таки, всё хорошо. Запускаете программу, и — о чудо, график рисуется. График рисуется, доходит до конца оставленного для него участка числовой оси и останавливается. Программа не вылетела, ничего страшного с ОС не произошло. Но произошло с данными. Он плавно доходит до координаты 6524, а затем вдруг его кусает бешеная блоха. Он скачет вверх и вниз абсолютно случайно и заполняет всё оставшееся полупространство справа кромешной тьмой.
На этом этапе вы удручены неудачей. Но всё же находитесь в выигрышной позицией по сравнению с си, ведь:
а) система не упала
б) вы уже локализовали итерацию, где возникает ошибка
в) вы можете пробовать новые исправления, не боясь что система вдруг упадёт
и г) — вы видите одно выведенное значение. Оно неверно. И возможно, уже на этом этапе предполагаете что это может быть, скажем, переполнение, либо ошибка чисел с плавающей точкой. Осталось понять где именно, чтобы не создавать тысячи экстендедов. Вы дописываете строчку кода чтобы остановить на номере 6524. Проходитесь дебаггером и замечаете подозрительную операцию очень маленького числа на очень большое, начиная с которой тянется ошибка, которая приводит к расходимости. Есть оказывается, в вашей функции такая особая точка, крохотный участок, где состояние имеет локальный минимум, но за пределами небольшого ограждения в потенциале начинается крутой спуск, приводящий к расходимости всего решения, и уводящий значения за пределы и дальше — в переполнение их типов. И естественно, численный метод немного промахивается на этом шаге мимо нужного локального минимума, из-за банальной нехватки точности для всего двух коэффициентов. И заметно это становится только здесь — в том узком месте, где даже незначительно отклонение отправляет итерации по дурной дорожке расходимости. Значит надо ещё разок переписать код, ввернув пару экстендед вариэйблс… Но опять какой-то косяк происходит, хотя график выглядит уже немного по другому. Я имею ввиду, конечно же, чёрную его часть. Брыкающуюся. Сколько ещё надо колдовать над кодом? Ну, поработать над реализацией работы с длинными числами. И обнаружить (после исправления этих двух) ещё один «протекающий»* коэффициент.
*протекающий — это потому что выход за пределы типа должен приводить к изменению бита слева от места его физического хранения, т. е. Нам нужен ещё один, либо несколько бит, которые не содержаться по нашей ссылке на объект такого типа. Ок, это заняло не так уж много времени.
перепишем на Python
Проверяем код, запускаем. Проблема не наблюдается. График чертится полностью, на всей области. Никаких скачков и черноты, никаких вылетов.
Волшебство интерпретатора помогает доброму питону увидеть, сколько памяти требует моё число, и выделить ровно столько, сколько ему будет достаточно.
Разработка закончена.
Внимательный читатель может заметить, что вскоре мы, учёные, можем захотеть выполнить этот код перебирая параметры, и поместив его во вложенный цикл. И тогда скорость выполнения Python нанесёт нам удар ножом в спину. Однако должен заметить, что код для других значений параметров мог иметь слабые места в других точках, на других рассчитанных коэффициентах. И тогда у человека есть выбор — либо оставить код на Python работать всю ночь ради одного графика. Либо оставить код на С++ работать пока он сходит сделать себе чай. Но когда он вернётся — есть колоссальный шанс увидеть зависание системы, или ошибки другого рода. Он сядет их исправлять для нового слабого места, потом потратит ещё один поход за чаем впустую, в итоге плюнет и начнёт писать более абстрактный код, который подойдёт его растущим требованиям к программе. Вот только написать абстрактный код, который следит за переменными в потоках самых разнообразных вычислений и даёт им больше памяти при необходимости — не забывая в итоге приводить типы к менее ёмким, чтобы не засорять память — это значит написать значительную часть интерпретатора Python.
Так что я с тех самых пор пользуюсь именно Python’ом почти для всех научных работ (одну я всё-таки сделал в спец пакете, и только ради красивой визуализации под случай), и даже разобрался в некоторых его тонкостях. Чего и вам желаю.
Используйте Python. Двигайте науку. Лучей добра.
Перед описанием ситуации, хочу оставить вторую преамбулу для тех программистов, которые не сталкивались с программированием в науке, и не имеют базовых представлений об этом. Осторожно, далее будет несколько мини-инфарктов:
— В коде просто обязано иметь много переменных. И многие из них должны иметь названия типа mu_e, mu_m, E, T_f, T и так далее — просто потому, что в формулы эти величины входят именно так, и учёным, которые в этой области занимаются исследованиями, могут потратить много времени вспоминая как расшифровывается название коэффициента в полном виде, зато безошибочно узнают букву, которой он обозначается во всех статьях, и скажет — «Ага! Это наша мю...»
— Код никто не защищает от неправильного ввода и нет никакого UX. Код пишется для того чтобы с ним работали учёные, аспиранты, и специально обученные студенты, и никогда — типичные_пользователи и хакеры.
— Технического задания у вас нет. Вы ведёте исследование, значит вы пишете один код, смотрите на результат, который он выдал, анализируете, и в зависимости от результатов формируется следующая задача, под которую будет дописываться существующий проект.
— Большая часть кода будет использоваться один, или два раза. Один — если аспирант, или научный сотрудник делает некий расчёт. Два — если студент делает некий расчёт, а затем демонстрирует работу кода своему научному руководителю. Соответственно, код почти никогда не абстрактен, и в нём может повстречаться костыль любого вида и размера.
На первое время хватит, начну свой рассказ. Я был слегка моложе, голова моя была чистой, сердце — горячим, а руки были холодными и чесались что-то делать. И так случилось, что одним из первых заданий в новой исследовательской группе (после смены руководителя), для меня стал разбор чужого кода, который по неясным его создателю причинам — не работал. Далее последуют три истории отладки и настройки программы. Одна из них слегка уходит в сторону размышлений, две другие про реально предпринятые действия, и столь же достоверно описывают среднестатистический опыт работы учёного с соответствующим языком программирования.
Задача была одна на всех — решение дифференциального уравнения (у которого не существует аналитического решения) с разными начальными условиями. Уравнение было большим и пёстрым — f(x) и x, описывающие состояние системы в некоем нестационарном процессе, входили в нём и в производные, и в производные их частного, и в экспоненты, и в логарифмы, и в разнообразные степенные функции. Некоторые коэффициенты при степенях неявно также зависели от них, и их необходимо было пересчитывать на каждом шаге. В дипломах и диссертациях такие уравнения встречаются довольно часто, и обычно в таких выражениях вводят замены, чтобы уместить в одну строку, но часто и замены не помогают — тогда уже расписывают на весь лист. Предсказать вид такой функции (он же — поведение системы) просто по виду задающего её уравнения — задача абсолютно неподъёмная для человека, пусть даже и человека с самыми развитыми способностями к анализу функций и самыми благими намерениями. Соответственно задача программы — тут же вывести этот самый вид, начертив график.
Решаются такие задачи просто: есть одна приблизительно известная точка, которая берётся за исходные данные. Затем от неё на каждой итерации цикла делается шаг в сторону, с расчётом новых значений из старых. В случае метода Эйлера величина шага пропорциональна производной в исходной точке. В случае метода Рунге-Кутты зависимость несколько сложнее, чтобы быть слегка точнее. Есть и другие методы, когда точка, строго говоря, вам неизвестна. Но вам известна y координата одной точки, и x координата другой. Тогда применяются слегка более извращённые методы, но сводятся они всегда к перебору всех возможных вариантов и множественному использованию Эйлера, или РК. Ну а найдя точки — надо тут же отправлять их на интерфейс в график. Звучит несложно, приступим.
Начнём с С++
После проверки всего кода, с множеством переменных с одно- и двухбуквенными названиями, вы приходите к выводу, что программа должна правильно рассчитать все коэффициенты, правильно найти значения для следующего шага, правильно перейти на следующий шаг, и правильно вывести данные на график. Следовательно, программа в целом должна отработать правильно, решаете вы, и запускаете.
Компьютер зависает и происходит какая-то очень нехорошая ошибка каким-то непонятным образом связанна с памятью.
Вы остаётесь с задачкой как это всё отладить. Очевидный путь — включить дебаггер и пройтись по потоку вычислений — в данном случае пока не подходит. Ведь вы не знаете сколько раз успел отработать цикл, вдруг миллион, два? Выходит, чтобы начать отладку, вам нужно найти точку вылета. То есть намечается целых четыре новых этапа: переписать код так, чтобы он выводил в лог номер цикла, чтобы заметить последний номер итерации, где происходит вылет, затем переписать код ещё раз, чтобы остановиться на этой итерации и выставить красную точечку в дебаггере, и пройтись по последней итерации дебаггером, надеясь на то что ошибка станет ясной. При этом
а) у вас нет гарантий что лог сохранится, в связи с падением ОС. Возможно придётся проводить работу визуально, разбив первый этап на на два: сначала выводить номер цикла в консоль и пытаться запомнить последнее число примерно, затем переписать код чтобы включить в цикл искусственный замедлитель после этой запомненной итерации — тогда на следующем тесте можно будет увидеть как числа перед падением начинают поступать медленнее, и успеть прочитать последнюю итерацию (скажем, 6524-ая), значит количество падений будет не 3, а 4.
б) у вас нет гарантий, что вы поймёте, или даже найдёте ошибку с первого раза. Скорее всего где-то перед самой ошибкой код начнёт уводить вас вглубь кроличьей норы — по извилистым лабиринтам стандартных библиотек, и хоть вы пройдётесь всего по одной ветке, вам потребуется изучить значительную часть дерева возможностей вокруг этой ветки, чтобы понять что вообще произошло, ведь функции в ней зачастую просто перебрасывают какие-то ошмётки данных от одних проверок к другим, для понимания многих из них вам понадобиться ещё и разбираться в физическом устройстве машинных операций, а понимание значения иных приходит к человеку только после попытки написать свою ОС. Более того, после того как вы решите что нашли в чём проблема, нет никаких гарантий, что вы правы. Начался длительный процесс отладки, и вам предстоит ещё много тестов… Вместо четырёх падений мы получаем n, где n>3.
в) необходимые для тестов падения, вообще говоря, вещь довольно стрёмная. Мало ли что произойдёт? А тем более, когда предстоит сделать это много раз.
Теперь вы сидите и думаете, а стоит ли вообще приступать к реализации такого плана? Он выглядит громоздко, и неизвестность вас слегка пугает — мало ли где зарылась ошибка.
Pascal/Delphi
Проверка кода. Вообще говоря она может занимать больше времени, но в данном конкретном случае вам нравится, что все 66 коэффициентов объявлены и заданы в начале, а не в любом месте кода, например перед использованием. Вы сверили коэффициентики, отбросили табличку и теперь просто сверили все формулы. Опять таки, всё хорошо. Запускаете программу, и — о чудо, график рисуется. График рисуется, доходит до конца оставленного для него участка числовой оси и останавливается. Программа не вылетела, ничего страшного с ОС не произошло. Но произошло с данными. Он плавно доходит до координаты 6524, а затем вдруг его кусает бешеная блоха. Он скачет вверх и вниз абсолютно случайно и заполняет всё оставшееся полупространство справа кромешной тьмой.
На этом этапе вы удручены неудачей. Но всё же находитесь в выигрышной позицией по сравнению с си, ведь:
а) система не упала
б) вы уже локализовали итерацию, где возникает ошибка
в) вы можете пробовать новые исправления, не боясь что система вдруг упадёт
и г) — вы видите одно выведенное значение. Оно неверно. И возможно, уже на этом этапе предполагаете что это может быть, скажем, переполнение, либо ошибка чисел с плавающей точкой. Осталось понять где именно, чтобы не создавать тысячи экстендедов. Вы дописываете строчку кода чтобы остановить на номере 6524. Проходитесь дебаггером и замечаете подозрительную операцию очень маленького числа на очень большое, начиная с которой тянется ошибка, которая приводит к расходимости. Есть оказывается, в вашей функции такая особая точка, крохотный участок, где состояние имеет локальный минимум, но за пределами небольшого ограждения в потенциале начинается крутой спуск, приводящий к расходимости всего решения, и уводящий значения за пределы и дальше — в переполнение их типов. И естественно, численный метод немного промахивается на этом шаге мимо нужного локального минимума, из-за банальной нехватки точности для всего двух коэффициентов. И заметно это становится только здесь — в том узком месте, где даже незначительно отклонение отправляет итерации по дурной дорожке расходимости. Значит надо ещё разок переписать код, ввернув пару экстендед вариэйблс… Но опять какой-то косяк происходит, хотя график выглядит уже немного по другому. Я имею ввиду, конечно же, чёрную его часть. Брыкающуюся. Сколько ещё надо колдовать над кодом? Ну, поработать над реализацией работы с длинными числами. И обнаружить (после исправления этих двух) ещё один «протекающий»* коэффициент.
*протекающий — это потому что выход за пределы типа должен приводить к изменению бита слева от места его физического хранения, т. е. Нам нужен ещё один, либо несколько бит, которые не содержаться по нашей ссылке на объект такого типа. Ок, это заняло не так уж много времени.
перепишем на Python
Проверяем код, запускаем. Проблема не наблюдается. График чертится полностью, на всей области. Никаких скачков и черноты, никаких вылетов.
Волшебство интерпретатора помогает доброму питону увидеть, сколько памяти требует моё число, и выделить ровно столько, сколько ему будет достаточно.
Разработка закончена.
Внимательный читатель может заметить, что вскоре мы, учёные, можем захотеть выполнить этот код перебирая параметры, и поместив его во вложенный цикл. И тогда скорость выполнения Python нанесёт нам удар ножом в спину. Однако должен заметить, что код для других значений параметров мог иметь слабые места в других точках, на других рассчитанных коэффициентах. И тогда у человека есть выбор — либо оставить код на Python работать всю ночь ради одного графика. Либо оставить код на С++ работать пока он сходит сделать себе чай. Но когда он вернётся — есть колоссальный шанс увидеть зависание системы, или ошибки другого рода. Он сядет их исправлять для нового слабого места, потом потратит ещё один поход за чаем впустую, в итоге плюнет и начнёт писать более абстрактный код, который подойдёт его растущим требованиям к программе. Вот только написать абстрактный код, который следит за переменными в потоках самых разнообразных вычислений и даёт им больше памяти при необходимости — не забывая в итоге приводить типы к менее ёмким, чтобы не засорять память — это значит написать значительную часть интерпретатора Python.
Так что я с тех самых пор пользуюсь именно Python’ом почти для всех научных работ (одну я всё-таки сделал в спец пакете, и только ради красивой визуализации под случай), и даже разобрался в некоторых его тонкостях. Чего и вам желаю.
Используйте Python. Двигайте науку. Лучей добра.