Julia language

Фармакокинетическая модель это часто упрощенное математическое описание того, как изменяется концентрация(количество) исследуемого вещества во времени в биологической системе. Упрощенное потому, что организм сложный механизм в котором абсорбция, распределение, метаболизм и выведение веществ (ADME) может происходить различными путями в зависимости, как от самого вещества, так и от текущего состояния организма и не редко с не очевидными обратными связями. Поэтому отразить и описать все возможные факторы, действующие на кинетику вещества практически невозможно и в большинстве случаев не целесообразно. Тем не менее, возможно представление, в котором вещество распределено в ограниченных областях организма, которые будут далее называться камерами, а переход вещества из одной камеры в другую описывается определенным уравнением.

Фармакокинетическое моделирование можно рассматривать как частный случай нелинейного моделирования. Для этого широко используются коммерческие пакеты, таких как: Phoenix WinNonlin, NONMEM, Monolix и т.д., которые позволяют получить оценки параметров модели разной (произвольной) сложности. Существуют также пакеты доступные для среды вычислений R Project (nlme, nls, saemix) и Julia (CurveFit.jl, LsqFit.jl, NLreg.jl), которые можно использовать для определения параметров фармакокинетических моделей. Ограничения пакетов R Project и Julia связаны с тем, что исследователь должен задать нелинейную функцию в явном виде, что не всегда возможно. Кроме того, не все указанные пакеты обладают достаточным функционалом для оценки модели со смешанными эффектами (случайными и фиксированными эффектами).

Несколько особняком стоит Pumas.jl, т.к. несмотря на серьезный заявленный функционал пакет не получил широкого распространения (возможно пока), вероятно этому (широкому распространению) мешает невнятная лицензия использования. Цель данной работы: описание основных подходов к фармакокинетическому моделированию, построение фармакокинетической модели в виде системы дифференциальных уравнений, подбор параметров фармакокинетической модели.

Модели можно классифицировать по количеству камер, способу поступления и интервалу дозирования вещества, механизму(функции) элиминации:

1. Количество камер:
   
   • Однокамерные
   • Двухкамерных
   • Трехкамерные (и более)
2. Способ поступления:
   
   • Внутривенно болюсно
   • Внутривенно инфузионно
   • Внутрь (per os)
3. Интервал дозирования:
   
   • Однократно
   • Многократно
4. Механизм элиминации:
   
   • Линейный
   • Согласно уравнению Михаэлиса-Ментен

Простейшая однокамерная модель представляет собой одну камеру, в которой распределено вещество, в этом случае считается, что вещество быстро распределяется: относительно быстро переносится из кровотока в ткани и обратно, т.е. ограничивающее влияние периферических камер незначительно и не учитывается. Представим, что вещество уже попало и распределилось в организме, тогда изменение его концентрации во времени можно представить в виде следующего дифференциального уравнения 2 (однокамерная модель для болюсного внутривенного введения), при решении которого будет получена функция зависимости концентрации от времени 3. Где: A количество вещества, K константа, определяющая переход из одной камеры в другую (в частном случае Kel константа элиминации), t время, отношение концентрации С и количества вещества A описывается выражением 1, где V объем распределения.

Более интересный случай внесосудистое введение. В этом случае вещество сперва должно попасть в центральную камеру из камеры, которая отражает желудочно-кишечный тракт (Vg). Данную модель можно представить в виде системы дифференциальных уравнений 4. Решение относительно количества вещества в центральной камере представлено уравнением 5. Здесь появляется дополнительная константа F биодоступность, которая отражает долю вещества, перешедшего в системный кровоток из ЖКТ. Используя выражение 1 можно вывести уравнение зависимости концентрации C от времени t уравнение 6.

Еще более сложный (приближенный к реальности) вариант двухкомпартментная модель. В этом случае кроме центральной камеры в модель включена периферическая камера дополнительный объем, при этом перенос из центральной камеры в периферическую и обратно может происходить с разными скоростями (к примеру, липофильное вещество может из кровотока довольно быстро переходить в жировую ткань, а обратно значительно медленнее). Соответствующая система двух дифференциальных уравнений для болюсного внутривенного введения выражение 7, для приема внутрь выражение 8. Система уравнений 7 составлена относительно концентрации в камерах, что также справедливо. Так как измерить количество вещества в камере очень затруднительно, то построение модели относительно концентрации позволяет связать модель и фактические измерения. В выражении 8 первое уравнение приведено для количества вещества A, а последующие для концентрации с использованием константы F (биодоступности) и объема распределения V, что не должно вызывать смущения.

Решение этого уравнения возможно, но громоздко. Поэтому решать его будем при помощи DiffEquations.jl в Julia. Для чего определим функцию pkf! как отражение вышеуказанной системы дифференциальных уравнений. Значения начальных параметров определяет начальное состояние системы для нас это доза вещества с которой начался процесс (в камере ЖКТ 3.0 условных единицы, в камерах 1, 2 ноль)

#Функцияfunction pkf!(du, u, p, t)  K = p[1]  K= p[2]  K= p[3]  K = p[4]  A = u[1]  C = u[2]  C = u[3]  Vf = p[5]  du[1] = - K  * A  du[2] =   K  * A / Vf   - K * C  - K * C  + K * C  du[3] =   K * C - K * Cend#Начальное значениеu0    = [3.0, 0.0, 0.0]#Параметрыp     = [0.8, 0.35, 0.20, 0.15, 0.8]#Отрезок времениtspan = (0.0, 50.0)#ODE problemprob  = ODEProblem(pkf!, u0, tspan, p)#Решениеsol   = solve(prob)#Графикиplot(sol)

При глубоком изучении может потребоваться построение моделей с несколькими камерами, с различными путями элиминации, а также с включением механизмов обработки событий. В таких случаях прямое решение может быть либо очень сложным, либо невозможным. В это время вышеуказанный подход позволяет получать решения уравнений с различными модификациями довольно быстро. К примеру, эту систему можно дополнить дополнительной камерой, которая будет накапливать вещество и в определенные моменты времени быстро выбрасывать его в начальную камеру имитируя энтерогепатическую рециркуляцию (Event Handling and Callback Functions). Возможно включение фармакодинамической части, влияющей на объем распределения определенной камеры или части описывающей кинетику растворения твердой лекарственной формы. Все эти модификации могут быть гибко и быстро учтены с помощью DiffEquations.jl и крайне сложно решаются аналитически.

В реальной жизни мы не знаем параметров модели, но можем располагать данными о концентрации исследуемых веществ, на основе которых можно подобрать наиболее подходящие параметры модели. Практические задачи фармакокинетического моделирования можно свести к определению параметров индивидуальных фармакокинетических моделей и их последующему анализу с помощью известных статистических методов и к определению популяционных фармакокинетических моделей с учетом изучения влияния различных интересующих факторов (пол, возраст, раса, сопутствующие заболевания, стадия течения заболевания и т.д.).

Для примера мы приведем способ подбора индивидуальной двухкомпартментной фармакокинетической модели для однократного внесосудистого применения лекарственного вещества.

Генерация индивидуальных данных с использованием решения выше:

x = [0.25, 0.5, 1.0, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0,  8.0, 12.0, 24.0, 36.0, 48.0, 60.0]y = hcat(sol.(x)...)'[:,2] .* exp.(rand(Normal(), length(x)) ./ 8)

Функция потерь необходима для оптимизации (минимизации). Обычно используется сумма квадратов отклонений реальных данных от модели, но могут быть использованы и другие. Чем меньше значение функции потерь, тем лучше подгонка (см. также методы максимального правдоподобия).

function loss_function(sol, x, y)   tot_loss = 0.0   if any((s.retcode != :Success for s in sol))     tot_loss = Inf   else     tot_loss = sum(value(L2DistLoss(), hcat(sol.(x)...)'[:,2], y))     #tot_loss = sum(value(LogitDistLoss(), hcat(sol.(x)...)'[:,2], y))     #tot_loss = sum(value(HuberLoss(mean(y)), hcat(sol.(x)...)'[:,2], y))   end   tot_lossend

Оптимизация выполняется с помощью Optim.jl, для чего должна быть определена cost_function (подробно тут: Parameter Estimation and Bayesian Analysis).

cost_function = build_loss_objective(prob, Tsit5(), f ->  loss_function(f, x, y))result = optimize(cost_function, [0.4, 0.1, 0.1, 0.1, 1.0], NelderMead(), Optim.Options(allow_f_increases = true, iterations = 5000, time_limit = 60.0))#Можно пробовать разные методы:#result = optimize(cost_function, [0.4, 0.1, 0.1, 0.1, 1.0], Newton(), Optim.Options(allow_f_increases = true, iterations = 5000, time_limit = 60.0))

Проверка решения:

prob   = ODEProblem(pkf!, u0, tspan, result.minimizer)sol    = solve(prob)plot!(sol,vars=(2))

Возможные затруднения в решении связаны с тем, что сложная модель может не сходиться или сходиться в локальных минимумах, в связи с чем оценка будет не верной. Что можно сделать?

1. Адекватно оценить сложность модели и при необходимости упростить.
2. Попробовать другие методы (Ньютона, LBFGS, Рой частиц).
3. Изменить начальные значения. Можно провести начальную оценку внемодельными методами (NCA), упрощенной моделью, методом Монте-Карло.
4. Задать интервалы параметров для оптимизации.

Послесловие.

Фармакокинетическое моделирование может выполняться с использованием двух основных подходов. Первый случай полные данные. При этом данные каждого субъекта позволяют построить индивидуальную модель (как в рассмотренном примере). Для этого требуетcя, что бы от каждого субъекта было получено достаточное количество качественных наблюдений (под качественными наблюдениями понимаются такие наблюдения, которые получены строго по плану исследования с минимальными и зарегистрированными погрешностями). Примером таких данных служат клинические исследования I фазы, когда для каждого субъекта исследования получен полный фармакокинетических профиль. Как правило, в таких исследованиях регистрируются более 95% данных от планируемого, отклонения по времени составляю не более 5 минут, а погрешности биоаналитики не превышают 15% (обычно менее 5%). После построения моделей уже индивидуальные параметры модели для каждого субъекта могут подвергаться статистическому анализу. Несмотря на хорошее качество данных, могут возникать ситуации, когда для определенных субъектов модель не сможет быть построена, что создаст проблемы для целостности анализа.

Второй случай разреженные данные. В этом случае всего несколько наблюдений (а иногда всего одно) в слабо регламентированное время могут быть получены от одного субъекта. Таким образом, индивидуальные профили не могут быть построены, но данное обстоятельство компенсируется большим количеством наблюдений. В этом случае данные могут быть стратифицированы, и модель может быть построена по обобщенным данным для каждой страты. Подобное популяционное моделирование может быть эффективно при большом количестве данных, когда индивидуальные параметры не представляют особого интереса. Данный подход не учитывает того, что некоторые наблюдения могут быть от одного субъекта и, следовательно, скоррелированы. Эта проблема решается применением нелинейных моделей со смешанными эффектами в которых можно учесть внутри- и меж-индивидуальную вариацию, а также ковариацию и применение сложных вариационно-ковариационных структур, но эта большая тема лежит за границами данной работы. Тем не менее, нельзя сказать, что в практике используется только нелинейное моделирование со смешанными эффектами. В зависимости от требований и обстоятельств могут быть применимы разные подходы к анализу, как простой анализ с фиксированными факторами, так и анализ параметров индивидуальных моделей.

Связка DiffEquations.jl и Optim.jl позволяет наглядно конструировать модель и проводить анализ. При этом возможно построение как индивидуальных, так и простых популяционных моделей.

Полный код:

using DifferentialEquations,  Plots, Optim, LossFunctions, Random, Distributions, DiffEqParamEstim#Functionfunction pkf!(du, u, p, t)  K = p[1]  K= p[2]  K= p[3]  K = p[4]  A = u[1]  C = u[2]  C = u[3]  Vf = p[5]  du[1] = - K  * A  du[2] =   K  * A / Vf   - K * C  - K * C  + K * C  du[3] =   K * C - K * Cend#Start valueu0    = [3.0, 0.0, 0.0]#Parametersp     = [0.8, 0.35, 0.20, 0.15, 0.8]#Time spantspan = (0.0, 50.0)#ODE problemprob  = ODEProblem(pkf!, u0, tspan, p)#Solutionsol   = solve(prob)#Plots#plot(sol)#plot(sol,vars=(1))plot(sol,vars=(2))#plot!(sol,vars=(3))#x = float.(append!(collect(0.25:0.25:3.75), collect(4:2:40)))x = [0.25, 0.5, 1.0, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0,  8.0, 12.0, 24.0, 36.0, 48.0, 60.0]y = hcat(sol.(x)...)'[:,2] .* exp.(rand(Normal(), length(x)) ./ 10)scatter!(x,y)function loss_function(sol, x, y)   tot_loss = 0.0   if any((s.retcode != :Success for s in sol))     tot_loss = Inf   else     tot_loss = sum(value(L2DistLoss(), hcat(sol.(x)...)'[:,2], y))     #tot_loss = sum(value(LogitDistLoss(), hcat(sol.(x)...)'[:,2], y))     #tot_loss = sum(value(HuberLoss(mean(y)), hcat(sol.(x)...)'[:,2], y))   end   tot_lossendcost_function = build_loss_objective(prob, Tsit5(), f ->  loss_function(f, x, y))result = optimize(cost_function, [0.4, 0.1, 0.1, 0.1, 1.0], NelderMead(), Optim.Options(allow_f_increases = true, iterations = 5000, time_limit = 60.0))prob   = ODEProblem(pkf!, u0, tspan, result.minimizer)sol    = solve(prob)plot!(sol,vars=(2))

Что почитать:

• Basic pharmacokinetics
• Useful Pharmacokinetic Equations
• Mathematical Expressions of the Pharmacokinetic andPharmacodynamic Models implementedin the PFIM software
• Introduction to PK/PD modelling with focus on PK and stochastic differentialequations
• Parameter Estimation in Nonlinear Mixed Effect Models Using saemix, an R Implementation of the SAEM Algorithm. Emmanuelle Comets, Audrey Lavenu, Marc Lavielle

скришнот из metal slug 3

2020 год это определенно год странностей. Мой код тоже часто включает в себя некоторые странные ошибки. И в данном посте я хочу показать вам несколько методов отладки кода на языке julia.

Я ни с какой стороны не профессионал в этом деле, да и это справедливо для всего, о чем я пишу в блоге, так что просто имейте это в виду Ну, на самом деле некоторые из вас платят за мою работу, так что технически я могу назвать себя профессиональным блогером, не так ли?

Во всяком случае, давайте не будем отвлекаться на эту мысль. Добро пожаловать в мой блог, если вы новичок, и добро пожаловать обратно в противном случае. Хорошо, что ваш компьютер запрашивает что-то с моего сервера.

Я предполагаю, что у вас есть некоторые базовые знания о Джулии. Следующие посты могут дать вам основы, если вы заинтересованы:

• Множественная диспетчеризация или что в Джулии особенного
• Работа с REPL
• Как работать с пакетами

Кроме того, нужно знание базового синтаксиса.

Пример кода

Прежде чем мы начнем отладку, я хочу продемонстрировать это на некотором коде. Он достаточно короткий, чтобы показать его здесь, и содержит по крайней мере одну ошибку.

В качестве примера возьмем задачку ProjectEuler problem #21. Можете попробовать решить сами. Тут будет начало реализации возможной наивной версии.

Задача заключается в следующем: мы ищем дружественные числа меньше 10 000. Дружественное число определяется как элемент дружественной пары
Пара двух целых чисел (a,b) дружна, если d(a) = b и d(b) = a, где d сумма делителей, так что d(4) = 1+2 = 3.

Дана дружная пара a = 220 и b = 284.
Мы могли бы начать с функции, которая просто берет пару и решает, является ли она дружественной.

function is_amicable(a, b)    sum_divisors(a) == b && sum_divisors(b) == aend

Джулия всегда возвращает выходные данные последнего выполненного выражения в функции. Это означает, что ключевое слово return не обязательно.

Затем нам понадобится функция sum_divisors

function sum_divisors(a)    result = 0    for i = 1:a        if a % i == 0            result += i        end    end    return resultend

которая вызывается так

julia> is_amicable(220, 284)false

Возможно, вы заметили ошибку, но если нет, то, вероятно, лучше не искать ее сейчас. Вместо этого следуйте по пути отладки.

Отладка с помощью Debugger.jl в REPL

Этот пост показывает вам два различных варианта отладки, и первый вариант может быть выполнен в REPL или в вашей IDE, то есть VSCode.

В этом разделе я объясню, как работать с отладчиком на REPL. (Debugger.jl)

julia> ] add Debuggerjulia> using Debugger

Вы можете глянуть в пост про менеджер пакетов, если что-то не ясно.

julia> @enter is_amicable(220, 284)In is_amicable(a, b) at REPL[7]:1 1  function is_amicable(a, b)>2      sum_divisors(a) == b && sum_divisors(b) == a 3  endAbout to run: (sum_divisors)(220)1|debug> 

Я набил @enter is_amicable(220, 284), чтобы получить этот вывод. Кстати, я только что скопировал две функции, которые я определил ранее, в REPL. С другой стороны, Вы можете создать для этого файл amicable.jl и использовать Revise и include (см. REPL and Revise.jl).

В случае файла номера строк, вероятно, более полезны.

Я вернусь через секунду...

julia> using Revisejulia> includet("amicable.jl")julia> using Debuggerjulia> @enter is_amicable(220, 284)In is_amicable(a, b) at /home/ole/Julia/opensources/blog/2020-10-27-basics-debugging/amicable.jl:1 1  function is_amicable(a, b)>2      sum_divisors(a) == b && sum_divisors(b) == a 3  endAbout to run: (sum_divisors)(220)1|debug> 

Готово. Хорошо, теперь как уже упоминалось, в конце мы собираемся запустить sum_divisors(220).

Последняя строка 1|debug> дает нам возможность исследовать дальше, прыгая по коду, в том числе и низкоуровневому, и много чего еще всякого интересного.
Можно посмотреть полный список команд: Debugger.jl commands

Вы также можете ввести ? в режиме отладчика и нажать клавишу enter, чтобы увидеть список команд

Давайте начнем с n шаг к следующей строке.

1|debug> nIn is_amicable(a, b) at /home/ole/Julia/opensources/blog/2020-10-27-basics-debugging/amicable.jl:1 1  function is_amicable(a, b)>2      sum_divisors(a) == b && sum_divisors(b) == a 3  endAbout to run: return false

Значит sum_divisors(220) != 284. Мы, вероятно, хотим перейти к вызову sum_divisors(220).

Мы всегда можем выпрыгнуть из сеанса отладки с помощью q, а затем начать все сначала
Начнем снова с @enter is_amicable(220, 284) и используем s для шага в функцию

1|debug> sIn sum_divisors(a) at /home/ole/Julia/opensources/blog/2020-10-27-basics-debugging/amicable.jl:5  5  function sum_divisors(a)> 6      result = 0  7      for i = 1:a  8          if a % i == 0  9              result += i 10          endAbout to run: 01|debug> 

А дальше продолжаем использовать n, но вы, вероятно, можете себе представить, что это займет некоторое время.

Какие еще инструменты у нас есть, чтобы проверить, что происходит?

Некоторые из вас могут подумать: Хорошо, мы должны, по крайней мере, выяснить, что мы возвращаем, и мы можем просто вызвать sum_divisors(220). Это, вероятно, правильно, но не показывает возможности отладчика. Давайте представим, что мы имеем доступ только к режиму отладчика и не можем просто вызвать функцию.

В общем, этот способ узнавания нового, скрывая то, что мы уже знаем, довольно эффективен.

Я думаю, что пришло время, чтобы представить силу точек останова.

Вместо того, чтобы ползти по программе строка за строкой, часто разумно перейти к определенной точке, запустив код до тех пор, пока эта точка не будет достигнута.

Вы можете сделать это с помощью bp add, а затем указать файл, номер строки и возможное условие. Вы можете увидеть все параметры с помощью ? в режиме отладки.

В наших интересах будет поставить bp add 12. После этого мы можем использовать команду c, которая расшифровывается как continue (до точки останова).

1|debug> cHit breakpoint:In sum_divisors(a) at /home/ole/Julia/opensources/blog/2020-10-27-basics-debugging/amicable.jl:5  8          if a % i == 0  9              result += i 10          end 11      end>12      return result 13  endAbout to run: return 504

Итак, теперь мы знаем, что оно возвращает 504 вместо 284. Теперь мы можем использовать `, чтобы перейти в режим Джулии. (Я знаю, что это вроде как запрещено нашими правилами, но время от времени это имеет смысл, и мы видим, что мы находимся в 1|julia>, а не в julia>, так что я думаю, что все в порядке...)

504-284 это не самый сложный расчет, но мы можем использовать julia, чтобы сделать это за нас, не выходя полностью из режима отладки, используя:

1|debug> `1|julia> 504-284220

Похоже, мы нашли ошибку. Мы добавляем само число к результату, но оно на самом деле не считается за множитель.

А это значит, что мы можем сделать:

function sum_divisors(a)    result = 0    #for i = 1:a    for i = 1:a-1        if a % i == 0            result += i        end    end    return resultend

чтобы избежать эту проблему.

Да, я знаю, что мы можем избежать большего количества чисел, чтобы быть быстрее

Мы можем выйти из режима вычислений с помощью backspace, а затем q, чтобы выйти из режима отладки. Запускаем

julia> is_amicable(220, 284)true

и видим, что мы выковыряли этот баг.

Давайте запустим его в последний раз в сеансе отладки и посмотрим на переменные. Снова перейдем к точке останова c и запустим

1|debug> w add i1] i: 2191|debug> w add a1] i: 2192] a: 220

Теперь мы видим переменные. Если мы снова нажмем c, то снова перейдем к точке разрыва (для очередного вычисления sum_divisors(284) == 220).
Мы можем снова использовать букву w, чтобы увидеть список переменных в области видимости:

1|debug> w1] i: 2832] a: 284

Есть еще несколько способов поиграть, например, шагнуть в код, показать низкоуровневый код и многое другое. Этого должно быть достаточно для знакомства.
В следующем разделе я хочу привести вам тот же пример с помощью редактора кода visual studio с расширением julialang.

Использование VSCode

Я думаю, что большинство разработчиков Julia используют VSCode IDE и, по крайней мере, иногда, vim, emacs или еще что-то такое неудобное Ладно, это, наверное, просто слишком неудобно для меня

Определенно пришло время переключиться на VSCode с Atom/Juno, поскольку расширение Julia теперь разработано для VSCode вместо Atom.

Поскольку это IDE, имеет смысл иметь более визуальный отладчик, чем тот, который описан в предыдущем разделе.

Он довольно прост в навигации, и по умолчанию вы получаете больше выходных данных.

Чтобы начать сеанс отладки, вы нажимаете на кнопку с ошибкой и воспроизводите знак слева, пока у вас открыт файл julia.
Я добавил последнюю строку is_amicable(220, 284), так как VSCode просто запускает программу.

Вы можете добавить точку останова, щелкнув слева от номера каждой строки.

Я сделал снимок экрана после того, как сделал эти шаги, и последним шагом было нажатие на кнопку отладки.

Через несколько секунд сеанс отладки приостанавливается по мере достижения точки останова. С левой стороны можно увидеть локальные переменные в этой позиции. Это этап после того, как я исправил ошибку, так что вы можете видеть, что возвращается правильный результат "284". Однако вы также получаете значение для a и i.

Короче, все то же, что мы делали раньше с нашими переменными, но там нам пришлось вручную добавлять их.

Теперь мы также можем вручную добавлять выражения для наблюдения. Это можно сделать в части Watch ниже Variables, которая находится за пределами скриншота. Довольно приятно иметь возможность добавлять точки останова одним щелчком мыши, а также иметь локальные переменные, показанные слева по умолчанию.

Вы можете спросить себя: Ну, на самом деле это не два способа отладки, не так ли? Это примерно одно и то же, только с другим графическим интерфейсом.

Это правда! Вот почему я сейчас перехожу к следующему разделу поста

Infiltrator.jl для скорости

Существует одна огромная проблема с отладчиком Julia, которая решается по-разному различными пакетами. Проблема в том, что отладчик работает в интерпретируемом режиме, что делает его очень медленным. Если вы отлаживали код C++, то знаете, что отладчик там тоже работает медленнее, чем выполнение, но для Джулии это, на мой взгляд, огромная проблема.

Можно перейти в скомпилированный режим с отладчиком, но это экспериментально, и, по крайней мере, для меня он никогда не останавливался на точке останова.

Некоторые другие пакеты пытаются исправить эту проблему, делая какую-то причудливую магию, но я лично большой поклонник Infiltrator.jl. Правда, некоторое время не было никаких обновлений, и у меня есть некоторые проблемы с ним, но мне нравится сама идея.

Это также один из тех проектов с менее чем 100 звездами. Я хочу подтолкнуть его к этой вехе, так что если вам нравится то, что вы видите в этом разделе, пожалуйста, щелкните им звездочку.

Infiltrator.jl идет совершенно другим путем. Прежде всего, вам нужно немного изменить свой код. Он предоставляет макрос @infiltrate. О боже, как я люблю это название

Макрос примерно такой же, как и предыдущая точка останова. Всякий раз, когда достигается нужная строка, открывается новый вид режима REPL. Это немного усложняет переключение между режимом отладки и обычным режимом запуска, так как вам нужно добавить или удалить макросы @infiltrate, но я думаю, что это нормально.

Я снова продемонстрирую это на примере разобранном выше. Подобного рода использование было в debugging ConstraintSolver.jl.

Я скопировал код сверху и просто добавил using Infiltrator и @infiltrate.

using Infiltratorfunction is_amicable(a, b)    sum_divisors(a) == b && sum_divisors(b) == aendfunction sum_divisors(a)    result = 0    for i = 1:a-1        if a % i == 0            result += i        end    end    @infiltrate    return resultendis_amicable(220, 284)

При запуске кода с include("amicable.jl") получаем:

Hit `@infiltrate` in sum_divisors(::Int64) at amicable.jl:14:debug> 

Это означает, что мы знаем, какая точка останова была достигнута, и видим тип переменной, которую мы назвали sum_divisors. Однако в отличие от Debugger.jl мы не видим кода.

Вы можете снова увидеть раздел справки с помощью ?

debug> ?  Code entered is evaluated in the current function's module. Note that you cannot change local  variables.  The following commands are special cased:    - `@trace`: Print the current stack trace.    - `@locals`: Print local variables.    - `@stop`: Stop infiltrating at this `@infiltrate` spot.  Exit this REPL mode with `Ctrl-D`, and clear the effect of `@stop` with `Infiltrator.clear_stop()`.

Существует не так уж много команд, поэтому мы можем просто попробовать их одну за другой:

debug> @trace[1] sum_divisors(::Int64) at amicable.jl:14[2] is_amicable(::Int64, ::Int64) at amicable.jl:4[3] top-level scope at amicable.jl:18[4] include(::String) at client.jl:457

Таким образом, мы пришли из is_amicable и можем видеть типы, а также имя файла и номер строки, что полезно при использовании multiple dispatch.

debug> @locals- result::Int64 = 284- a::Int64 = 220

мы можем видеть локальные переменные, которые похожи на те, которые мы видели в представлении переменных VSCode.

Кроме того, мы можем просто вычислять выражения прям в этом режиме. Для Infiltrator.jl нет необходимости использовать `, чтобы переключиться на вычисления.

debug> a == 220true

Вы можете использовать @stop, чтобы больше не останавливаться на этой вехе, и Infiltrator.clear_stop(), чтобы очистить эти остановки.

Давайте не будем использовать @stop сейчас, а вместо этого перейдем к следующей точке @infiltrate с помощью CTRL-D:

Hit `@infiltrate` in sum_divisors(::Int64) at amicable.jl:14:debug> 

таким образом, мы находимся в той же точке останова, но со вторым вызовом. К сожалению, невозможно использовать клавишу со стрелкой вверх, чтобы перейти через историю команд, которые мы использовали, так что нам нужно снова ввести @locals, если мы хотим их увидеть.

Я открыл такую тему и попытался разрешить ее сам, но мне это не удалось. Было бы здорово, если бы это когда-нибудь было реализовано, потому что я думаю, что было бы очень полезно отлаживать быстрее таким образом.

Давайте рассмотрим сравнение двух различных способов в следующем разделе.

Выводы

Мы посмотрели на Debugger. jl, который дает вам всю информацию, которая может понадобиться в вашем REPL.
Поэтому он не зависит от редактора.

Следующий инструмент, который я упомянул, был дебагер в VSCode который является в основном просто графическим интерфейсом для Debugger.jl. Вероятно, его удобнее использовать людям, которые любят работать с IDE. Хотя в Debugger.jl могут быть некоторые опции, которые недоступны в графическом интерфейсе, как это часто бывает.

Оба этих инструмента дают то преимущество, что вы можете шаг за шагом переходить через свой код и исследовать все, что захотите. Вы можете взглянуть на низкоуровневый код (по крайней мере, в Debugger.jl). Это, вероятно, то, что каждый ожидает сделать с отладчиком. Проблема просто в том, что он слишком медленный во многих случаях использования, например, когда вы хотите отладить свой собственный пакет с 1000 строками кода.

В таком случае Infiltrator.jl это путь, по крайней мере для меня, и пока скомпилированный режим Debugger.jl работает недостаточно хорошо. У него есть и другие недостатки, так как не бывает чтоб все и сразу, но я думаю, что он часто превосходит использование println, поскольку можно распечатать все, что в данный момент интересует в данной точке останова, и увидеть все локальные переменные за один раз.

Спасибо за то, что дочитали и особая благодарность моим 10 покровителям!

Я буду держать вас в курсе Twitter OpenSourcES.

источник изображения

Возможно, вы видели предыдущий пост, где были предоставлены визуализации первых 1000 цифр и . Он возник в результате небольшого спора о том, лучше ли , чем . По этому поводу идут бесконечные дебаты, и я подумал, что могу пошутить по этому поводу. В этом посте я хочу показать, как создать визуализации, и надеюсь, что вы захотите попробовать удивительный пакет Luxor.jl после прочтения. Вчера я начал читать туториал, и это потрясающе! В прошлый раз визуализация делалась на Javascript, и я подумал, что этот аккуратный маленький проект сойдет, чтобы начать изучать Луксор. Как уже упоминалось в let me be your mentor: я думаю, что очень важно иметь такие маленькие проекты, чтобы освоить новый инструмент.

Основная идея

Я хотел воссоздать визуализацию, которую видел в Numberphile от Мартина Крживинского.

Там был круг (который, вполне ассоциируется и с и с ) разделенный на 10 сегментов, по одному для каждой цифры. Цифры нашего иррационального числа представляются кривыми внутри этого круга, так что 3.1415 (я начинаю с 14) это кривая от сегмента 1 до сегмента 4, а затем обратно к 1, потом до 5 и так далее. Каждый раз мы перемещаемся немного по часовой стрелке в сегменте так, что 14 создает различные кривые (в зависимости от текущего положения, в котором мы находимся).

Потом надобавляем всякие фичи. Мы должны начать чувствовать себя комфортно с Луксором. Важно: не надо искать математическую интерпретацию это просто небольшой проект визуализации ;)

Я знаю, вам интересно, как должен выглядеть конечный результат:

Начинаем

using Luxorfunction vis()    @png begin        sethue("black")        circle(O, 50, :stroke)        setdash("dot")        circle(O, 70, :stroke)        sethue("darkblue")        circle(O, 10, :fill)    end 500 200 "./start.png"end

вызываем vis() и создаем файл start.png который будет выглядеть как-то так:

Давайте быстренько пройдемся по командам:

@png beginend width height "filename.png"

просто хороший макрос. :)

sethue задает цвет и принимает либо строку, как показано выше или цвет пакета из Colors. Он устанавливает цвет для следующих команд рисования до тех пор, пока вы не выберете другой. То же самое верно и при установке ширины линии с помощью setline, или при установке размера шрифта, или при других общих настройках.

Команды рисования, такие как circle, обычно принимают некоторые параметры и заканчиваются параметром действия, таким как :stroke или :fill.

О это буква "О", а не число "0". :) Она представляет собой начало координат и является краткой формой для Point(0, 0). В Луксоре начало находится в центре полотна. В качестве второго параметра должен быть задан радиус.

Давайте сначала нарисуем внешний круг и добавим цифры:

radius = 100@png begin    background("black")    sethue("white")    circle(O, radius, :stroke)    for i in 0:9         = 2*0.1*i+0.1*        mid = Point(            radius*sin(),            -radius*cos(),        )        label(string(i), :N, mid)    endend 700 300 "./first_step.png"

Первая часть должна быть достаточно простой.

 = 2*0.1*i+0.1*

возможно, это не идеально написано (кроме того, я мог бы использовать :D). 2*0.1*i начинает с северного положения, а затем для следующего i происходит перемещение на . Я добавляю "0.1 ", потому что хочу переходить к середине каждого сегмента. Может быть, следует написать 0.5/10*2. Затем мы просто поворачиваем наш холст и двигаясь чуть выше радиуса, рисуем метки. На самом деле такое можно проделать в Luxor, используя rotate и translate. Но я решил сделать вручную, так как мне все равно это пригодится позже. В общем формула такова:

Такое преобразование поворачивает плоскость на и производит трансляцию на x,y. Поскольку я перевожу только на y, мне не нужно первое тождество. Помните, что y увеличивается, когда идет вниз.

В настоящее время есть две проблемы:

• на самом деле нам не нужен круг, нам нужны дуги (сегменты) для каждой цифры
• подписи не читаются

Команда label принимает три значения: текст, вращение и положение, где вращение может быть записано как :N,: E,: S,: W для севера, востока, юга, запада или как угол (в радианах). :N есть . Поэтому мы хотим начать с , а потом добавлять текущий угол поворота. Кроме того, смещение было бы здорово, если бы оно не доставало непосредственно до окружности или не подходило слишком близко к ней. Здесь мы могли бы увеличить радиус или использовать ;offset в команде label.

Для первой задачи нам нужна функция arc2r, которая принимает три аргумента
c1, p1, p2 + действие: c1 это центр окружности, а p1 и p2 точки на окружности, между которыми должен быть показан сегмент. По умолчанию выбрано направление по часовой стрелке.

Мы определяем следующую функцию, чтобы получить и соответствующую точку более простым способом:

function get_coord(val, radius)     = 2*0.1*val    return Point(        radius*sin(),        -radius*cos(),    )end

а потом:

background("black")for i in 0:9    from = get_coord(i, radius)    to = get_coord(i+1, radius)    randomhue()     = 2*0.1*i+0.1*    mid = Point(        radius*sin(),        -radius*cos(),    )    label(string(i), -/2+, mid; offset=15)    move(from)    arc2r(O, from, to, :stroke)end

Я использовал randomhue, чтобы получить случайный цвет. Мы исправим это в следующий раз :)
Также я переставлял порядок Label и arc2r и поставил move, так как в противном случае линии рисуются от метки дуги. Это происходит потому, что arc продолжает текущий путь.

Выглядит намного лучше! Давайте возьмем несколько хороших цветов из Colorschemes.jl.

Я использовал схему rainbow, начиная с 7-го цвета :D. Вы, возможно, захотите испытать другие цветовые схемы, так как здесь цвета не так легко различить, но мне все равно почему-то нравится именно она.

using ColorSchemescolors = ColorSchemes.rainbow[7:end]

и затем

sethue(colors[i+1])

помните, что индексация массивов в Julia начинается с единицы.

Каковы следующие шаги?

• Добавление строк
• Рефакторинг кода
• Оживление процесса
• Добавление точек
• Добавление гистограммы сверху

Я думаю, что визуально привлекательно иметь круг посередине, где мы можем добавить символ (или ) позже.
Поэтому мы не можем провести прямые линии от одного сегмента к другому. Для этого я использую квадратичные кривые Безье.

Давайте сначала получим цифры числа Пи:

max_digits = 10digits = setprecision(BigFloat, Int(ceil(log2(10) * max_digits+10))) do    return parse.(Int, collect(string(BigFloat(pi))[3:max_digits+2]))end

это дает нам первые 10 цифр после десятичной точки числа Пи. Для этого мне нужно установить точность BigFloat. Довольно интересно, что пи не является жестко закодированной константой в Джулии. Оно вычислено таким образом, что я в принципе могу получить любую точность, какую захочу. Точность должна быть задана в количестве битов, так что необходимо выполнить небольшое вычисление. Я добавил +10 в конце, чтобы быть уверенным :D

Чтобы нарисовать квадратичную кривую Безье, нам нужны три точки. Начало, конец и контрольная точка. В качестве контрольной точки я выбираю точку на внутреннем круге, который просто также разделен на десять сегментов, и выбираю сегмент, который находится посередине между текущей цифрой from_val и следующей цифрой to_val.

Я должен уточнить, что я имею в виду под серединой: средняя точка между 0 и 4 должна быть 2, но между 8 и 0 она должна быть 9. Она определяется кратчайшим путем от одного сегмента к другому, а потом берется середина.

Кроме того, у меня на самом деле нет 10 дискретных сегментов, это просто для понимания. Я могу использовать среднюю точку 1,23 или что-то в этом роде. Это используется, потому что мы меняем нашу начальную и конечную позиции на основе текущей позиции, которую мы находимся в нашем массиве цифр.

Я надеюсь, что все станет яснее, ели взглянуть на код:

small_radius = 70for i in 1:max_digits-1    from_val = digits[i]    to_val = digits[i+1]    sethue(colors[from_val+1])    f = from_val+(i-1)/max_digits    t = to_val+i/max_digits    from = get_coord(f, radius)    to = get_coord(t, radius)    # get the correct mid point for example for 0-9 it should be 9.5 and not 4.5    mid_val = (f+t)/2    mid_control = get_coord(mid_val, small_radius)    if abs(f-t) >= 5        mid_control = get_coord(mid_val+5, small_radius)    end    pts = Point[from, mid_control, mid_control, to]    bezpath = BezierPathSegment(pts...)    drawbezierpath(bezpath, :stroke, close=false)end

Думаю, уже выглядит достаточно хорошо. Цвета линий подгоняются под цвета из под цифр. Итак, в какой-то момент мы переходим от 9 к 2. Вместо этого я хотел бы посмотреть, куда мы идем и откуда идем. Это можно сделать с помощью blend и setblend. Это линейная смена цвета "от" и "до", так что на самом деле не по кривой, но я думаю, что она достаточно хороша.

setblend(blend(from, to, colors[to_val+1], colors[from_val+1]))

Это похоже на sethue поэтому нам нужно задать его в какой-то момент, прежде чем мы вызовем drawbezierpath.

Давайте добавим еще несколько цифр и немного уменьшим ширину линии: setline(0.1)

Ладно я думаю что внутренний радиус немного велик:

small_radius = 40

Затем мы можем добавить в середине, прежде чем немного очистить код, чтобы создать нашу первую анимацию.

Luxor.jl не поддерживает латексные стринги LaTeXStrings.jl это облом, но мы можем использовать UnicodeFun.jl.

using UnicodeFuncenter_text = to_latex("\\pi")

и промеж циклов ставим:

sethue("white")fontsize(60)text(center_text, Point(-2, 0), valign=:middle, halign=:center)

Мне кажется Point(-2, 0) более центральная, чем Point(0, 0) или O.

Анимация

Я хотел бы получить gif из конвейера визуализации таким образом, чтобы в каждом кадре добавлялась новая линия.

В Луксоре это можно сделать с помощью функции animate, которая берет несколько сцен и их номера кадров. Это также обеспечит немного большую структуру кода.

У нас может быть сцена для устойчивого фона и одна для линий.

Прежде чем мы напишем функцию, давайте определим очень короткую анимацию, чтобы увидеть, как это делается.

function draw_background(scene, framenumber)    background("black")endfunction circ(scene, framenumber)    setdash("dot")    sethue("white")    translate(-200, 0)    @layer begin         translate(framenumber*2, 0)        circle(O, 50, :fill)    endendfunction anim()    anim = Movie(600, 200, "test")    animate(anim, [        Scene(anim, draw_background, 0:200),        Scene(anim, circ, 0:200),    ],    creategif = true,    pathname = "./test.gif"    )end

Сначала мы создаем Movie с width, height и name.
Затем мы вызываем animate с помощью созданного Movie и списка scenes, а затем функции и диапазон кадров, начинающихся с 0.

Происходит вызов draw_background(сцена, 0) и circ(scene, 0) для первого кадра. Сцена может содержать некоторые аргументы, которые мы будем использовать для нашей анимации. Остальное в основном так же, как и раньше, просто мы можем, конечно, использовать переменную framenumber.

Теперь я разделю все это дело на функции и определю переменные, такие как цифры, которые мы хотим визуализировать, чтобы нам было легче визуализировать или другие вещи.

using Luxor, ColorSchemesusing UnicodeFunfunction get_coord(val, radius)     = 2*0.1*val    return Point(        radius*sin(),        -radius*cos(),    )endfunction draw_background(scene, framenumber)    background("black")    radius = scene.opts[:radius]    colors = scene.opts[:colors]    center_text = scene.opts[:center_text]    for i in 0:9        from = get_coord(i, radius)        to = get_coord(i+1, radius)        sethue(colors[i+1])         = 2*0.1*i+0.1*        mid = Point(            radius*sin(),            -radius*cos(),        )        label(string(i), -/2+, mid; offset=15)        move(from)        arc2r(O, from, to, :stroke)    end    sethue("white")    fontsize(60)    text(center_text, Point(-2, 0), valign=:middle, halign=:center)endfunction dig_line(scene, framenumber)    radius = scene.opts[:radius]    colors = scene.opts[:colors]    center_text = scene.opts[:center_text]    bezier_radius = scene.opts[:bezier_radius]    max_digits = scene.opts[:max_digits]    digits = scene.opts[:digits]    setline(0.1)    for i in 1:min(framenumber, max_digits-1)        from_val = digits[i]        to_val = digits[i+1]        f = from_val+(i-1)/max_digits        t = to_val+i/max_digits        from = get_coord(f, radius)        to = get_coord(t, radius)        # get the correct mid point for example for 0-9 it should be 9.5 and not 4.5        mid_val = (f+t)/2        mid_control = get_coord(mid_val, bezier_radius)        if abs(f-t) >= 5            mid_control = get_coord(mid_val+5, bezier_radius)        end        pts = Point[from, mid_control, mid_control, to]        bezpath = BezierPathSegment(pts...)        # reverse the color to see where it is going        setblend(blend(from, to, colors[to_val+1], colors[from_val+1]))        drawbezierpath(bezpath, :stroke, close=false)    endendfunction anim()    anim = Movie(700, 300, "test")    radius = 100    bezier_radius = 40    colors = ColorSchemes.rainbow[7:end]    max_digits = 1000    center_text = to_latex("\\pi")    digits_arr = setprecision(BigFloat, Int(ceil(log2(10) * max_digits+10))) do        return parse.(Int, collect(string(BigFloat(pi))[3:max_digits+2]))    end    args = Dict(:radius => radius,        :bezier_radius => bezier_radius,        :colors => colors, :max_digits => max_digits,        :digits => digits_arr, :center_text => center_text    )    animate(anim, [        Scene(anim, draw_background, 0:max_digits+50, optarg=args),        Scene(anim, dig_line, 0:max_digits+50, optarg=args),    ],    creategif = true,    pathname = "./pi_first.gif"    )end

Единственное, что я еще не объяснил, это optarg в функции Scene и получение его с помощью radius = scene.opts[:radius].

Мы как бы потеряли возможность создавать простые образы. Поэтому я создал структуру

struct PNGScene    opts::Dict{Symbol, Any}end

и использую некоторые аргументы в функции anim, которую я переименую в viz :D

Тогда я могу использовать что-то вроде:

scene = PNGScene(args)@png begin    draw_background(scene, max_digits)    dig_line(scene, max_digits)end 700 300 "./$fname.png"

Не волнуйтесь, в конце есть репка, где вы можете увидеть весь код целиком. Просто немного сложно описать здесь изменения.

Может, мне стоило снять видео? :D

Добавление точки Фейнмана

Мы визуализировали соединение цифр с цифрами с помощью кривых, но если бы у нас встретилось что-то вроде 555 в цифрах, мы бы видели только линию, идущую в направлении центра и обратно (или, может быть, мы видим две в зависимости от наших максимальных цифр, разрешения и т. д.)

Вместо этого мы можем показать дополнительную точку всякий раз, когда это происходит. Это можно получить благодаря аргументу функции show_dots, что вы можете найти в моем коде ;)

Я только что проверил длину последовательности, и когда она больше 1, я рисую круг, где это происходит, и цвет это цифра после этой последовательности. Большой круг в сегменте 9 это так называемая точка Фейнмана, где цифра 9 появляется 6 раз в позиции 762.

Добавление гистограмм

Последняя вещь в моем списке получить гистограмму на каждом сегменте, чтобы показать, случаются ли некоторые комбинации пар чаще, чем другие.

Для этого я использую функцию poly с четырьмя точками. В идеале, она должна быть ограничена двумя дугами, а не двумя линиями, но я оставляю это читателю :)

Тау

Да, можно было бы в принципе сгенерировать случайное число с 1000 цифрами и получить аналогичный результат...

Простое число

В двух словах: использование нашей функции для визуализации большилства элементов не так разумно, но так или иначе может получится что-то интересное.

При этом в качестве числовой последовательности используются последние цифры простых чисел. Я визуализировал простые числа меньше 100 000. Честно говоря, соединительные линии немного бесполезны, так как большую часть времени (если мы игнорируем первые несколько простых чисел: все время) возможны только четыре цифры. Это создает своего рода беспорядок в середине.

Тем не менее, гистограммы становятся все интереснее, я думаю:

Это ясно показывает, что не все пары одинаково вероятны. Особенно, если у нас есть простое число с последней цифрой x, то всегда менее вероятно, что последняя цифра также заканчивается на x по сравнению с одним из трех других вариантов.

Давайте сосредоточимся на гистограммах и визуализируем простые числа под 10 000 000:

Узор сохраняется.

Код

Окай, тут у нас репка

Я хотел бы создать что-то вроде штучек, из 3b1b.

По крайней мере, небольшие простые версии с некоторыми удобными функциями визуализации :)

Спасибо за чтение и особая благодарность моим 10 покровителям!

Я буду держать вас в курсе событий на Twitter OpenSourcES и на более личном:
Twitter Wikunia_de

Прошло уже достаточно времени с релиза Javis v0.2, что обсуждалось в соответствующем посте. Там я дал представление о потенциальном будущем этого графического пакета. Мы наконец-то выпустили v0.3, и будущее стало стандартом по умолчанию.

Просто перечислять все изменения, которые мы ввели, вероятно, было бы довольно скучно, так что лучше я создам аккуратную анимацию, а по пути буду объяснить некоторые удивительные штуки про ряды Фурье.

Я думаю, что довольно много людей слышали о рядах Фурье, но не совсем понимали их. И я совершенно уверен, что все еще нахожусь в их числе. Не знаю, будет ли целесообразно использовать этот пост в качестве вступления для тех, кто вообще не в теме, но я постараюсь сделать все возможное. К слову, завидую я этим людям, потому что ряды Фурье это нечто чрезвычайно увлекательное.

Этот блог в основном предназначен для читателей, которые интересуются Javis и Julia, поэтому я начну с изменений, которые появились в v0.3.

Изменения в v0.3

Наконец-то у нас есть версия способная транспонировать матрицы, или я предпочитаю говорить об общих возможностях морфинга, поскольку транспонирование работает только иногда.
Надо будет по этой теме подробней отписаться. (Ну, скажем, в ноябре надеюсь)

Что-то, к сожалению, не попало в v0.3, как например объединение объектов в слои, но работа идет полным ходом. Есть уже концептуальный план, осталось найти время, чтобы все закодить, попутно внося возможные изменения.
Еще было обещано поменять имена некоторых ключевых функций, так что надеюсь, что вы либо не использовали Javis v0.2, либо что вы адаптируетесь и довольны изменениями, которые мы провели.

Давайте просто перечислим некоторые из них:

Action и SubAction теперь называются соответственно Object и Action.
Это связано с тем, что Action раньше больше описывал сам объект и мог фактически перемещать его. Мы решили убрать такую возможность чтобы упростить кодовую базу и сделать различие гораздо более четким.

Далее мы упростили BackgroundAction до Background. Больше никаких комментариев по этому поводу

Наконец, мы удалили Translation, Rotation и Scaling и ввели anim_translate, anim_rotate/anim_rotate_around, а также anim_scale, потому что работа с функциями вместо структур кажется более естественной, и название указывает, что мы попутно делаем анимацию, а не статический перевод.

Последнее важное изменение состояло в том, чтобы избавиться от функции javis и странного дизайна, вызывавшего постоянный беспорядок с этими странными массивами.

Теперь просто определяем видео, объекты и действия строка за строкой и, в конце концов, визуализируем видео:

using Javisfunction ground(args...)    background("black")    sethue("white")endvideo = Video(800, 400)Background(1:100, ground)ball = Object((args...)->circle(O, 100, :fill), Point(-500, 0))rolling = Action(anim_translate(Point(1000, 0)))act!(ball, rolling)render(video; pathname="rolling_ball.gif")

Думаю, код предельно прозрачен, кроме, может быть, нескольких причуд ...

Функция circle принимает центральную точку и радиус, а также действие (здесь :fill). Почему я определяю центр в центре? Я мог бы определить его в Point(-500, 0), но мне нравится думать о том, что объекты центрируются, а затем, перемещаются вокруг оси. Мне кажется, что здесь есть преимущество: объект легко и правильно масштабируется, а также это облегчает вращение.

На самом деле невозможно увидеть вращение круга (в настоящее время мы его не вращаем). Как насчет того, чтобы добавить что-то к кругу, чтобы мы могли различать вращается ли он?

using Javisfunction ground(args...)    background("black")    sethue("white")endfunction my_circle(args...)    circle(O, 100, :fill)    sethue("black")    line(O, Point(100, 0), :stroke)endvideo = Video(800, 400)Background(1:100, ground)ball = Object(my_circle, Point(-500, 0))translating = Action(anim_translate(Point(1000, 0)))rotating = Action(anim_rotate(0.0, 2*2))act!(ball, [translating, rotating])render(video; pathname="rolling_ball_2.gif")

Таки катится. Ладно, возможно, я немного отвлекся, но позже нам это пригодится.

Думаю, в v0.2 было бы довольно трудно сделать линию внутри круга с длиной зависящей от времени. На самом деле, это было возможно с v0.2.2, но я еще не писал об этом, так что это время перемен

using Javisfunction ground(args...)    background("black")    sethue("white")endfunction my_circle(args...; line_length=0)    circle(O, 100, :fill)    sethue("black")    setline(2)    line(O, Point(line_length, 0), :stroke)endvideo = Video(800, 400)Background(1:100, ground)ball = Object(my_circle, Point(-500, 0))translating = Action(anim_translate(Point(1000, 0)))rotating = Action(anim_rotate(0.0, 2*2))changing_len = Action(change(:line_length, 0 => 100))act!(ball, [translating, rotating, changing_len])render(video; pathname="rolling_ball_3.gif")

Теперь, наконец, можно просто изменять значения непосредственно анимированным способом. Вот мы вроде как и собрали строительные блоки для Фурье

Ряды Фурье

Простыми словами, мы хотим из кучи простых штук воссоздать все что угодно. Например, можно комбинировать и суммировать синусоидальные волны для аппроксимации любой функции.

Как насчет ступенчатой функции?

using Javis, Colorsfunction ground(args...)    background("black")    sethue("white")    translate(-args[1].width/2+50, 0)    scale(700, 150)endfunction wave(xs, ys, opacity, color)    setline(1.5)    setopacity(opacity)    sethue(color)    points = Vector{Point}(undef, length(xs))    for (x, y, i) in zip(xs, ys, 1:length(xs))        points[i] = Point(x, y)    end    prettypoly(points, :stroke, ()->())endfunction term(x, i)    return 4/ * sin.(2*(2i-1)*x)/(2i-1)endfunction sum_term(x, k)    k == 0 && return zeros(length(x))    return sum(term(x, i) for i in 1:k)endnframes = 300frames_per_wave = 40video = Video(800, 600)Background(1:nframes, ground)x = 0.0:0.001:1.0k = 6colors = [RGB(0.0, 1.0, 0.4), RGB(0, 1.0, 1.0), RGB(1.0, 0.25, 0.25), RGB(1.0, 1.0, 0.0), RGB(1.0, 0.5, 1.0), RGB(0.75, 0.75, 1.0)]waves = Object[]for i = 1:k    frames = frames_per_wave*(i-1)+1:nframes    push!(waves, Object(frames, (args...; y=term(x,i)) -> wave(x, y, 0.5, colors[i])))    act!(waves[end], Action(5:frames_per_wave,                    change(:y, term(x,i) => sum_term(x, i))))endsum_wave = Object(1:nframes, (args...; y=zeros(length(x)))->wave(x, y, 1.0, "white"))for i = 1:k    act!(sum_wave, Action(frames_per_wave*(i-1)+1:frames_per_wave*i,                         change(:y, sum_term(x, i-1) => sum_term(x, i))))endrender(video; pathname="images/fourier_1D.gif")

Мы, в основном, пытаемся объединить простые функции (синусоидальные волны) путем их сложения, чтобы создать квадратную волну.

Важно отметить, что в квадратной функции нет ничего особенного. Подход обобщается на любой пример. Просто нужно насуммировать довольно много функций. Для идеального результата бесконечно много.

А теперь, бегом в 2D

На самом деле, хочется нарисовать что-то в 2D это гораздо интереснее. Кстати, вся эта идея использования Javis для визуализации рядов Фурье пришла от Риккардо Чиоффи, которого я люблю называть первым пользователем Javis. Он сам реализовал идею, без нашего ведома, пока мы не увидели его пост об этом на JuliaLang-дискурсе. Посмотрите, пожалуйста!

Имплементация была на предыдущей версии Javis, поэтому код больше не работает. Я подожду здесь пару минут, пока вы будете наслаждаться его потрясающей визуализацией.

Давайте начнем с некоторых строительных блоков, которые имеют смысл для понимания рядов Фурье в 2D. Я начну с кругов и комплексных чисел.

В общем, я пройдусь по тому материалу, который Грант Сандерсон прекрасно объяснил на своем YouTube-канале 3blue1brown Что такое ряды Фурье?. В конце этого поста есть еще несколько ссылок

Мы уже немного покатали круги, но давайте посмотрим, как работать с ними, используя комплексные числа. Я здесь имею в виду использование формулы

где t наша переменная в диапазоне от 0 до 1, а i мнимая единица.

Таким образом, эта формула дает комплексное число для некоторого t. Мы можем его визуализировать, просто откладывая действительную компоненту нашего комплексного числа на оси x и комплексную компоненту на оси y.

В этой визуализации круг рисуется до определенной точки в каждом кадре и завершает круг с t=1.00

Картинку нужно попридержать в уме до конца поста.

Далее мы начнем с некого многоугольника, который мы хотим нарисовать используя круги, вращающиеся вокруг кругов, поскольку круги это синусоидальные волны двумерного пространства.

Я думаю, что это лучше всего объяснить визуально, поэтому я покажу вам, чего мы пытаемся достичь, а затем покажу, как этого достичь.

Что мы хотим воссоздать

Вот своего рода то, чего мы хотим достичь, но позже мы рассмотрим более сложные формы. Для этого мы будем использовать больше кругов, хоть тогда и труднее понять, что происходит.

using Javis, FFTW, FFTViewsfunction ground(args...)    background("black")    sethue("white")endfunction circ(; r = 10, vec = O, action = :stroke, color = "white")    sethue(color)    circle(O, r, action)    my_arrow(O, vec)    return vecendfunction my_arrow(start_pos, end_pos)    start_pos  end_pos && return end_pos    arrow(        start_pos,        end_pos;        linewidth = distance(start_pos, end_pos) / 100,        arrowheadlength = 7,    )    return end_posendfunction draw_line(    p1 = O,    p2 = O;    color = "white",    action = :stroke,    edge = "solid",    linewidth = 3,)    sethue(color)    setdash(edge)    setline(linewidth)    line(p1, p2, action)endfunction draw_path!(path, pos, color)    sethue(color)    push!(path, pos)    draw_line.(path[2:end], path[1:(end - 1)]; color = color)endfunction get_points(npoints, options)    Drawing()    shape = poly([Point(-200, 0), Point(250, 70), Point(165, -210)]; close=true)    points = [Javis.get_polypoint_at(shape, i / (npoints-1)) for i in 0:(npoints-1)]    return pointsendfunction poly_color(points, action; color=nothing)    color !== nothing && sethue(color)    poly(points, action)endc2p(c::Complex) = Point(real(c), imag(c))remap_idx(i::Int) = (-1)^i * floor(Int, i / 2)remap_inv(n::Int) = 2n * sign(n) - 1 * (n > 0)function animate_fourier(options)    npoints = options.npoints    nplay_frames = options.nplay_frames    nruns = options.nruns    nframes = nplay_frames + options.nend_frames    points = get_points(npoints, options)    npoints = length(points)    println("#points: $npoints")    # optain the fft result and scale    x = [p.x for p in points]    y = [p.y for p in points]    fs = fft(complex.(x, y)) |> FFTView    # normalize the points as fft doesn't normalize    fs ./= npoints    npoints = length(fs)    video = Video(options.width, options.height)    Background(1:nframes, ground)    Object((args...)->poly_color(points, :stroke; color="green"))    circles = Object[]    npoints = 5    for i in 1:npoints        ridx = remap_idx(i)        push!(circles, Object((args...) -> circ(; r = abs(fs[ridx]), vec = c2p(fs[ridx]))))        if i > 1            # translate to the tip of the vector of the previous circle            act!(circles[i], Action(1:1, anim_translate(circles[i - 1])))        end        ridx = remap_idx(i)        act!(circles[i], Action(1:nplay_frames, anim_rotate(0.0, ridx * 2 * nruns)))    end    trace_points = Point[]    Object(1:nframes, (args...) -> draw_path!(trace_points, pos(circles[end]), "red"))    render(video, pathname = options.filename)endfunction main()    gif_options = (        npoints = 1000, # rough number of points for the shape => number of circles        nplay_frames = 400, # number of frames for the animation of fourier        nruns = 2, # how often it's drawn        nend_frames = 0,  # number of frames in the end        width = 800,        height = 500,        shape_scale = 0.8, # scale factor for the logo        tsp_quality_factor = 40,        filename = "images/fourier_tri_5.gif",    )    animate_fourier(gif_options)endmain()

Здесь мы попытаемся нарисовать зеленый треугольник, используя пять кругов. Они аппроксимируют треугольник, как вы можете видеть, красной фигурой. У каждого круга имеется прикрепленный к нему вектор. Начнем с внутренней неподвижной окружности, где вектор фиксирован.

Вектор указывает на центр треугольника, что имеет смысл, поскольку это лучшая точка, чтобы убедиться, что баланс вращающихся вокруг нее кругов не нарушен. Я также вернусь к этому позже в этом длинном-длинном посте.

Как насчет следующего круга в цепочке? Он не движется, но его центр находится там, куда указывал последний круг. Кроме того, вектор внутри теперь вращается. Один полный оборот по часовой стрелке эквивалентен одному проходу по контуру рисунка.

Вы могли бы задаться вопросом, почему это по часовой стрелке, а не против? Это связано с тем, что ось y положительна вниз и отрицательна вверх. Впрочем, это не имеет значения.

А вот следующие три круга двигаются. Вектор внутри каждого последующего в цепочке имеет ту же скорость, что и предыдущий круг, но работает против часовой стрелки. Два других делают два полных круга за один проход по рисунку. Один идет по часовой стрелке, другой против.

Это вроде как все, что нужно. Ну, по крайней мере, это дает визуализацию того, чего мы хотим достичь, и вы уже можете догадаться, что добавление большего количества кругов делает приближение лучше.

Как это посчитать

В этом разделе я хочу немного поговорить о формуле, которую некоторые из вас, возможно, видели раньше, и вычислить параметры пяти кругов, показанных выше.

Мы уже видели раньше, так что k это единственная новая вещь в экспоненциальной части. Из анимации круга выше видно, что k в основном, просто определяет скорость, а знак k направление вращения.

Это означает, что у нас есть невращающийся круг с k=0 и вращающиеся круги с . Единственное, что осталось, это .

определяет радиус круга, например, 2.5 или 0.5. Кроме того, они задают начальное вращение круга. Это показывается с помощью внутреннего вектора. Мы уже немного посмотрели на комплексные числа, и при определении как комплексного числа у нас есть необходимые два компонента, заданные абсолютным значением и углом комплексного числа.

Такое объяснение отлично от того, что Вы можете найти на других сайтах, поскольку есть разные способы выражения данной концепции. Моя версия эквивалентна тому, как Грант Сандерсон объяснил это в своем видео Что такое ряды Фурье?.

Давайте посмотрим на компоненты для n=1:

что может быть упрощено с помощью :

Что произойдет, если мы возьмем среднее значение этого показателя, когда позволим t перейти от 0 к 1 и нарисуем полные круги?

Прежде всего, мы можем вычислить среднее значение суммы, суммируя средние значения.

Если мы усредним положения окружности, то получим центр, который находится в точке (0,0) для всех из них. Поэтому усреднение по всем возможным t между 0 и 1 (или, по крайней мере, хорошее подмножество всех возможностей хорошо, скажем, несколько сотен значений вместо бесконечного множества)
мы получаем . Чудненько, да?

Давайте повторим это еще раз: мы усредняем по слагаемым, и поскольку все слагаемые, кроме , являются кругами, они усредняются до и не имеют отношения к делу.

Это была также интуиция, которую мы получили из результата, показанного с пятью кругами, отчаянно пытающимися приблизиться к треугольнику.

Ну подождите, а как мы получим и все остальные константы ?

На самом деле мы делаем точно то же самое и просто умножаем все уравнение на , чтобы получить сдвинутую версию. Помните, что такое, что мы можем отменить показатель степени с помощью этого умножения. Умножая на , мы вычеркиваем компонент, который у нас есть для . Мы можем сделать то же самое для всех констант.

Давайте сделаем хорошую визуализацию для этого:

Немного сложно для восприятия, так что остановимся поподробней:

1. Начнем с зеленого треугольника
2. Выбираем 100 точек на этом треугольнике и усредняем его
3. Создаем окружность с центром в начале координат и указываем на среднее значение с помощью вектора
4. Умножаем треугольник на с и t зависящим от положения, в котором мы находимся на треугольнике.
5. Создаем новые 100 точек и снова усредняем их
6. Создаем новый круг
7. Перемещаем вновь созданный круг на вершину предыдущего круга

И так повторяем, чтобы получить все остальные константы и, следовательно, круги, которые нам нужно добавить.

Работать сподручней с пакетами

Вместо того чтобы умножать и усреднять собственноручно, мы можем использовать для этого пакет. Я надеюсь, что смогу убедить всех, кто не был убежден раньше, что ряды Фурье это прекрасно, и конечно, с ними связано гораздо больше крутых штук, чем я объяснил здесь или понял сам

Есть два хороших пакета, которые я использовал в своей анимации треугольника Фурье показанном раньше:

• FFTW
• FFTViews

Первый из них выполняет тяжелые вычисления, и делает это быстро, поэтому первая буква "F" в названии означает "быстрый".

С другой стороны, FFTViews просто облегчает работу. Как вы уже видели, мы имеем дело с отрицательными и положительными значениями для k.

FFTW создает вектор со всеми этими константами , но FFTViews позволяет индексировать эту информацию более естественно, например так c[-3], даже несмотря на то, что Julia это язык, с индексацией начинающейся с единицы.

А потом вводится вектор точек как вектор комплексных чисел. Дополнительный шаг, который нужно выполнить, это усреднить результат, поскольку это не делается алгоритмом (вероятно, это будет один из способов сделать его быстрее )

Под средним я подразумеваю, что нужно разделить результат на количество точек.

Я хотел бы закрыть этот пост, показав окончательный результат и немного того, что было необходимо для этого. После этого я, как обычно, выложу ссылку на код.

The final animation

В этом примере я также использовал задачу коммивояжера (TSP), чтобы объединить все буквы в одну непрерывную форму.

Короткая заметка:
Если вас интересует проблема коммивояжера или оптимизация в целом я настоятельно рекомендую проверить мои ранние посты:

• Solving the TSP using mixed integer programming
• Doing the same thing but with neighborhoods
• Building a constraint solver (in several steps)

Ладно вернемся к посту:

Нам нужно использовать задачу коммивояжера, чтобы избежать больших скачков от одной буквы к другой, поскольку рисунок непрерывен. Потенциально можно избавиться от всех скачков вместе, имея версию ряда Фурье с разрывом.

В коде, который я прикрепил ниже, я также использую одну переменную для количества точек и кругов, что не самая лучшая идея, и я изменил ее для рисования треугольника, но в этом случае я хотел ограничить время работы эвристического решателя TSP.

Выводы

Прежде чем я закончу этот пост рекламой для себя, позвольте мне немного подытожить то, что мы узнали за последние 20 минут.

1. Основные изменения в Javis v0.3
2. Рисование простых фигур и применение к ним действий
3. Что можно сделать с рядами Фурье в 1D и 2D
4. Формула и объяснение ряда Фурье
5. Интуиция, как получить необходимые константы

Ссылки

• Видео-объяснение от 3blue1brown But what is a Fourier series?
• Видео-кодинг от The Coding Train Coding Challenge: #125 Fourier Series
• Пост на betterexplained.com An Interactive Guide To The Fourier Transform
• Рисовалка надписи на гитхабе

Далее идут ссылочки на реквизиты и персональные благодарности за патронаж, так что не скупимся и поддерживаем товарища Оле Крюгера ведущего интересный блог.