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

Худший код тот, что написан до меня

Причина, по которой [разработчики] считают старый код бардаком, заключается в кардинальном, фундаментальном законе программирования: читать код труднее, чем писать его.

Джоэл Спольски Грабли, на которые не стоит наступать

Вот почему вы их боитесь

Сравнение зелёных и коричневых языков

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

Java, C, C++, C#, Python, PHP, JavaScript, Swift, Perl, Ruby, Assembly, R, Objective-C, SQL

Зелёный язык: язык, который вы с большей вероятностью будете использовать в новом проекте.

Go, Rust, TypeScript, Kotlin, Julia, Dart, Scala и Haskell

Страшные коричневые языки

Любимые зелёные языки

Ещё один недостаток человеческого характера заключается в том, что все хотят строить и никто не хочет заниматься обслуживанием.

Курт Воннегут

Устранение предвзятости

Цикл хайпа языков программирования

Так как Python один из лучших языков программирования, мы в EDISON часто используем его в сложных интересных проектах.


Мы разработали приложения и сайты Московского ювелирного завода.

Полное тестирование новой версии сайта было осуществлено на Python и Django.

Дзен Python VS жадности Julia

Мы ненасытны: мы хотим большего.

Нам нужен язык с открытым исходным кодом, со свободной лицензией. Мы хотим скорость Си с динамикой Ruby. Мы хотим, чтобы язык был гомоиконическим, с настоящими макросами, такой как Lisp, но с очевидными, знакомыми математическими сущностями, такой как Matlab. Мы хотим настолько же удобное для общего программирования, как Python, такое же простое для статистики, как R, такое же естественное для обработки строк, как Perl, такое же мощное для линейной алгебры, как Matlab, и способное объединять все эти возможности под одной оболочкой. Нечто простое в освоении, но при этом радующее самых серьёзных хакеров. Мы хотим, чтобы язык был интерактивным, и мы хотим, чтобы он был компилируемым.

За что разработчики любят Julia

Многосторонность

Скорость

Сообщество

Преобразование кода

Библиотеки

Динамическая и статическая типизация

Данные: инвестируйте в небольшие вещи, имеющие большой потенциал

Подытожим: знание Julia может стать конкурентным преимуществом

Признавая свою непомерную алчность, мы всё также жаждем заполучить всё. Около двух с половиной лет назад мы решили создать язык нашей жадности. Он не завершён, но пришло время для версии 1.0 созданный нами язык называется Julia. Он уже удовлетворяет 90% наших неблагодарных требований, и теперь нужны неблагодарные требования остальных, чтобы развиваться дальше. Итак, если вы также непомерно жадный и излишне требовательный программист, мы хотим, чтобы вы попробовали Julia.

Фармакокинетическая модель это часто упрощенное математическое описание того, как изменяется концентрация(количество) исследуемого вещества во времени в биологической системе. Упрощенное потому, что организм сложный механизм в котором абсорбция, распределение, метаболизм и выведение веществ (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

Сейчас мне хочется поделиться своими выводами сделанными после нескольких бесед, в которых я участвовал на JuliaCon 2020.

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

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

В течение многих лет я бы с этим согласился, но после JuliaCon 2020, я думаю, мы можем с уверенностью заявить, что

Джулия готова идти в производство!

Позвольте мне теперь привести список ключевых (на мой взгляд) презентаций, сделанных во время JuliaCon 2020, которые заставили меня сделать этот вывод.

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

Создание микросервисов и приложений

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

Еще будет интересно:

• Доклад по характеристикам приложений, где Кристоффер Карлссон покажет как создать исполняемые файлы, которые могут быть запущены на машинах, на которых не установлена Julia.
• Презентация Julia for scripting, в ходе которой [Фредрик Экре]() обсуждает лучшие практики использования Julia в контекстах, где вам нужно много раз выполнять короткие фрагменты кода.
• Доклад Genie.jl, в котором Адриан Сальчану показывает зрелый, стабильный, производительный и многофункциональный фреймворк веб-разработки на Julia.

Управление зависимостями

Пара докладов Pkg.update() и What's new in Pkg показывает, что в настоящее время Julia имеет лучшие в своем классе функциональные возможности для управления зависимостями корпоративного уровня для ваших проектов. Список предоставляемых функций настолько велик, что здесь трудно перечислить их все.

Позвольте мне только упомянуть один конкретный инструмент в этой экосистеме BinaryBuilder.jl, который позволит взять программное обеспечение, написанное на компилируемых языках, таких как C, C++, Fortran, Go или Rust, и построить предварительно скомпилированные артефакты, которые могут быть легко использованы из пакетов Julia (что означает, что не нужна никакая компиляция на стороне клиента, когда вы устанавливаете пакеты, имеющие такие зависимости).

Интеграция с внешними библиотеками

Естественной темой, связанной с управлением зависимостями, является интеграция Julia с внешними инструментами. Эта область функциональности действительно зрелая. Вот список докладов, которые охватывают эту тему:

• Сделайте ваш Julia-код более быстрым и совместимым с не-Julia-кодом
• Обертывание библиотеки C++ с помощью CxxWrap.jl
• Julia and C++: технический обзор CxxWrap.jl
• Представляем макрос @ccall
• Интеграция Julia в R с JuliaConnectoR
• Интеграция Julia и Javascript с помощью расширения Node.js

Здесь стоит добавить, что Джулия уже много лет отлично интегрируется с Python, см. JuliaPy.

Хорошим всеобъемляющим примером выполнения некоторой реальной работы в Julia, требующей интеграции, является создание многоканальной беспроводной настройки динамиков, где показано, как легко сшивать всякие штуки вместе (и, в частности, с помощью ZMQ.jl, Opus.jl, PortAudio.jl, и DSP.jl).

Еще один интересный доклад, демонстрирующий возможности интеграции, это JSServe: Websites & Dashboards в Julia, который показывает высокопроизводительный фреймворк для легкого объединения интерактивных графиков, markdown, виджетов и простого HTML/Javascript в Jupyter / Atom / Nextjournal и на веб-сайтах.

Инструментарий разработчика

В паре замечательных докладов Juno 1.0 и Using VS Code рассказывается, что текущая поддержка IDE для Джулии в VS Code первоклассна. У вас есть все для полного счастья: анализ кода (статический и динамический), отладчик, рабочие области, интеграция с ноутбуками Jupyter и удаленные возможности.

Управление рабочими процессами в машинном обучении

Я не хочу охватывать множество различных алгоритмов ML, доступных в Julia изначально, так как их просто слишком много (и если чего-то не хватает, вы можете легко интегрировать его см. раздел возможности интеграции выше).

Однако в дополнение к конкретным моделям вам нужны фреймворки, позволяющие управлять рабочими процессами ML. В этой области есть две интересные презентации: MLJ: A machine learning toolbox for Julia, и AutoMLPipeline: A ToolBox for Building ML Pipelines. По моему опыту, такие инструменты имеют решающее значение, когда вы хотите перейти с вашими ML-моделями из песочницы дата-саянтиста в реальное производственное использование.

Выводы

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

• модели корпоративного уровня для создания микросервисов, приложений и
  надежные инструменты управления зависимостями,
• очень гибкие и мощные возможности для интеграции Julia с существующими кодовыми базами, которые не были написаны на Julia,
• отличный инструмент разработчика в VSCode,
• зрелые пакеты, которые помогут вам создать производственный код для развертывания ML-решений,

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

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

Готова Ли Джулия для прода? Вопросы и ответы с Богумилом Каминским

Статья профессора Каминского вызвала немного ажиотажа на Hacker News. Некоторые комментаторы высказывали сомнения в том, что Джулия может считаться готовой к производству из-за, в частности, документации, пакетов, инструментов и поддержки.

InfoQ решило поговорить с профессором Камински, чтобы лучше понять его позицию.

InfoQ: Не могли бы вы описать ваш быкграунд и вашу связь с Julia?

Bogumi Kamiski: Я профессор Варшавской школы экономики SGH, Польша, работаю в области операционных исследований и имитационного моделирования.
Я нахожусь в топ-5% участников языка Julia по количеству коммитов, внес значительный вклад в экосистему данных Julia и, в частности, один из основных мейнтейнеров DataFrames.jl.
Я занимаю второе место по ответу на тег [julia] в StackOverflow. Прежде чем перейти на полный рабочий день в академию, я более 10 лет руководил командой из нескольких сотен разработчиков и аналитиков, развертывающих проекты BI/DWH/data science для крупнейших польских корпораций и учреждений.

InfoQ: Основная аргументация в вашей статье, по-видимому, заключается в том, что экосистема Julia теперь достигла уровня зрелости, который делает ее готовой к производству. Не могли бы вы подробнее остановиться на этом моменте? Что препятствовало внедрению Джулии в производство и каковы наиболее значительные достижения, которые устранили эти блокирующие моменты?

Kamiski: Здесь очень важно определить, что я понимаю под готовностью Джулии к проду.

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

В прошлом это было серьезной проблемой. Основные пакеты, как правило, часто меняли свой API; учебные пособия, созданные год назад, теперь не будут работать без получения обновлений и т. д. Это было естественное состояние, когда развивались язык и экосистема. Теперь я вижу, что эти вещи значительно изменились, особенно для ядра языка Julia, но аналогичные вещи происходят и в экосистеме пакетов.

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

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

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

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

Почему это так важно? Ну, вы можете легко "отправить" проект Julia и ожидать, что любой человек в любой среде должен быть в состоянии относительно легко запустить его. Конечно, будут угловатые моменты; но мой опыт ежедневного использования Linux и Windows 10 заключается в том, что все должно просто работать на обеих платформах.

Если вы делаете проект в Julia, вы можете ожидать, что либо есть пакет, который делает то, что вы хотите, либо вы можете легко использовать код, написанный на C или Python, и заставить его работать. В своем посте я хотел подчеркнуть, что мы находимся именно в этом состоянии. В качестве примера, основанного на чем-то, над чем я работал на этой неделе, у Джулии есть отличный пакет LightGraphs.jl для работы с графиками, но мои сотрудники используют Python и предпочитают использовать igraph. Вот пример кода с использованием графика (взято из учебника для графика в Python):

import igraph as igg = ig.Graph()g.add_vertices(3)g.add_edges([(0,1), (1,2)])g.add_edges([(2, 0)])g.add_vertices(3)g.add_edges([(2, 3), (3, 4), (4, 5), (5, 3)])g.pagerank()

Теперь вы спросите, как бы выглядел эквивалент в Julia. А вот так:

using PyCallig = pyimport("igraph")g = ig.Graph()g.add_vertices(3)g.add_edges([(0,1), (1,2)])g.add_edges([(2, 0)])g.add_vertices(3)g.add_edges([(2, 3), (3, 4), (4, 5), (5, 3)])g.pagerank()

Как видите, он идентичен. Конечно, не во всех случаях это так просто, так как, например, словари в Julia и Python имеют разные синтаксисы, но это общее правило о том, как все работает. У вас даже есть вкладка-завершение и прямой доступ к строкам документов.

Что это означает на практике? Если вы делаете проект, вы не застреваете с мыслью: "могу ли я использовать Джулию, так как через три месяца, возможно, мне понадобится что-то в проекте, чего еще нет в Джулии"? Но скорее вы знаете: "я могу относительно безопасно использовать Julia, так как в настоящее время многие пакеты общего назначения уже существуют, и даже если чего-то не хватает, я могу просто использовать это с другого языка, и это будет относительно безболезненно, независимо от того, является ли это C/Python/R/...".

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

Теперь позвольте мне уточнить, что я не вижу в определении термина "готовый к производству". Я не рассматриваю Джулию как язык, который лучше всего подходит для любого вида проекта. Каждый язык программирования имеет свою нишу, и ниша Julia это высокопроизводительные вычисления / наука о данных (или как вы это там называете). Если вам нужны компактные бинарные файлы наверняка Julia не является лучшим вариантом. Если вы хотите разрабатывать приложения для Android тут тоже посоветую смотреть в другом направлении.

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

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

Kamiski: Процитирую то, что я написал в начале своего поста:

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

Это основа моего поста, то есть заниматься наукой о данных, но не на "вашем заднем дворе" или "академической исследовательской лаборатории", а скорее в "деловой среде". Как я уже объяснял выше, ни я, ни, вероятно, кто-либо из разработчиков, участвующих в Julia, не утверждают, что это "универсам" для любого вида девелоперского проекта.

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

Теперь, что касается таких вещей, как документация, пакеты, инструменты и поддержка конечно, это должно и может быть улучшено. И я согласен, что более зрелые экосистемы, такие как R/Python/Java, в среднем имеют здесь лучший охват. Например, как соразработчик DataFrames.jl, я могу сказать вам, что большинство последних PR связаны с документацией. Но я бы не стал недооценивать сообщество Джулии здесь. В общем, если у вас есть какой-либо вопрос и вы разместите его на SO, Julia Discourse или Julia Slack, вы можете получить ответ обычно в течение нескольких минут, самое большее часов. Вот типичная история такого рода: люди обычно очень отзывчивы, и ошибки исправляются довольно быстро.

Если бы меня попросили назвать главный сомнительный момент касательно Джулии, то это было бы наличие достаточного количества людей, квалифицированных в этом языке в общем сообществе разработчиков. Я могу понять владельцев продуктов / менеджеров проектов и их чувство риска того, что они не смогут найти достаточно людей для работы над своими проектами, как только они возьмут на себя смелость начать работать на Julia. Однако здесь я убежден, что ситуация быстро улучшается. В прошедшем JuliaCon 2020 приняли участие более 20 000 участников. Кроме того, есть много ресурсов, уже доступных бесплатно в интернете.

InfoQ: Помимо вопроса о том, готова ли Джулия к проду или нет, или для каких областей она предназначена, каковы, на ваш взгляд, основные сильные стороны языка? Рассматриваете ли вы его как замену Python, R или любому другому языку, по крайней мере в области научных вычислений и науки о данных?

Kamiski: Я думаю, что здесь снова лучше всего процитировать последний опрос разработчиков Julia, слайды с 8 по 11. Я бы сосредоточился на трех главных вещах из слайда 8:

• Скорость здесь ситуация относительно проста. Возьмите любой зрелый пакет, такой как TensorFlow или PyTorch, который требует производительности; они написаны в основном на C++. А Python это всего лишь тонкая оболочка вокруг ядра C++. Теперь возьмем Flux.jl или Knet.jl. Они по существу реализованы в чистом виде. Итак, суть такова: если вам нужно, чтобы ваш код работал быстро, в то же время используя преимущества языка высокого уровня, то Julia это естественный выбор. Кроме того, как объяснялось выше, если есть внешняя библиотека, которая очень быстра, и вы хотите использовать ее, обычно это относительно легко сделать.

  
  

• Простота использования есть множество проявлений этого аспекта, но мой опыт показывает, что когда кто-то получает представление о принципах Джулии, он очень хорошо сравнивается с точки зрения синтаксиса и дизайна, например, с R/Python при выполнении вычислений. Язык был разработан вокруг поддержки этих видов задач в лучшем виде. И это не только синтаксис это также выбор того, когда вещи будут происходить явно, а когда нет (например, с трансляцией). Исходя из моего опыта, это делает код Джулии простым в обслуживании, а хорошо написанный код в значительной степени самодокументируется.

  
  

• Код является открытым и может быть изменен этот аспект имеет два измерения. Прежде всего, большинство пакетов Julia лицензированы MIT, что часто очень приветствуется в корпоративных средах. Во-вторых поскольку большинство пакетов написано на Julia, то если вам не нравится, как что-то работает вы просто модифицируете это сами (что гораздо проще, чем в R/Python, где, скорее всего, то, что вам нужно изменить, написано, например, на C, C++, Fortran).

  
  

Люди часто спрашивают меня, вижу ли я Джулию в качестве замены R / Python. И я думаю об этом так::

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

  
  

• Если у вас есть новая задача, не требующая критической производительности просто используйте язык, который вы знаете лучше всего (и если для вас часто актуален пункт 1 как для меня вы можете смело выбрать для этого Julia).

  
  

• Если у вас много устаревшего кода R/Python и вы им довольны просто придерживайтесь его и помните, что если у вас есть некоторые критически важные для производительности части, они могут быть относительно легко переписаны на Julia и интегрированы обратно с вашей исходной кодовой базой (я сделал много таких проектов).

  
  

InfoQ: Как вы видите эволюцию Джулии?

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

Если бы я назвал здесь некоторые основные изменения, то это были бы:

• Улучшение поддержки многопоточности. Я думаю, что это действительно актуально для основных случаев использования Julia. Так-то поддержка есть, но все же здесь многое можно улучшить. Точно так же следует ожидать улучшения в обработке GPU/TPU.

  
  

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

  
  

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

  
  

• Рост сообщества я считаю, что число людей, заинтересованных в Джулии, значительно увеличивается, о чем свидетельствует, например, количество участников JuliaCon2020. Это означает, что: а) в будущем будет легче нанимать качественных разработчиков Julia, если они понадобятся, и б) это создаст положительную обратную связь для качества экосистемы Julia, поскольку все больше людей сообщают о проблемах и участвуют в них. Опять же, как сопровождающий DataFrames.jl я наблюдаю этот сдвиг: люди, которые никогда не были в "ядре" разработки пакета, открывают вопросы / делают PR и обсуждают функциональные возможности в социальных сетях.

  
  

Если вы заинтересовались языком Julia, все доклады с JuliaCon 2020 доступны на YouTube.

Продолжаем рубрику тем для первого свидания. На сегодняшней повестке дня упрощение схем для квантовых программ методами ZX-исчисления.

Далее от имени автора

Прежде всего, я хотел бы поблагодарить моих наставников Роджера Ло и Цзиньго Лю за то, что они курировали меня во время Google Summer of Code 2020. В этом проекте GSoC я разрабатываю новый пакет Julia ZXCalculus.jl, который реализует ZX-исчисление, и интегрирую его в качестве инструмента упрощения схем и квантовых программ в YaoLang.jl DSL (предметно-ориентированный язык) для Yao.jl. Yao.jl современный квантовый симулятор со множеством продвинутых функций, такими как автоматическое дифференцирование и символьные вычисления. Как пользователь "Yao.jl", я рад возможности внести свой вклад в его развитие. В этом блоге я подытожу результаты проделанной работы.

Квантовые схемы и ZX-исчисление

Квантовое программирование использует принципы квантовой механики для вычислений. Для описания квантовых алгоритмов мы часто используем модель квантовой схемы, которая является аналогом модели классической логической схемы. Квантовые схемы состоят из набора основных квантовых элементов, таких как вентили Паули, Адамара, T, CNOT.

В общем, будет много эквивалентных квантовых схем, представляющих одну и ту же операцию. С точки зрения квантового оборудования, нужно найти схему, которая минимизирует использование аппаратных ресурсов. К сожалению, поиск наиболее желаемой схемы является очень сложной задачей в плане вычислительной сложности, так как проверка эквивалентности между двумя квантовыми схемами является coQMA-трудной (квантовый аналог coNP) [1], [2]. Несмотря на это, у нас все еще есть некоторые эвристические эффективные алгоритмы для упрощения квантовых схем. В этом блоге я представлю мощный инструмент для этой проблемы ZX-исчисление.

Квантовая схема это граф, представляющий тензорные и матричные произведения. Например

Кроме того, квантовые схемы можно рассматривать как особый тип тензорной сети. В ZX-исчислении мы вводим два вида тензоров: Z-пауков и X-пауков. Тензорные сети, состоящие из этих тензоров, называются ZX-диаграммами. Поскольку все основные квантовые элементы могут быть представлены ZX-диаграммами, мы можем преобразовать все квантовые схемы в ZX-диаграммы.

Более того, по определению Z-пауков и X-пауков, существуют основные эквивалентные правила, с помощью которых можно упростить ZX-диаграммы.

После упрощения ZX-диаграмм мы должны преобразовать их обратно в схемы. В работе [3] разработаны алгоритмы для приведенной выше схемы.

Во время GSoC 2020 я собираюсь добиться реализации на Джулии ZX-исчисления. В следующих разделах я кратко объясню причины и методы этого проекта.

Цель GSoC проекта

Основная цель этого проекта заключается в реализации ZX-исчисления и связанного с ним алгоритма на языке Julia. И мы выпустим пакет Julia ZXCalculus.jl. Существует Python-реализация ZX-исчисления, PyZX. Как и PyZX, ZXCalculus.jl должна предоставлять API для перевода на ZX-диаграммы, извлечения схем из ZX-диаграмм, упрощения и визуализации. Большинство функций PyZX будут реализованы в ZXCalculus.jl. Кроме того, мы интегрируем ZXCalculus.jl в квантовый компилятор YaoLang.jl и реализуем механизм упрощения схем на уровне компиляции.

В YaoLang.jl можно определить гибридные квантовые программы с классической информацией (например, классические параметры, классические потоки управления). Эти квантовые программы будут преобразованы в промежуточное представление (Yao IR), которое содержит квантовые схемы наряду с формой SSA (Статическая форма единого назначения) классической программы. Эти схемы будут упрощены с помощью ZXCalculus.jl. Затем YaoIR будет скомпилирован в бэкенд-инструкции, такие как QASM или инструкции симулятора Yao.

Для достижения наилучшей производительности при компиляции и обеспечения возможности дальнейшей настройки процедуры упрощения необходима чистая реализация ZX-исчисления на языке Julia.

Методы

Метод извлечения схем является каноническим методом алгоритмов упрощения схем на основе ZX-исчисления. Я тут буду в основном обсуждать его реализую.

Структуры данных

Для реализации этих алгоритмов нам необходимо реализовать структуры данных для хранения ZX-диаграмм и правила изменения ZX-диаграмм. ZX-диаграмму можно рассматривать как мультиграф с дополнительной информацией о его вершинах. Каждая вершина может быть представлена чем-то из Z-пауков, X-пауков или H-коробок. В исходных ZX-диаграммах допускаются открытые ребра. Для простоты мы добавили 2 специальных вида вершин, входную границу и выходную границу для записи этих ребер. Мы используем представление списка смежности для хранения мультиграфов. Список смежности хранится в виде "Dict", ключами которого являются идентификаторы вершин.

struct Multigraph{T<:Integer} <: AbstractMultigraph{T}    adjlist::Dict{T, Vector{T}}end

Поскольку нам нужно переписать мультиграфы с помощью правил ZX-исчисления, будет удобнее фиксировать идентификаторы вершин. Поэтому мы используем Dict вместо Array. Еще бывают фазы для Z- и X-пауков. Нам нужен еще один словарь для хранения этих фаз. Структура данных ZX-диаграммы аналогична следующей:

struct ZXDiagram{T<:Integer, P} <: AbstractZXDiagram{T, P}    mg::Multigraph{T}    st::Dict{T, SpiderType.SType}    ps::Dict{T, P}    ...end

В статье предложены графоподобные ZX-диаграммы. Они включают в себя только Z-пауков и имеют разные типы граней: адамаровские и неадамаровские. Мы также можем использовать мультиграфы для представления графоподобных ZX-диаграмм и различные кратности ребер для представления различных типов ребер. Мы определяем ZXGraph для графоподобных ZX-диаграмм так:

struct ZXGraph{T<:Integer, P} <: AbstractZXDiagram{T, P}    mg::Multigraph{T}    st::Dict{T, SpiderType.SType}    ps::Dict{T, P}    ...end

С этого момента мы создали структуры данных, которые нам нужны.

Правила преобразований

Чтобы упростить ZX-диаграммы, нам нужно реализовать правила перезаписи в ZX-исчислении. Существует два шага: сопоставление правил и переписывание ZX-диаграмм с соответствующими правилами. Мы использовали для этого благословенную джулианскую множественную диспетчеризацию. API-интерфейсы выглядят как-то так:

match(r::AbstractRule, zxd::AbstractZXDiagram)rewrite!(r::AbstractRule, zxd::AbstractZXDiagram{T, P}, matches::Vector{Match{T}}) where {T, P}

Здесь match вернет все совпадающие вершины, которые будут сохранены в структуре

struct Match{T<:Integer}    vertices::Vector{T}end

А rewrite! попытается использовать правило, чтобы переписать ZX-диаграмму на всех совпадающих вершинах. Так как некоторые совпадающие вершины могут стать неудовлетворяющими соответствующему правилу после переписывания некоторых вершин, нам нужно проверить, доступно ли еще правило для совпадающих вершин.

check_rule(r::AbstractRule, zxd::AbstractZXDiagram{T, P}, vs::Vector{T}) where {T, P}

Для упрощения ZX-диаграммы с помощью правила мы определили simplify! так:

function simplify!(r::AbstractRule, zxd::ZXDiagram)    matches = match(r, zxd)    while length(matches) > 0        rewrite!(r, zxd, matches)        matches = match(r, zxd)    end    return zxdend

На этапе упрощения [3], мы должны конвертировать ZX-диаграмму в графоподобный вариант с правилами i1, i2, h, и f. После, упрощаем графоподобную ZX-диаграмму локальным дополнительным правилом и правилом замещения. Мы можем выполнить эти шаги с помощью вышеуказанных функций. Упрощенные графоподобные ZX-диаграммы будут просто маленькими скелетами по сравнению с первоначально большой схемой. Единственное, что остается это извлечение схем из ZX-диаграмм.

Извлечение схем

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

Существует особенность для любой ZX-диаграммы квантовой схемы. Мы можем определить на ней g-поток. С помощью g-потока мы можем узнать порядок всех пауков. Правила, которые мы использовали для упрощения ZX-диаграммы, сохранят это хорошее свойство. Это дает нам возможность извлекать схемы из упрощенных ZX-диаграмм.

Кроме того, любая графоподобная ZX-диаграмма локально представима по CNOT + H. Мы можем извлечь эти схемы CNOT с гауссовским устранением над $F_2$. Вместе с остальными H-, CZ-гейтами и Z-вращениями мы получаем схему, эквивалентную исходной.

Demo

Демонстрационная схема взята из приложения к статье [3]. Все коды можно найти в разделе examples\ex1.jl у ZXCalculus.jl.

# this is the original circuit.zxd = generate_example()ZXplot(zxd)

# convert a ZX-diagram to a graph-like ZXdiagramzxg = ZXGraph(zxd)ZXplot(zxg, linetype = "curve")

# simplify with local complementary rule simplify!(Rule{:lc}(), zxg)ZXplot(zxg, linetype = "curve")

# simplify with pivoting rulesimplify!(Rule{:p1}(), zxg)ZXplot(zxg, linetype = "curve")

# removing all Paulis adjancent to boundariesreplace!(Rule{:pab}(), zxg)ZXplot(zxg, linetype = "curve")

# extract circuitcir = circuit_extraction(zxg)ZXplot(cir)

Наконец, мы получили упрощенную схему.

• [1]: "NON-IDENTITY-CHECK" IS QMA-COMPLETE
• [2]: Exact Non-identity check is NQP-complete
• [3]: Graph-theoretic Simplification of Quantum Circuits with the ZX-calculus
• [4]: Reducing T-count with the ZX-calculus

Метапрограммирование, квантовая компиляция и оптимизация схем при квантовой компиляции

В прошлом месяце я в основном работал над интеграцией с YaoLang.jl для ZXCalculus.jl. Это был мой первый опыт работы с метапрограммированием Джулии на практике. Я хочу поблагодарить своего наставника Роджера Ло за то, что он научил меня основным понятиям и полезным методам метапрограммирования.

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

Как Julia работает

Все коды Джулии по сути строки (Strings), которые хранятся на диске. Когда мы запускаем коды Julia, Джулия сначала разбирает этот код на AST (абстрактные синтаксические деревья). AST будет храниться в виде выражений в структуре данных Expr в Julia. На этом уровне синтаксического анализа мы называем эти выражения IR (промежуточные представления) поверхностного уровня.

julia> s = "1 + 1 * 2""1 + 1 * 2"julia> ex = Meta.parse(s):(1 + 1 * 2)julia> dump(ex)Expr  head: Symbol call  args: Array{Any}((3,))    1: Symbol +    2: Int64 1    3: Expr      head: Symbol call      args: Array{Any}((3,))        1: Symbol *        2: Int64 1        3: Int64 2

В этом примере ex это AST. Его head это :call, что означает, что это вызов функции. Он вызывает функцию + с аргументами 1 и еще одним Expr.

Затем опускаемся на уровень. На этом уровне макросы будут расширены, а "синтаксический сахар" Джулии будет преобразован в вызовы функций. Например, a[i] будет заменено на getindex(a, i). После серии преобразований, IR поверхностного уровня будет преобразован в SSA (Статическая форма единого назначения). IR также называется "пониженным" IR. В SSA IR каждая переменная может быть назначена только один раз. В Julia мы можем использовать макрос @code_lowed, чтобы увидеть SSA IR объекта Expr:

julia> @code_lowered 2 + 3CodeInfo(1  %1 = Base.add_int(x, y)      return %1)

Джулия сделает вывод типа на SSA IR и оптимизирует его. А затем преобразует его в LLVM коды. Мы можем использовать макросы @code_typed и @code_llvm, чтобы увидеть эти IRs.

julia> @code_typed 2 + 3CodeInfo(1  %1 = Base.add_int(x, y)::Int64      return %1) => Int64julia> @code_llvm 2 + 3;  @ int.jl:53 within '+'; Function Attrs: uwtabledefine i64 @"julia_+_15307"(i64, i64) #0 {top:  %2 = add i64 %1, %0  ret i64 %2}

Наконец, LLVM преобразует эти коды в собственные машинные коды.

julia> @code_native 2 + 3        .text;  @ int.jl:53 within '+'        pushq   %rbp        movq    %rsp, %rbp        leaq    (%rcx,%rdx), %rax        popq    %rbp        retq        nopw    (%rax,%rax); 

Вот пикча из JuliaCon 2018, которая демонстрирует, как работает компилятор Julia.

Как работает YaoLang.jl

Цель YaoLang.jl построить удобный квантовый компилятор для гибридных квантово-классических программ в Julia. То есть, используя только несколько макросов и добавляя их к собственным функциям Julia, можно определить квантовые программы. В YaoLang.jl функция, украшенная макросом device, будет рассматриваться как функция с квантовыми операциями. В этих функциях макросы @ctrl, @measure, @expect и "синтаксический сахар" locs => gate доступны для определения квантовых операций. Например, представим схему для квантового преобразования Фурье n кубитов:

@device function qft(n::Int)    1 => H    for k in 2:n        @ctrl k 1 => shift(2 / 2^k)    end    if n > 1        2:n => qft(n - 1)    endend

Подобно процедурам компиляции Julia, макрос device будет анализировать функцию в IR поверхностного уровня в YaoLang.jl. Затем все макросы и синтаксический сахар для квантовых операторов будут заменены вызовами функций. Эти вызовы функций будут помечены меткой :quantum. Теперь IR поверхностного уровня будет преобразован в пониженный SSA IR. В YaoLang.jl SSA IR будет храниться в структуре данных YaoIR.

Остальные части это оптимизация YaoIR и преобразование из него в коды аппаратного уровня. ZXCalculus.jl предназначен для оптимизации квантовых схем и должен быть интегрирован на уровне оптимизации.

Интеграция ZXCalculus.jl

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

1. преобразовать его в ZX-диаграмму
2. Упростить ZX-диаграмму
3. Перевести упрощенную ZX-диаграмму обратно в YaoIR

Второй шаг уже реализован в ZXCalculus.jl. Нам нужно только реализовать преобразование между ZXDiagram и YaoIR.

YaoIR в ZXDiagram

Поскольку каждый YaoIR является SSA, мы можем просмотреть все утверждения, чтобы получить информацию о вентилях и их местоположениях. Для этого окинем взором самую обширную локацию по числу кубитов. Чтобы построить соответствующую ZX-диаграмму, построим пустую схему и последовательно протолкнем в нее гейты при обходе YaoIR. А код такой:

function ZXDiagram(ir::YaoIR)    if ir.pure_quantum        n = count_nqubits(ir)        circ = ZXDiagram(n)        stmts = ir.body.blocks[].stmts        for stmt in stmts            # push gates        end        return circ    endend

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

ZXDiagram в YaoIR

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

Их можно выудить из информации о расположении ZXDiagram. Так мы узнаем, как пауки сортируются от входа к выходу и высмотрим расположение кубитов для каждого паука. Если паук имеет степень 2, он представляет собой однокубитовый гейт. В противном случае он представляет собой вентиль с несколькими кубитами. Пройдя всех пауков от входа к выходу, мы можем получить последовательность квантовых вентилей. И тогда мы сможем построить новый YaoIR с этой последовательностью.

Параметры оптимизации

В "ZXCalculus.jl" существует несколько методов упрощения схем. И мы предлагаем реализовать другие методы упрощения, которые не основаны на ZX-исчислении. Необходимо предоставить пользователю возможность выбирать, какие методы оптимизации будут применяться.

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

@device optimizer = [opt...] function my_circuit(args...)    ...end

Так optimizer вполне представим как подмножество [:zx_clifford, zx_teleport]. И мы добавим больше методов в будущем.

Примеры

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

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

using YaoLang@device function test_cir()    5 => H    5 => shift(0.0)    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    5 => shift($(1/4*))    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    4 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 1 4 => X    4 => shift($(7/4*))    1 => shift($(1/4*))    @ctrl 1 4 => X    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 3 5 => X    5 => shift($(1/4*))    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 3 5 => X    4 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 3 4 => X    4 => shift($(7/4*))    5 => H    3 => shift($(1/4*))    @ctrl 3 4 => X    5 => shift(0.0)    @ctrl 4 5 => X    5 => H    5 => shift(0.0)    @ctrl 3 5 => X    5 => shift($(7/4*))    @ctrl 2 5 => X    5 => shift($(1/4*))    @ctrl 3 5 => X    5 => shift($(7/4*))    @ctrl 2 5 => X    3 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 2 3 => X    3 => shift($(7/4*))    5 => H    2 => shift($(1/4*))    @ctrl 2 3 => X    5 => shift(0.0)    @ctrl 3 5 => X    5 => H    5 => shift(0.0)    @ctrl 2 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    5 => shift($(1/4*))    @ctrl 2 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    2 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 1 2 => X    2 => shift($(7/4*))    5 => H    1 => shift($(1/4*))    @ctrl 1 2 => X    5 => shift(0.0)    @ctrl 2 5 => X    @ctrl 1 5 => Xendcir = test_cir()@device optimizer = [:zx_teleport] function teleport_cir()    5 => H    5 => shift(0.0)    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    5 => shift($(1/4*))    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    4 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 1 4 => X    4 => shift($(7/4*))    1 => shift($(1/4*))    @ctrl 1 4 => X    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 3 5 => X    5 => shift($(1/4*))    @ctrl 4 5 => X    5 => shift($(7/4*))    @ctrl 3 5 => X    4 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 3 4 => X    4 => shift($(7/4*))    5 => H    3 => shift($(1/4*))    @ctrl 3 4 => X    5 => shift(0.0)    @ctrl 4 5 => X    5 => H    5 => shift(0.0)    @ctrl 3 5 => X    5 => shift($(7/4*))    @ctrl 2 5 => X    5 => shift($(1/4*))    @ctrl 3 5 => X    5 => shift($(7/4*))    @ctrl 2 5 => X    3 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 2 3 => X    3 => shift($(7/4*))    5 => H    2 => shift($(1/4*))    @ctrl 2 3 => X    5 => shift(0.0)    @ctrl 3 5 => X    5 => H    5 => shift(0.0)    @ctrl 2 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    5 => shift($(1/4*))    @ctrl 2 5 => X    5 => shift($(7/4*))    @ctrl 1 5 => X    2 => shift($(1/4*))    5 => shift($(1/4*))    @ctrl 1 2 => X    2 => shift($(7/4*))    5 => H    1 => shift($(1/4*))    @ctrl 1 2 => X    5 => shift(0.0)    @ctrl 2 5 => X    @ctrl 1 5 => Xendtp_cir = teleport_cir()

С помощью пакета YaoArrayRegister.jl, мы можем вычислить матрицу для каждой схемы.

using YaoArrayRegistermat = zeros(ComplexF64, 32, 32)for i = 1:32    st = zeros(ComplexF64, 32)    st[i] = 1    r0 = ArrayReg(st)    r0 |> cir    mat[:,i] = r0.stateendtp_mat = zeros(ComplexF64, 32, 32)for i = 1:32    st = zeros(ComplexF64, 32)    st[i] = 1    r1 = ArrayReg(st)    r1 |> tp_cir    tp_mat[:,i] = r1.stateend

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

sum(abs.(mat - tp_mat) .> 1e-14) == 0

Итого

Во время второго этапа кодирования я реализовал преобразование между "ZXDiagram" и "YaoIR", что обеспечивает интеграцию ZXCalculus".jl с YaoLang.jl. Кроме того, документация теперь доступна здесь. И во время испытания ZXCalculus".jl с помощью YaoLang.jl и YaoArrayRegister.jl я нашел несколько ошибок в реализации извлечения контуров и фазовой телепортации. Эти ошибки уже исправлены.

На следующем этапе я буду работать над компиляцией кодов OpenQASM в коды YaoIR. Это позволит нам читать схемы из кода OpenQASM. А также протестирую результативность ZXCalculus.jl на некоторых эталонных схемах и сравню его с PyZX.

Упрощение квантовой схемы

Предположим, что у нас есть квантовая схема с 24 гейтами, как описано выше. Мы можем определить эту схему с помощью YaoLang.jl используя макрос @device. YaoLang.jl это компилятор для гибридных квантово-классических программ, которые очень практичны в нынешнюю эпоху NISQ (noisy intermediate-scale quantum). Кроме того, YaoLang.jl интегрирован с ZXCalculus.jl.

julia> using YaoLang;julia> @device optimizer = [:zx_teleport] function demo_circ_simp()           1 => shift($(7/4))           1 => H           1 => Rx($(/4))           4 => H           @ctrl 1 4 => Z           @ctrl 4 1 => X           1 => H           4 => H           1 => T           4 => shift($(3/2))           4 => X           1 => H           4 => S           4 => X           2 => S           @ctrl 2 3 => X           2 => H           @ctrl 2 3 => X           2 => T           3 => S           2 => H           3 => H           3 => S           @ctrl 2 3 => X       enddemo_circ_simp (generic circuit with 1 methods)

Можно добавить аргумент optimizer = [opts...] в макрос @device, чтобы упростить эту схему во время компиляции. В настоящее время существует только два этапа оптимизации: :zx_clifford для упрощения Клиффорда [1] и :zx_teleport для фазовой телепортации [2]. Например, с помощью optimizer = [:zx_teleport] компилятор вызовет алгоритм фазовой телепортации в ZXCalculus.jl для упрощения схемы.

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

julia> using YaoLang.Compilerjulia> gate_count(demo_circ_simp)Dict{Any,Any} with 8 entries:  "YaoLang.Rx(3.141592653589793)"     => 2  "YaoLang.H"                         => 6  "YaoLang.Rx(0.7853981633974483)"    => 1  "YaoLang.shift(4.71238898038469)"   => 1  "YaoLang.shift(1.5707963267948966)" => 4  "YaoLang.shift(0.7853981633974483)" => 1  "@ctrl YaoLang.Z"                   => 1  "@ctrl YaoLang.X"                   => 4

Мы можем использовать YaoArrayRegister.jl чтобы применить эту упрощенную схему к квантовому состоянию:

julia> using YaoArrayRegister;julia> circ_teleport = demo_circ_simp()demo_circ_simp (quantum circuit)julia> r = rand_state(4);julia> r |> circ_teleportArrayReg{1, Complex{Float64}, Array...}    active qubits: 4/4

Можно также загрузить схемы из OpenQASM. OpenQASM это квантовая инструкция. Коды OpenQASM можно запускать на квантовых устройствах IBM. А квантовые схемы могут храниться в виде кодов OpenQASM. Я использовал пакет Джулии RBNF.jl (парсит код в ограниченную форму Бэкуса-Наура) для перевода OpenQASM кодов в AST, с последующим преобразованием его в YaoIR, промежуточное представление для смешанной квантово-классической программы YaoLang.jl. Это делает возможным чтение схем из OpenQASM в ZXCalculus.jl через YaoLang.jl.

using YaoLang: YaoIR, is_pure_quantumusing ZXCalculuslines = readlines("gf2^8_mult.qasm")src = prod([lines[1]; lines[3:end]])ir = YaoIR(@__MODULE__, src, :qasm_circ)ir.pure_quantum = is_pure_quantum(ir)circ = ZXDiagram(ir)pt_circ = phase_teleportation(circ)

Здесь мы загрузили схему в виде ZXDiagram из файла .qasm, который можно найти тут. И мы использовали алгоритм фазовой телепортации, чтобы упростить его. Мы видим, что Т-число в цепи уменьшилось с 448 до 264.

julia> tcount(circ)448julia> tcount(pt_circ)264

Приведенные выше примеры показали, как ZXCalculus.jl работает в качестве механизма упрощения схем в YaoLang.jl. А теперь давайте посмотрим, что скрывается за сценой.

Низкоуровневые интерфейсы ZXCalculus

В ZX-исчислении мы будем иметь дело с ZX-диаграммами и мультиграфами с некоторой дополнительной информацией. Каждая вершина ZX-диаграммы называется пауком. Есть два типа пауков, Z-паук и X-паук. Каждый паук связан с числом, называемым фазой. По нотации Дирака Z-паук и X-паук представляют собой следующие матрицы ранга 2.

ZX-диаграммы можно рассматривать как особый тип тензорной сети. Как и квантовые схемы. А квантовые схемы могут быть преобразованы в ZX-диаграммы по следующим правилам.

Желтая коробка, H-box, это простая нотация следующих пауков в ZX-исчислении

Для представления общих ZX-диаграмм мы определили структуру ZXDiagram в ZXCalculus.jl. Мы можем построить ZX-диаграмму, которая представляет собой пустую квантовую схему с n кубитами с помощью ZXDiagram(n). Например, если мы хотим упростить приведенную выше схему с помощью ZXCalculus.jl вручную, то сначала мы построим ZX-диаграмму схемы из 4 кубитов.

using ZXCalculuszxd = ZXDiagram(4)

Затем мы добавляем некоторые элементы к ZX-диаграмме, которую мы только что построили. Просто используем push_gate! и push_ctrl_gate!

push_gate!(zxd, Val(:Z), 1, 7//4)push_gate!(zxd, Val(:H), 1)push_gate!(zxd, Val(:X), 1, 1//4)push_gate!(zxd, Val(:H), 4)push_ctrl_gate!(zxd, Val(:CZ), 4, 1)push_ctrl_gate!(zxd, Val(:CNOT), 1, 4)push_gate!(zxd, Val(:H), 1)push_gate!(zxd, Val(:H), 4)push_gate!(zxd, Val(:Z), 1, 1//4)push_gate!(zxd, Val(:Z), 4, 3//2)push_gate!(zxd, Val(:X), 4, 1//1)push_gate!(zxd, Val(:H), 1)push_gate!(zxd, Val(:Z), 4, 1//2)push_gate!(zxd, Val(:X), 4, 1//1)push_gate!(zxd, Val(:Z), 2, 1//2)push_ctrl_gate!(zxd, Val(:CNOT), 3, 2)push_gate!(zxd, Val(:H), 2)push_ctrl_gate!(zxd, Val(:CNOT), 3, 2)push_gate!(zxd, Val(:Z), 2, 1//4)push_gate!(zxd, Val(:Z), 3, 1//2)push_gate!(zxd, Val(:H), 2)push_gate!(zxd, Val(:H), 3)push_gate!(zxd, Val(:Z), 3, 1//2)push_ctrl_gate!(zxd, Val(:CNOT), 3, 2)

Теперь давайте нарисуем ZX-диаграмму, которую мы построили. Инструмент визуализации ZXCalculus.jl в настоящее время предоставляется в YaoPlots.jl.

using YaoPlotsplot(zxd)

Заодно поднимем алгоритмы упрощения clifford_simplification [1] и phase_teleportation [2]

ex_zxd = clifford_simplification(zxd);pt_zxd = phase_teleportation(zxd);plot(ex_zxd)plot(pt_zxd)

Схема после упрощения Клиффорда

Схема после фазовой телепортации

Алгоритм фазовой телепортации может уменьшить число Т-вентилей квантовой схемы. Воспользуемся tcount чтобы показать количество Т-вентилей. В этом примере число Т уменьшилось с 4 до 2.

julia> tcount(zxd)4julia> tcount(pt_zxd)2

Эти алгоритмы используют правила ZX-исчисления для упрощения ZX-диаграмм. Эти правила определяют, как ZX-диаграммы могут быть преобразованы. Вот некоторые основные правила для ZXDiagram.

В работе [1] определен особый тип ZX-диаграммы графоподобная ZX-диаграмма. Мы используем ZXGraph, чтобы представить ее в ZXCalculus.jl. И вот некоторые правила для ZXGraph

Можно применить правила на ZX-диаграмме вручную. Для этого мы предоставляем различные API.

Функция match найдет все вершины соответствующие заданному правилу. Применение rewrite! переписывает ZX-диаграмму на некоторых совпадающих вершинах. replace! функция просто сопоставляет и переписывает все совпадающие вершины один раз. simplify! функция будет сопоставлять и переписывать ZX-диаграмму с правилом до тех пор, пока не останется сопоставимых вершин.

В clifford_simplification, мы сначала преобразуем заданную ZX-диаграмму в графоподобную ZX-диаграмму

zxg = ZXGraph(zxd)plot(zxg)

Затем, мы упрощаем ее с помощью правил :lc, :p1, и :pab

simplify!(Rule{:lc}(), zxg)simplify!(Rule{:p1}(), zxg)replace!(Rule{:pab}(), zxg)plot(zxg)

Наконец, мы извлекаем новую схему из упрощенной графоподобной ZX-диаграммы

ex_circ = circuit_extraction(zxg)plot(ex_circ)

Почему именно ZXCalculus.jl?

Вышеприведенные алгоритмы уже реализованы в пакете Python PyZX. PyZX это полнофункциональная библиотека для манипулирования крупномасштабными квантовыми схемами и ZX-диаграммами. Он предоставляет множество удивительных возможностей визуализации и поддерживает различные формы квантовых схем, включая QASM, Quipper и Quantomatic.

Так почему же мы разработали ZXCalculus.jl? Это связано с тем, что ZXCalculus.jl является не только полнофункциональной библиотекой для ZX-исчисления, но и одним из механизмов упрощения схем для YaoLang.jl. Легкая нативная реализация ZX-исчисления необходима, так как зависимость от пакета Python сделает работу компилятора медленнее и сложнее.

Benchmarks

Мы провели сравнительный анализ алгоритма фазовой телепортации на 40 схемах с различным числом вентилей (от 57 до 91642). ZXCalculus.jl имеет ускорение от 6 до 50 раз в этих примерах (время выполнения ZXCalculus.jl масштабируется до 1 для каждой схемы на этом рисунке). Эти тесты выполняются на ноутбуке с процессором Intel i7-10710U и 16 ГБ оперативной памяти. Код для бенчмарков можно найти здесь.

В большинстве примеров T-число оптимизированных схем, производимых ZXCalculus.jl, совпадает с PyZX. Однако в 2 примерах ZXCalculus.jl имеет больше T-count, чем PyZX. Это может быть вызвано различными стратегиями упрощения между ZXCalculus.jl и PyZX. Мы продолжим исследовать это в будущем.

Кроме того, YaoLang.jl поддерживает гибридные квантово-классические программы. Можно оптимизировать гибридные квантово-классические программы с помощью ZXCalculus.jl.

Резюме и дальнейшей работы

В течение GSoC 2020 я в основном выполнил следующие работы.

• Представление и манипулирование ZX-диаграммами с высокой производительностью.
• Реализация двух алгоритмов упрощения на основе ZX-исчисления.
• Добавление визуализации в ZX-схемы в YaoPlots.jl.
• Интеграция ZXCalculus.jl с YaoLang.jl.
• Добавление поддержки OpenQASM в YaoLang.jl.

Есть еще кое-что, что нужно отполировать.

• Поиск лучшей стратегии упрощения, чтобы получить более низкие значения T.
• Полная поддержка визуализации "ZXGraph" (сценарий построения графика может потерпеть неудачу на некоторых "ZXGraph" с фазовыми гаджетами).
• Преобразование на ZX-схемы для тензорных сетей без YaoLang.jl.
• Преобразование между "YaoIR" и "ZXDiagram" может привести к тому, что схема отличается от глобальной фазы. Мы должны записать эту глобальную фазу в более поздней версии.

Кроме того, я буду продолжать работать над YaoLang.jl с Роджером Ло, Чтобы поддерживать больше методов упрощения схем (методы сопоставления шаблонов, методы на основе Quon и т.д.).

Благодарности

Я хочу поблагодарить своих наставников, Роджера Ло и Цзиньго Лю. Без их помощи я не смог бы осуществить этот проект. ZXCalculus.jl очень вдохновлен PyZX. Спасибо Алекс Киссинджер и Джон ван де Ветеринг, авторам PyZX. Они дали мне полезные советы по алгоритму фазовой телепортации и рассмотрели бенчмарки между PyZX и ZXCalculus.jl. Благодарю Google за проведение Google Summer of Code, которое способствует развитию сообщества с открытым исходным кодом.

[1]: Graph-theoretic Simplification of Quantum Circuits with the ZX-calculus
[2]: Reducing T-count with the ZX-calculus

скришнот из 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
• Рисовалка надписи на гитхабе

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

Введение

Бенчмарки языков программирования на задаче вычисления ядра матрицы

$ ldc2 -O --boundscheck=off --ffast-math --mcpu=native --boundscheck=off mathdemo.d && ./mathdemoTime taken for c log: 0.324789 seconds.Time taken for d log: 2.30737 seconds.

Условия тестирования

$ julia --versionjulia version 1.4.1$ dmd --versionDMD64 D Compiler v2.090.1$ ldc2 --versionLDC - the LLVM D compiler (1.18.0):  based on DMD v2.088.1 and LLVM 9.0.0$ chpl --versionchpl version 1.22.0

Компиляция

chpl script.chpl kernelmatrix.chpl --fast && ./script

ldc2 script.d kernelmatrix.d arrays.d -O5 --boundscheck=off --ffast-math -mcpu=native && ./script

julia script.jl

Реализации

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

Chapel

proc calculateKernelMatrix(K, data: [?D] ?T){  var n = D.dim(0).last;  var p = D.dim(1).last;  var E: domain(2) = {D.dim(0), D.dim(0)};  var mat: [E] T;  var rowPointers: [1..n] c_ptr(T) =    forall i in 1..n do c_ptrTo(data[i, 1]);  forall j in guided(1..n by -1) {    for i in j..n {      mat[i, j] = K.kernel(rowPointers[i], rowPointers[j], p);      mat[j, i] = mat[i, j];    }  }  return mat;}

D

auto calculateKernelMatrix(alias K, T)(K!(T) kernel, Matrix!(T) data){  long n = data.ncol;  auto mat = Matrix!(T)(n, n);  foreach(j; taskPool.parallel(iota(n)))  {    auto arrj = data.refColumnSelect(j).array;    foreach(long i; j..n)    {      mat[i, j] = kernel(data.refColumnSelect(i).array, arrj);      mat[j, i] = mat[i, j];    }  }  return mat;}

Julia

function calculateKernelMatrix(Kernel::K, data::Array{T}) where {K <: AbstractKernel,T <: AbstractFloat}  n = size(data)[2]  mat = zeros(T, n, n)  @threads for j in 1:n      @views for i in j:n          mat[i,j] = kernel(Kernel, data[:, i], data[:, j])      end  end  return Symmetric(mat, :L)end

struct DotProduct <: AbstractKernel end@inline function kernel(K::DotProduct, x::AbstractArray{T, N}, y::AbstractArray{T, N}) where {T,N}  ret = zero(T)  m = length(x)  @inbounds @simd for k in 1:m      ret += x[k] * y[k]  end  return retend

Использование памяти

Command being timed: "./script"User time (seconds): 113190.32System time (seconds): 6.57Percent of CPU this job got: 1196%Elapsed (wall clock) time (h:mm:ss or m:ss): 2:37:39Average shared text size (kbytes): 0Average unshared data size (kbytes): 0Average stack size (kbytes): 0Average total size (kbytes): 0Maximum resident set size (kbytes): 5761116Average resident set size (kbytes): 0Major (requiring I/O) page faults: 0Minor (reclaiming a frame) page faults: 1439306Voluntary context switches: 653Involuntary context switches: 1374820Swaps: 0File system inputs: 0File system outputs: 8Socket messages sent: 0Socket messages received: 0Signals delivered: 0Page size (bytes): 4096Exit status: 0

Command being timed: "./script"User time (seconds): 106065.71System time (seconds): 58.56Percent of CPU this job got: 1191%Elapsed (wall clock) time (h:mm:ss or m:ss): 2:28:29Average shared text size (kbytes): 0Average unshared data size (kbytes): 0Average stack size (kbytes): 0Average total size (kbytes): 0Maximum resident set size (kbytes): 20578840Average resident set size (kbytes): 0Major (requiring I/O) page faults: 0Minor (reclaiming a frame) page faults: 18249033Voluntary context switches: 3833Involuntary context switches: 1782832Swaps: 0File system inputs: 0File system outputs: 8Socket messages sent: 0Socket messages received: 0Signals delivered: 0Page size (bytes): 4096Exit status: 0

Command being timed: "julia script.jl"User time (seconds): 49794.85System time (seconds): 30.58Percent of CPU this job got: 726%Elapsed (wall clock) time (h:mm:ss or m:ss): 1:54:18Average shared text size (kbytes): 0Average unshared data size (kbytes): 0Average stack size (kbytes): 0Average total size (kbytes): 0Maximum resident set size (kbytes): 7496184Average resident set size (kbytes): 0Major (requiring I/O) page faults: 794Minor (reclaiming a frame) page faults: 38019472Voluntary context switches: 2629Involuntary context switches: 523063Swaps: 0File system inputs: 368360File system outputs: 8Socket messages sent: 0Socket messages received: 0Signals delivered: 0Page size (bytes): 4096Exit status: 0

Оптимизация производительности

• Статическая диспетчеризация функций ядра вместо использования полиморфизма. Это означает, что при передаче функции ядра используется параметрический (времени компиляции) полиморфизм, а не динамический (времени исполнения), при котором диспетчеризация с виртуальными функциями влечет за собой накладные расходы.
• Использование представлений/ссылок, вместо копирования данных в многопоточном режиме, имеет большое значение.
• Распараллеливание вычислений имеет огромное значение.
• Знание того, что массив является основным для строки/столбца, и использование этого в вычислениях имеет огромное значение.
• Проверки границ и оптимизации компилятора дают огромную разницу, особенно в Chapel и D.
• Включение SIMD в D и Julia внесло свой вклад в производительность. В D это было сделано с помощью флага -mcpu=native, а в Julia это было сделано с помощью макроса @ simd.

Качество жизни

Интерактивность

real0m0.545suser0m0.447ssys0m0.101s

real0m5.980suser0m5.787ssys0m0.206s

Документация и примеры

Поддержка многомерных массивов

Мощность языка

Конкурентность и параллельное программирование

Стандартная библиотека

Менеджеры и экосистема пакетов

Интеграция с Cи

Сообщество

Резюме