28.11.2018 18:04
Yermack 3 1200 Источник


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


Строки


Создаются строковые переменные с помощью обрамления символов в двойные кавычки, одинарными же задаются одиночные символы типа char. Конкатенация строк выполняется с помощью умножения " * ":


cats = 4
s1 = "How many cats ";
s2 = "is too many cats?";
s3 = s1*s2

Out[]: "How many cats is too many cats?"

То есть возведение в степень действует так:


s1^3
Out[]: "How many cats How many cats How many cats "

С другой стороны, на строки распространяется индексация:


s1[3]
Out[]: 'w': ASCII/Unicode U+0077 (category Ll: Letter, lowercase)

s1[5:13]
Out[]: "many cats"

s1[13:-1:5]
Out[]: "stac ynam"

s2[end]
Out[]: '?': ASCII/Unicode U+003f (category Po: Punctuation, other)

Строки можно создавать из строк и других типов:


s4 = string(s3, " - I don't know, but ", cats, " is too few.")
# или так s4 = "$s3 - I don't know, but $cats is too few."

Out[]: "How many cats is too many cats? - I don't know, but 4 is too few."

В арсенале есть много всяких полезных функций, например поиск элемента:


findfirst( isequal('o'), s4 )
Out[]: 4

findlast( isequal('o'), s4 )
Out[]: 26

findnext( isequal('o'), s4, 7 ) # ищем совпадение после седьмого элемента
Out[]: 11

А вот так можно запросто реализовать Шифр Цезаря


caesar(X, n) = prod( [x += n for x in X] )

str3 = "хатифнатты ведут порочный образ жизни";

caesar(str3, 3)
Out[]: "шгхлчргххю#еизцх#тсусърюм#сдугк#йлкрл"

str4 = caesar(str3, -32)
Out[]: "ХАТИФНАТТЫ\0ВЕДУТ\0ПОРОЧНЫЙ\0ОБРАЗ\0ЖИЗНИ"

"Latin letters go before Cyrillic" < str4 < str3
Out[]: true

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


Символьные вычисления


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



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


m1 = [1 1 "a"; 1 0 1]
Out[]: 2?3 Array{Any,2}:
            1  1   "a"
            1  0  1 

m3 = [1 2 "ln(3)"; 2 1 0]
Out[]: 2?3 Array{Any,2}:
            1  2   "ln(3)"
            2  1  0

m1+m3
Out[]: MethodError: no method matching +(::String, ::String) ...

Для начала, запустим свои шаловливые ручки в базовые операторы:


import Base: *, -, +

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


+(a::String, b::String) = a * "+" * b
Out[]: + (generic function with 164 methods)

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


m1+m3
Out[]: 2?3 Array{Any,2}:
            2  3   "a+ln(3)"
            3  1  1

Теперь подпортим умножение, а именно определим случай умножения строки на число:


function *(a::String, b::Number) 
    if b == 0
        0
    else
        if b == 1
            a
        else
            "$b" * "*(" * a * ")"
        end
    end
end
Out[]: * (generic function with 344 methods)

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


*(a::String, b::Number) = (b==0) ? 0 : (b==1) ? a : "$b" * "*(" * a * ")"
Out[]:

Подвергнем подобным экзекуциям остальные случаи и получим что-то вроде:


import Base: *, -, +
+(a::String, b::String) = a * "+" * b
*(a::String, b::String) = a * "*(" * b * ")"
*(a::String, b::Number) = (b==0) ? 0 : (b==1) ? a : "$b" * "*(" * a * ")"
*(b::Number, a::String) = (b==0) ? 0 : (b==1) ? a : "$b" * "*(" * a * ")"
+(a::String, b::Number) = (b==0) ? a : a * "+" * "$b"
+(b::Number, a::String) = (b==0) ? a : a * "+" * "$b"
-(b::Number, a::String) = (b==0) ? "-" *  a : "$b" * "-" * a
-(a::String, b::Number) = (b==0) ? a : a * "-" * "$b"

# m1 = [1 1 "a"; 1 0 1]
m2 = [2 0 2; 2 "b" 2; 2 2 2]

m1*m2
Out[]: 2?3 Array{Any,2}:
           "2*(a)+4"   "b+2*(a)"   "2*(a)+4"
            4           2           4

А посчитаем-ка детерминант! Но так как встроенная функция устроена сложно, наших методов недостаточно, чтоб скормить ей матрицу со строками — получим ошибку. Значит реализуем на скорую руку свой с использованием символа Леви-Чивиты


? = zeros(Int, 3,3,3)
?[1,2,3] = ?[2,3,1] = ?[3,1,2] = 1
?[3,2,1] = ?[1,3,2] = ?[2,1,3] = -1
?
Out[]: 3?3?3 Array{Int64,3}:
[:, :, 1] =
 0   0  0
 0   0  1
 0  -1  0

[:, :, 2] =
 0  0  -1
 0  0   0
 1  0   0

[:, :, 3] =
  0  1  0
 -1  0  0
  0  0  0

Формулу смешанного произведения


$ \left[\vec{a}\vec{b}\vec{c}\right]=\sum_{i,j,k=1}^{3}\varepsilon_{ijk}a^ib^jc^k $


можно использовать для нахождения определителя матрицы 3х3, понимая a, b, c как первую, вторую и третью строку (столбец) соответственно


detrmnt(arr) = sum( ?[i,j,k]*arr[1,i]*arr[2,j]*arr[3,k] for i in 1:3, j in 1:3, k in 1:3 )

detrmnt(["a" 2 "b"; 1 0 1; "c" 2 "d"])
Out[]: "2*(c)+2*(b)+2*(-1*(a))+-2*(d)"

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


Rx = [1 0 0 0; 0 "cos(?)" "sin(?)" 0; 0 "-sin(?)" "cos(?)" 0; 0 0 0 1];
Rz = ["cos(?)" "sin(?)" 0 0; "-sin(?)" "cos(?)" 0 0; 0 0 1 0; 0 0 0 1];
T = [1 0 0 0; 0 1 0 0; 0 0 1 0; "x0" "y0" "z0" 1];

Rx*Rz
Out[]: MethodError: no method matching zero(::String) ...

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


function *( a::Array{Any,2}, b::Array{Any,2} )
    if size(a)[2] == size(b)[1]
        res = Array{Any}(undef, size(a)[1], size(b)[2] )
        fill!(res, 0)

        for i = 1:size(a)[1], j = 1:size(b)[2], k = 1:size(a)[2]
            res[i,j] = res[i,j] + a[i,k]*b[k,j]
        end
        res
    else
        error("Matrices must have same size mult")
    end
end
Out[]: * (generic function with 379 methods)

Теперь можно смело умножать:


X = Rx*Rz*T
Out[]: 4?4 Array{Any,2}:
 "cos(?)"             "sin(?)"            0          0
 "cos(?)*(-sin(?))"   "cos(?)*(cos(?))"    "sin(?)"  0
 "-sin(?)*(-sin(?))"  "-sin(?)*(cos(?))"   "cos(?)"  0
 "x0"                 "y0"                 "z0"      1

Попридержим эту матрицу в уме, а пока начнем первые шаги в


Метапрограммирование


Котировки


Джулия поддерживает метапрограммирование. Это похоже на символьное программирование, где мы имеем дело с выражениями (например, 6 * 7) в отличие от значений (например, 42). Загоняя операции в интерпретатор, мы сразу получаем результат, что, как мы убедились, можно обойти с помощью строк:


x = "6*7"
Out[]: "6*7"

Строки можно превратить в выражения с помощью parse(), а потом вычислить, используя eval():


eval(Meta.parse(ans))
Out[]: 42

К чему все эти сложности? Фокус в том, что, когда мы имеем выражение, мы можем его модифицировать различными интересными способами:


x = replace(x, "*" => "+")
eval(Meta.parse(x))
Out[]: 13

Чтоб избежать возни со строками предусмотрен оператор хмурого лица :()


y = :(2^8-1)
Out[]: :(2^8-1)

eval(y)
Out[]: 255

Можно "зацитировать" выражение, функцию, блок кода...


quote
  x = 2 + 2
  hypot(x, 5)
end

Out[]: quote
    #= In[13]:2 =#
    x = 2 + 2
    #= In[13]:3 =#
    hypot(x, 5)
end

:(function mysum(xs)
    sum = 0
    for x in xs
      sum += x
    end
  end)

Out[]: :(function mysum(xs)
      #= In[14]:2 =#
      sum = 0
      #= In[14]:3 =#
      for x = xs
          #= In[14]:4 =#
          sum += x
      end
  end)

А теперь отпарсим посчитанную ранее матрицу (лучше начать новый сеанс после предыдущего раздела, мои некрасивые перегрузки запросто могут конфликтовать с подключаемыми пакетами):


X = [ "cos(?)"             "sin(?)"            0          0;
 "cos(?)*(-sin(?))"   "cos(?)*(cos(?))"    "sin(?)"  0;
 "-sin(?)*(-sin(?))"  "-sin(?)*(cos(?))"   "cos(?)"  0;
 "x0"                 "y0"                 "z0"      1; ]

for i = 1:size(X,1), j = 1:size(X,2)
    if typeof(X[i,j]) == String
        X[i,j] = Meta.parse(X[i,j]) 
    end
end
X

Out[]: 4?4 Array{Any,2}:
 :(cos(?))                 :(sin(?))              0           0
 :(cos(?) * -(sin(?)))     :(cos(?) * cos(?))      :(sin(?))  0
 :(-(sin(?)) * -(sin(?)))  :(-(sin(?)) * cos(?))   :(cos(?))  0
 :x0                       :y0                     :z0        1

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


? = 20*pi/180
x0 = 4
y0 = -0.5
z0 = 1.3

Xtr = [ eval(x) for x in X]

Out[]: 4?4 Array{Real,2}:
  0.939693   0.34202   0         0
 -0.321394   0.883022  0.34202   0
  0.116978  -0.321394  0.939693  0
  4         -0.5       1.3       1

Эта матрица производит поворот вокруг осей X и Z на $20^{\circ}$ и перенос на $\vec{R}=\{x0,y0,z0\}$


x = [-1 -1 1 1 -1 -1 1 1 -1 -1]; 
y = [-1 -1 -1 -1 -1 1 1 1 1 -1];
z = [1 -1 -1 1 1 1 1 -1 -1 -1]
R = [ x' y' z' ones( length(x) ) ]
plot(x', y', z', w = 3)


R2 = R*Xtr
plot!( R2[:,1], R2[:,2], R2[:,3], w = 3)


Fruit of the Expression Tree


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


x = "миллион"
print("$x $x $x алых роз...")

Out[]: миллион миллион миллион алых роз...

С котировками (цитатами) — та же история:


x = :(6*7)
y = :($x + $x)
Out[]: :(6 * 7 + 6 * 7)

eval(y)
Out[]: 84

The Root of all Eval


eval() не только возвращает значение выражения. Попробуем зацитировать объявление функции:


ex = :( cats() = println("Meow!") )
cats()

Out[]: UndefVarError: cats not defined

А теперь оживим её:


eval(ex)
Out[]: cats (generic function with 1 method)

cats()
Out[]: Meow!

Используя интерполяцию, мы можем построить определение функции «на лету»; на самом деле, мы можем сразу сделать целый ряд функций.


for name in [:dog, :bird, :mushroom]
  println(:($name() = println($("I'm $(name)!"))))
end

Out[]: dog() = begin
        #= In[27]:2 =#
        println("I'm dog!")
    end
bird() = begin
        #= In[27]:2 =#
        println("I'm bird!")
    end
mushroom() = begin
        #= In[27]:2 =#
        println("I'm mushroom!")
    end

for name in [:dog, :bird, :mushroom]
  eval(:($name() = println($("I'm $(name)!"))))
end

dog()
Out[]: I'm dog!

mushroom()
Out[]: I'm mushroom!

Это может оказаться чрезвычайно полезным при обертке API (скажем, из библиотеки C или через HTTP). API-интерфейсы часто определяют список доступных функций, поэтому вы можете захватить их и создать всю оболочку сразу! См. Примеры Clang.jl, TensorFlow.jl или базовая линейная алгебра.


Original sin


Вот более практичный пример. Рассмотрим следующее определение функции sin(), основанное на рядах Тейлора:


$ sin(x) = \sum_{k=1}^{\infty} \frac{(-1)^k}{(1+2k)!} x^{1+2k} $


mysin(x) = sum((-1)^k/factorial(1+2*k) * x^(1+2k) for k = 0:5)
mysin(0.5), sin(0.5)

Out[]: (0.4794255386041834, 0.479425538604203)

using BenchmarkTools
@benchmark mysin(0.5)

Out[]: BenchmarkTools.Trial: 
  memory estimate:  112 bytes
  allocs estimate:  6
  --------------
  minimum time:     1.105 ?s (0.00% GC)
  median time:      1.224 ?s (0.00% GC)
  mean time:        1.302 ?s (0.00% GC)
  maximum time:     9.473 ?s (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     10

Сейчас это намного медленнее, чем могло бы быть. Причина в том, что мы перебираем k, что относительно дорого. Явная запись гораздо быстрее:


mysin(x) = x - x^3/6 + x^5/120 # + ...

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


plus(a, b) = :($a + $b)
plus(1, 2)

Out[]: :(1 + 2)

С plus() мы можем делать много интересных вещей, например, символьная сумма:


reduce(+, 1:10)
Out[]: 55

reduce(plus, 1:10)
Out[]: :(((((((((1 + 2) + 3) + 4) + 5) + 6) + 7) + 8) + 9) + 10)

eval(ans)
Out[]: 55

reduce(plus, [:(x^2), :x, 1])
Out[]: :((x ^ 2 + x) + 1)

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


k = 2
:($((-1)^k) * x^$(1+2k) / $(factorial(1+2k)))
Out[]: :((1 * x ^ 5) / 120)

Теперь можно клепать элементы ряда как на конвейере:


terms = [:($((-1)^k) * x^$(1+2k) / $(factorial(1+2k))) for k = 0:5]
Out[]: 6-element Array{Expr,1}:
 :((1 * x ^ 1) / 1)         
 :((-1 * x ^ 3) / 6)        
 :((1 * x ^ 5) / 120)       
 :((-1 * x ^ 7) / 5040)     
 :((1 * x ^ 9) / 362880)    
 :((-1 * x ^ 11) / 39916800)

И суммировать их:


reduce(plus, ans)
Out[]: :((((((1 * x ^ 1) / 1 + (-1 * x ^ 3) / 6) + (1 * x ^ 5) / 120) + (-1 * x ^ 7) / 5040) + (1 * x ^ 9) / 362880) + (-1 * x ^ 11) / 39916800)

:(mysin(x) = $ans)

Out[]: :(mysin(x) = begin
          #= In[52]:1 =#
          (((((1 * x ^ 1) / 1 + (-1 * x ^ 3) / 6) + (1 * x ^ 5) / 120) + (-1 * x ^ 7) / 5040) + (1 * x ^ 9) / 362880) + (-1 * x ^ 11) / 39916800
      end)

eval(ans)
mysin(0.5), sin(0.5)
Out[]: (0.4794255386041834, 0.479425538604203)

@benchmark mysin(0.5)

Out[]: BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     414.363 ns (0.00% GC)
  median time:      416.328 ns (0.00% GC)
  mean time:        431.885 ns (0.00% GC)
  maximum time:     3.352 ?s (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     201

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

Комментарии (3)


  1. kirtsar
    28.11.2018 23:17
    #19429560

    Про символьные вычисления можно ещё вспомнить Reduce.


  1. sairus777
    29.11.2018 03:31
    #19429964

    Хороший пример метапрограммирования — в пакете StructArrays.jl, где очень элегантно решена проблема SoA / AoS


  1. Yermack Автор
    29.11.2018 08:04
    #19430238

    Почему-то из вставки «Исходный код» пропала julia, и перестала работать подсветка, поэтому сейчас всё поменял на ```cpp```