text:
function allfeatures()
    Feature[
        @feature((scrabble_score(word) == j for j in 1:26), "has scrabble score $j")
        @feature((letter_tallies(word)[c - 'a' + 1] == j for j in 1:5, c in 'a':'z'), "contains $j '$c's")
        @feature((letter_tallies(word)[c - 'a' + 1] >= j for j in 1:5, c in 'a':'z'), "contains at least $j '$c's")
        @feature((length(word) >= j && word[j] == c for c in 'a':'z', j in 1:12), "contains '$c' at index $j")
        @feature((length(word) >= j && word[end - j + 1] == c for c in 'a':'z', j in 1:6), "contains '$c' at index $j from end")
        @feature((num_unique_vowels(letter_tallies(word)) == j for j in 1:5), "has $j unique vowels")
        @feature((num_unique_consonants(letter_tallies(word)) == j for j in 1:15), "has $j unique consonants")
        @feature((num_unique_letters(letter_tallies(word)) == j for j in 1:26), "has $j unique letters")
        @feature((vowel_pattern(word, (true, false)) || vowel_pattern(word, (false, true))), "alternates consonant vowel")
        @feature((vowel_pattern(word, p) for p in ((true, false),
                                                   (false, true))), "has vowel/consonant pattern $p")
        @feature((num_double_letters(word) == j for j in 1:3), "contains $j double letters")
        @feature((num_repeats(letter_tallies(word), i) == j for j in 1:5, i in 2:4), "has $j letters repeated at least $i times each")
        @feature((num_repeats_strict(letter_tallies(word), i) == j for j in 1:5, i in 2:4), "has $j letters repeated exactly $i times each")
        @feature((num_repeats(letter_tallies(word), i) >= j for j in 1:5, i in 2:4), "has at least $j letters repeated at least $i times each")
        @feature((num_repeats_strict(letter_tallies(word), i) >= j for j in 1:5, i in 2:4), "has at least $j letters repeated exactly $i times each")
        @feature(contains_repeated_consonant(letter_tallies(word)), "contains a repeated consonant")
        @feature(contains_repeated_vowel(letter_tallies(word)), "contains a repeated vowel")
        @feature(contains_day_of_week(word), "contains a day of the week abbreviation")
        @feature(is_hill(word), "is a hill word")
        @feature(is_valley(word), "is a valley word")
        @feature(word == reverse(word), "is a palindrome")
        @feature(is_pyramid(letter_tallies(word)), "is a pyramid word")
        @feature((num_alpha_bigrams(word) == j for j in 0:10), "has $j alphabetical bigrams")
        @feature((num_alpha_bigrams(word) >= j for j in 1:10), "has at least $j alphabetical bigrams")
        @feature((num_reverse_alpha_bigrams(word) == j for j in 0:10), "has $j reverse alphabetical bigrams")
        @feature((num_reverse_alpha_bigrams(word) >= j for j in 1:10), "has at least $j reverse alphabetical bigrams")
        @feature((num_sequential_bigrams(word) == j for j in 0:10), "has $j sequential bigrams")
        @feature((num_sequential_bigrams(word) >= j for j in 1:10), "has at least $j sequential bigrams")
        @feature((num_reverse_sequential_bigrams(word) == j for j in 0:10), "has $j reverse sequential bigrams")
        @feature((num_reverse_sequential_bigrams(word) >= j for j in 1:10), "has at least $j reverse sequential bigrams")
        @feature((num_cardinal_directions(word) == j for j in 1:6), "has $j cardinal direction abbreviations (NESW)")
        @feature((num_cardinal_directions(word) >= j for j in 1:6), "has at least $j cardinal direction abbreviations (NESW)")
        @feature(ismatch(GREEK_REGEX, word), "contains a greek letter")
        @feature((longest_vowel_streak(word) == j for j in 2:5), "has $j vowels in a row")
        @feature((longest_vowel_streak(word) >= j for j in 2:5), "has at least $j vowels in a row")
        @feature((longest_consonant_streak(word) == j for j in 2:5), "has $j consonants in a row")
        @feature((longest_consonant_streak(word) >= j for j in 2:5), "has at least $j consonants in a row")
        @feature((num_morse_bits(letter_tallies(word)) == j for j in 1:10), "has $j morse bits ('i's and 't's)")
        @feature((num_morse_bits(letter_tallies(word)) >= j for j in 1:5), "has at least $j morse bits ('i's and 't's)")
        @feature(ismatch(ENTIRELY_ELEMENTS_REGEX, word), "can be completely broken down into chemical element symbols")
        @feature(ismatch(ENTIRELY_STATES_REGEX, word), "can be completely broken down into US state abbreviations")
        @feature((num_state_abbreviations(word) == j for j in 1:5), "contains $j US state abbreviations")
        @feature(in(word, MA_BELL_EXCHANGES_SET), "is a Ma Bell recommended telephone exchange name")
        @feature((has_transaddition(BitsTally(word), c) for c in 'a':'z'), "has a transaddition with letter '$c'")
        @feature((has_transdeletion(BitsTally(word), c) for c in 'a':'z'), "has a transdeletion with letter '$c'")
        @feature(has_transaddition(BitsTally(word)), "has a 1-letter transaddition")
        @feature(has_transdeletion(BitsTally(word)), "has a 1-letter transdeletion")
        @feature(in(word, UNICODE_CHARACTERS), "is the name of a unicode character")
        ]
end

function scrabble_score(word::String)
    score = 0
    for c in word
        score += scrabble_score(c)
    end
    score
end

scrabble_score(char::Char) = SCRABBLE_SCORES[char]

const SCRABBLE_SCORES = Dict{Char, Int}(
    'e' => 1,
    'a' => 1,
    'i' => 1,
    'o' => 1,
    'n' => 1,
    'r' => 1,
    't' => 1,
    'l' => 1,
    's' => 1,
    'u' => 1,
    'd' => 2,
    'g' => 2,
    'b' => 3,
    'c' => 3,
    'm' => 3,
    'p' => 3,
    'f' => 4,
    'h' => 4,
    'v' => 4,
    'w' => 4,
    'y' => 4,
    'k' => 5,
    'j' => 8,
    'x' => 8,
    'q' => 10,
    'z' => 10
    )

const ALPHABET = 'a':'z'
const VOWELS = UInt8[c - 'a' + 1 for c in "aeiouy"]
const CONSONANTS = UInt8[c - 'a' + 1 for c in "bcdfghjklmnpqrstvwxyz"]
const VOWELS_SET = Set(ALPHABET[VOWELS])
const CONSONANTS_SET = Set(ALPHABET[CONSONANTS])
const ELEMENT_DATA = readdlm(joinpath(Pkg.dir("Collective"), "data", "elements.tsv"), '\t', String, skipstart=1)
const ELEMENTAL_SYMBOLS = lowercase.(strip.(ELEMENT_DATA[:,2]))
const STATES_DATA = readdlm(joinpath(Pkg.dir("Collective"), "data", "states.tsv"), '\t', String, skipstart=1)
const STATE_ABBREVIATIONS = strip.(lowercase.(STATES_DATA[:,2]))
const WORDS = wordlist(open(joinpath(Pkg.dir("Collective"), "data", "113809of.fic")))
const UNICODE_CHAR_DATA = readdlm(joinpath(Pkg.dir("Collective"), "data", "unicode_characters.tsv"), '\t', String, skipstart=1)
const UNICODE_CHARACTERS = Set(cleanup.(UNICODE_CHAR_DATA[:,1]))

const bitstallies = Set(BitsTally.(WORDS))

has_transaddition(t::BitsTally, c::Char) = (t + c) in bitstallies
has_transdeletion(t::BitsTally, c::Char) = (t - c) in bitstallies

function has_transaddition(t::BitsTally)
    for letter in 'a':'z'
        if has_transaddition(t, letter)
            return true
        end
    end
    return false
end

function has_transdeletion(t::BitsTally)
    for letter in 'a':'z'
        if has_transdeletion(t, letter)
            return true
        end
    end
    return false
end

has_anagram(t::BitsTally) = t in bitstallies

parenwrap(s) = "($s)"

const SINGLE_ELEMENT_REGEX = Regex("$(join((parenwrap(s) for s in ELEMENTAL_SYMBOLS), '|'))")
const ENTIRELY_ELEMENTS_REGEX = Regex("^($(join((parenwrap(s) for s in ELEMENTAL_SYMBOLS), '|')))*\$")

const SINGLE_STATE_REGEX = Regex("$(join((parenwrap(s) for s in STATE_ABBREVIATIONS), '|'))")
const ENTIRELY_STATES_REGEX = Regex("^($(join((parenwrap(s) for s in STATE_ABBREVIATIONS), '|')))*\$")

const GREEK_REGEX = r"(alpha)|(beta)|(gamma)|(delta)|(epsilon)|(zeta)|(eta)|(theta)|(iota)|(kappa)|(lambda)|(mu)|(nu)|(omicron)|(pi)|(rho)|(sigma)|(tau)|(upsilon)|(phi)|(chi)|(psi)|(omega)"

const MA_BELL_EXCHANGES_SET = Set(lowercase.(readdlm(joinpath(Pkg.dir("Collective"), "data", "ma_bell_exchanges.tsv"), '\t', String)))

isconsonant(char) = char in CONSONANTS_SET
isvowel(char) = char in VOWELS_SET

function num_unique_vowels(tallies)
    c = 0
    for x in VOWELS
        c += (tallies[x] > 0)
    end
    c
end

function num_unique_consonants(tallies)
    c = 0
    for x in CONSONANTS
        c += (tallies[x] > 0)
    end
    c
end

function num_unique_letters(tallies)
    c = 0
    for t in tallies
        c += (t > 0)
    end
    c
end

function vowel_pattern(word, pattern)
    i = 1
    for c in word
        if pattern[i]
            if !isvowel(c)
                return false
            end
        else
            if !isconsonant(c)
                return false
            end
        end
        if i == length(pattern)
            i = 1
        else
            i += 1
        end
    end
    true
end

function num_double_letters(word)
    c = 0
    for i in 1:(length(word) - 1)
        if word[i] == word[i+1]
            c += 1
        end
    end
    c
end

function letter_tallies(word)
    tallies = zeros(UInt8, 26)
    for c in word
        tallies[convert(UInt8, c - 'a' + 1)] += 1
    end
    tallies
end

function num_repeats_strict(tallies, n)
    c = 0
    for t in tallies
        c += (t == n)
    end
    c
end

function num_repeats(tallies, n)
    c = 0
    for t in tallies
        c += (t >= n)
    end
    c
end

function contains_repeated_consonant(tallies)
    for x in CONSONANTS
        if tallies[x] > 1
            return true
        end
    end
    false
end

function contains_repeated_vowel(tallies)
    for x in VOWELS
        if tallies[x] > 1
            return true
        end
    end
    false
end

contains_day_of_week(word) = ismatch(r"(sat)|(sun)|(mon)|(tue)|(wed)|(thu)|(fri)", word)

"""
Letters are alpha, then reverse alpha
"""
function is_hill(word)
    has_rise = false
    has_fall = false
    rising = true
    for i in 1:(length(word) - 1)
        if word[i + 1] > word[i]
            if !rising
                return false
            end
            has_rise = true
        elseif word[i + 1] < word[i]
            rising = false
            has_fall = true
        end
    end
    has_rise && has_fall
end

"""
Letters are reverse alpha, then alpha
"""
function is_valley(word)
    has_rise = false
    has_fall = false
    rising = false
    for i in 1:(length(word) - 1)
        if word[i + 1] < word[i]
            if rising
                return false
            end
            has_fall = true
        elseif word[i + 1] > word[i]
            rising = true
            has_rise = true
        end
    end
    has_rise && has_fall
end

function is_sequential(a::AbstractArray)
    for i in 1:(length(a) - 1)
        a[i+1] == a[i] + 1 || return false
    end
    true
end

"""
Letter tally is 1, 2, 3, etc.
"""
function is_pyramid(tallies)
    nonzero_tallies = Int[t for t in tallies if t > 0]
    is_sequential(sort!(nonzero_tallies))
end

function num_alpha_bigrams(word)
    count = 0
    for i in 1:(length(word) - 1)
        if word[i] < word[i + 1]
            count += 1
        end
    end
    count
end

function num_reverse_alpha_bigrams(word)
    count = 0
    for i in 1:(length(word) - 1)
        if word[i] > word[i + 1]
            count += 1
        end
    end
    count
end

function num_sequential_bigrams(word)
    count = 0
    for i in 1:(length(word) - 1)
        if word[i] + 1 == word[i+1]
            count += 1
        end
    end
    count
end

function num_reverse_sequential_bigrams(word)
    count = 0
    for i in 1:(length(word) - 1)
        if word[i] - 1 == word[i+1]
            count += 1
        end
    end
    count
end

function num_cardinal_directions(word)
    num_n = 0
    num_e = 0
    num_s = 0
    num_w = 0
    for c in word
        num_n += c == 'n'
        num_e += c == 'e'
        num_s += c == 's'
        num_w += c == 'w'
    end
    num_n + num_e + num_s + num_w
end

function longest_vowel_streak(word)
    current_streak = 0
    longest_streak = 0
    for c in word
        if isvowel(c)
            current_streak += 1
            if current_streak > longest_streak
                longest_streak = current_streak
            end
        else
            current_streak = 0
        end
    end
    longest_streak
end

function longest_consonant_streak(word)
    current_streak = 0
    longest_streak = 0
    for c in word
        if isconsonant(c)
            current_streak += 1
            if current_streak > longest_streak
                longest_streak = current_streak
            end
        else
            current_streak = 0
        end
    end
    longest_streak
end

num_morse_bits(tallies) = tallies['t' - 'a' + 1] + tallies['i' - 'a' + 1]

function num_elemental_symbols(word)
    i = 0
    foreach(m -> i += 1, eachmatch(SINGLE_ELEMENT_REGEX, word, true))
    i
end

function num_state_abbreviations(word)
    i = 0
    foreach(m -> i += 1, eachmatch(SINGLE_STATE_REGEX, word, true))
    i
end


text:
import Escher: @api
using Base.Test

comparefields(x, y) =
    try
       Bool[getfield(x, f) == getfield(y, f) for f in union(fieldnames(x), fieldnames(y))] |> all
    catch
        false
    end

@api test1 => TestType1 begin
    arg(a::AbstractFloat)
end
@test test1(3) == TestType1(3.0)

@api testsub => (TestSubType <: Number) begin
    arg(a::AbstractFloat)
end
@test testsub(3) == TestSubType(3.0)

@api test2 => TestType2 begin
    typedarg(a::AbstractFloat)
    arg(b::Int)
    curry(x::Int)
end

@test_throws MethodError test2(2, 3, 4) # Because 2 should be float
@test comparefields(test2(2.0, 3, 4), TestType2(2.0, 3, 4))
@test comparefields(test2(2.0, 3)(4), TestType2(2.0, 3, 4))

@api test3 => TestType3 begin
    typedarg(a::AbstractFloat)
    arg(b::Int)
    curry(x::Int)
    kwarg(p::Any="black")
    typedkwarg(q::AbstractFloat=1.0)
end

@test_throws MethodError test3(2, 3, 4) # Because 2 should be float
@test_throws TypeError test3(2.0, 5, 4, q=1) # q should be float

@test comparefields(test3(2.0, 3, 4), TestType3(2.0, 3, 4, "black", 1.0))
@test comparefields(test3(2.0, 3)(4), TestType3(2.0, 3, 4, "black", 1.0))
@test comparefields(test3(2.0, 3, 4, p=:white), TestType3(2.0, 3, 4, :white, 1))
@test comparefields(test3(2.0, 3, q=5.0)(4), TestType3(2.0, 3, 4, "black", 5.0))

abstract A
immutable B<:A end

#@show macroexpand(:(
@api test4 => (Test4{P <: A} <: Tile) begin
    typedarg(a::AbstractArray=Integer[])
    arg(b::P)
    curry(c::Int)
end
#))
@test comparefields(test4(B())(2), test4(Integer[], B(), 2))

ext: jl


text:
export Operator
module Operator

using CUDA
using GemmKernels
using GemmKernels.Tiling

# -------------------------------------
# Default definition for padded layouts
# -------------------------------------

for f in (:fragtype_a, :fragtype_b, :fragtype_accum, :load_a, :load_b, :load_c, :store_d)
    @eval @inline $f(op, ::Type{Layout.Padded{L, P}}, args...) where {L, P} = $f(op, L, args...)
end

# ----
# WMMA
# ----

struct WMMAOp{M, N, K} end

@inline shape(::Type{WMMAOp{M, N, K}}) where {M, N, K} = (M = M, N = N, K = K)

# convert_index_func: function used to transpose the index in case of a row-major layout
for (layout_type, wmma_layout_type, convert_index_func) in [
                                        (Layout.AlignedColMajor, WMMA.ColMajor, identity),
                                        (Layout.AlignedRowMajor, WMMA.RowMajor, x -> reverse(Tuple(x)))
                                       ]
    @eval begin
        @inline fragtype_a(::Type{WMMAOp{16, 16, 16}}, ::Type{$layout_type{Float16}}) = WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixA}
        @inline fragtype_b(::Type{WMMAOp{16, 16, 16}}, ::Type{$layout_type{Float16}}) = WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixB}
        @inline fragtype_accum(::Type{WMMAOp{16, 16, 16}}, ::Type{$layout_type{Float32}}) = WMMA.Fragment{16, 16, 16, 8, Float32, WMMA.Unspecified, WMMA.Accumulator}

        @inline function load_a(::Type{WMMAOp{M, N, K}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(Float16)
            return WMMA.load_a(ptr, size(workspace, 1), $wmma_layout_type, conf)
        end

        @inline function load_b(::Type{WMMAOp{M, N, K}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(Float16)
            return WMMA.load_b(ptr, size(workspace, 1), $wmma_layout_type, conf)
        end

        @inline function load_c(::Type{WMMAOp{M, N, K}}, ::Type{$layout_type{Float32}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(Float32)
            return WMMA.load_c(ptr, size(workspace, 1), $wmma_layout_type, conf)
        end

        @inline function store_d(::Type{WMMAOp{M, N, K}}, ::Type{$layout_type{Float32}}, workspace, frag, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}

            linear_base = linearise($convert_index_func(tile.base), size(workspace))
            linear_offset = linearise($convert_index_func(tile.offset), size(workspace))

            ptr = pointer(workspace, linear_base) + (linear_offset - 1) * sizeof(Float32)
            WMMA.store_d(ptr, frag, size(workspace, 1), $wmma_layout_type, conf)
        end
    end
end

function mma(::Type{WMMAOp{M, N, K}}, a_frag, b_frag, c_frag) where {M, N, K}
    conf = WMMA.Config{M, N, K, Float32}
    return WMMA.mma(a_frag, b_frag, c_frag, conf)
end

# -----------
# WMMAComplex
# -----------

struct WMMAComplexOp{M, N, K} end

@inline shape(::Type{WMMAComplexOp{M, N, K}}) where {M, N, K} = (M = M, N = N, K = K)

# convert_index_func: function used to transpose the index in case of a row-major layout
for (layout_type, wmma_layout_type, convert_index_func) in [
                                        (Layout.SplitColMajor, WMMA.ColMajor, identity),
                                        (Layout.SplitRowMajor, WMMA.RowMajor, x -> reverse(Tuple(x))),
                                       ]
    @eval begin
        @inline fragtype_a(::Type{WMMAComplexOp{16, 16, 16}}, ::Type{$layout_type{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixA}}
        @inline fragtype_b(::Type{WMMAComplexOp{16, 16, 16}}, ::Type{$layout_type{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, $wmma_layout_type, WMMA.MatrixB}}
        @inline fragtype_accum(::Type{WMMAComplexOp{16, 16, 16}}, ::Type{$layout_type{Float32}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 8, Float32, WMMA.Unspecified, WMMA.Accumulator}}

        @inline function load_a(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{16, 16, 16, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            return (WMMA.load_a(pointer(workspace, ind), size(workspace)[1], $wmma_layout_type, conf),
                    WMMA.load_a(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], $wmma_layout_type, conf))
        end

        @inline function load_b(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float16}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{16, 16, 16, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            return (WMMA.load_b(pointer(workspace, ind), size(workspace)[1], $wmma_layout_type, conf),
                    WMMA.load_b(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], $wmma_layout_type, conf))
        end

        @inline function load_c(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float32}}, workspace, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            return (WMMA.load_c(pointer(workspace, ind), size(workspace)[1], $wmma_layout_type, conf),
                    WMMA.load_c(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], $wmma_layout_type, conf))
        end

        @inline function store_d(::Type{WMMAComplexOp{M, N, K}}, ::Type{$layout_type{Float32}}, workspace, frag, tile::Tile) where {M, N, K}
            conf = WMMA.Config{M, N, K, Float32}
            ind = linearise($convert_index_func(tile.index), (size(workspace)[1], size(workspace)[2]))

            WMMA.store_d(pointer(workspace, ind), frag[1], size(workspace)[1], $wmma_layout_type, conf)
            WMMA.store_d(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), frag[2], size(workspace)[1], $wmma_layout_type, conf)
        end
    end
end

using LLVM

multiply_fp16(a::Float16, b::Float16) =
    Base.bitcast(Float16,
                 LLVM.Interop.@asmcall(
                                       "{mul.f16 \$0,\$1,\$2;}",
                                       "=h,h,h",
                                       false,
                                       Int16,
                                       Tuple{Int16, Int16},
                                       Base.bitcast(Int16, a),
                                       Base.bitcast(Int16, b)
                                      )
                )

@inline function mma(::Type{WMMAComplexOp{M, N, K}}, a_frag, b_frag, c_frag) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}

    c_re = c_frag[1]
    c_im = c_frag[2]

    c_re = WMMA.mma(a_frag[1],  b_frag[1], c_re, conf)
    c_re = WMMA.mma(multiply_fp16.(a_frag[2], Float16(-1)), b_frag[2], c_re, conf)

    c_im = WMMA.mma(a_frag[1], b_frag[2], c_im, conf)
    c_im = WMMA.mma(a_frag[2], b_frag[1], c_im, conf)

    return (c_re, c_im)
end

# --------
# WMMADual
# --------

struct WMMADualOp{M, N, K} end

@inline shape(::Type{WMMADualOp{M, N, K}}) where {M, N, K} = (M = M, N = N, K = K)

@inline fragtype_a(::Type{WMMADualOp{16, 16, 16}}, ::Type{Layout.SplitColMajor{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, WMMA.ColMajor, WMMA.MatrixA}}
@inline fragtype_b(::Type{WMMADualOp{16, 16, 16}}, ::Type{Layout.SplitColMajor{Float16}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 16, Float16, WMMA.ColMajor, WMMA.MatrixB}}
@inline fragtype_accum(::Type{WMMADualOp{16, 16, 16}}, ::Type{Layout.SplitColMajor{Float32}}) = NTuple{2, WMMA.Fragment{16, 16, 16, 8, Float32, WMMA.Unspecified, WMMA.Accumulator}}

@inline function load_a(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float16}}, workspace, tile::Tile) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    return (WMMA.load_a(pointer(workspace, ind), size(workspace)[1], WMMA.ColMajor, conf),
            WMMA.load_a(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], WMMA.ColMajor, conf))
end

@inline function load_b(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float16}}, workspace, tile::Tile) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    return (WMMA.load_b(pointer(workspace, ind), size(workspace)[1], WMMA.ColMajor, conf),
            WMMA.load_b(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], WMMA.ColMajor, conf))
end

@inline function load_c(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float32}}, workspace, tile::Tile) where {M, N, K}
    conf = WMMA.Config{M, N, K, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    return (WMMA.load_c(pointer(workspace, ind), size(workspace)[1], WMMA.ColMajor, conf),
            WMMA.load_c(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), size(workspace)[1], WMMA.ColMajor, conf))
end

@inline function store_d(::Type{WMMADualOp{M, N, K}}, ::Type{Layout.SplitColMajor{Float32}}, workspace, frag, tile::Tile) where {M, N, K}
    conf = WMMA.Config{M, N, K, Float32}
    ind = linearise(tile.index, (size(workspace)[1], size(workspace)[2]))

    WMMA.store_d(pointer(workspace, ind), frag[1], size(workspace)[1], WMMA.ColMajor, conf)
    WMMA.store_d(pointer(workspace, ind + size(workspace)[1] * size(workspace)[2]), frag[2], size(workspace)[1], WMMA.ColMajor, conf)
end

@inline function mma(::Type{WMMADualOp{M, N, K}}, a_frag, b_frag, c_frag) where {M, N, K}
    conf = WMMA.Config{16, 16, 16, Float32}

    c_re = c_frag[1]
    c_du = c_frag[2]

    c_re = WMMA.mma(a_frag[1],  b_frag[1], c_re, conf)

    c_du = WMMA.mma(a_frag[1], b_frag[2], c_du, conf)
    c_du = WMMA.mma(a_frag[2], b_frag[1], c_du, conf)

    return (c_re, c_du)
end

end


text:
#!/usr/bin/julia
#=
9x9 multiplication table in Julia
CC0, Wei-Lun Chao <bluebat@member.fsf.org>, 2018.
=#
# ./mt9x9.jl || julia mt9x9.jl

using Printf
for i in 1:3:9, j = 1:9, k = [i, i+1, i+2]
    @printf("%dx%d=%2d", k, j, k*j)
    print(k!=i+2 ? "\t" : j==9 ? "\n\n" : "\n")
end


text:
module TestCas
using ProjectiveGeometricAlgebra3d.CAS
using .CAS: BasisBlade, rev, dual, geomul, wedge, anti_wedge
using Accessors
using Test

function e(inds...;gendim, coeff=1)
    mask = CAS.create_mask(1:gendim, inds...)
    BasisBlade(coeff, mask)
end

@testset "BasisBlade" begin
    @testset "rev" begin
        e = BasisBlade(1, [false])
        e1 = BasisBlade(1, [true])
        e2 = BasisBlade(1, [false, true])
        e12 = BasisBlade(1, [true, true])
        e14 = BasisBlade(1, [true, false, false, true])
        e123 = BasisBlade(1, [true, false, false, true, true])
        e145 = BasisBlade(1, [true, false, false, true, true])
        e2456 = BasisBlade(1, [false, true, false, true, true, true])
        @test rev(e) == e
        @test rev(e1) == e1
        @test rev(e2) == e2
        @test rev(e12) == @set e12.coeff *= -1
        @test rev(e14) == @set e14.coeff *= -1
        @test rev(e123) == @set e123.coeff *= -1
        @test rev(e145) == @set e145.coeff *= -1
        @test rev(e2456) == e2456

    end

    b1 = BasisBlade(5, [true, true, ]) #e12
    b2 = BasisBlade(3, [true, true]) #scalar
    b = wedge(b1, b2)
    @test b.coeff == 0
    g = geomul(b1, b2, [1,1])
    @test g.coeff == -15
    @test g.mask == [false,false]
    g = geomul(b1, b2, [2,-1])
    @test g.coeff == 30
    @test g.mask == [false,false]

    b1 = BasisBlade(1, [true, false, true ]) # e13
    b2 = BasisBlade(1, [true, false, false]) # e1
    w = wedge(b1, b2)
    @test w.coeff == 0
    g = geomul(b1, b2, [1,1,1])
    @test g.coeff == -1
    @test g.mask == [false, false, true]
    g = geomul(b1, b2, [10,11,12])
    @test g.coeff == -10
    @test g.mask == [false, false, true]
    g = geomul(b1, b2, [0,11,12])
    @test g.coeff == 0

    b1 = BasisBlade(1, [true, false, true ]) #e13
    b2 = BasisBlade(1, [false, true, false]) #e2
    w = wedge(b1, b2)
    @test w.coeff == -1
    @test w.mask == [true,true,true]
    w = wedge(b2, b1)
    @test w.coeff == -1
    g = geomul(b1, b2, [1,1,1])
    @test g.coeff == -1
    @test g.mask == [true,true,true]

    b1 = BasisBlade(5, [true, true, ]) #e12
    b2 = BasisBlade(3, [false, false]) #scalar
    b = wedge(b1, b2)
    @test b.coeff == 15
    @test b.mask == [true,true]
end

@testset "dual 4d" begin
    @test dual(e(;gendim=4))    == e(1,2,3,4, gendim=4)

    @test dual(e(1;gendim=4))   == e(2,3,4, gendim=4)
    @test dual(e(2;gendim=4))   == e(1,3,4, gendim=4, coeff=-1)
    @test dual(e(3;gendim=4))   == e(1,2,4, gendim=4)
    @test dual(e(4;gendim=4))   == e(1,2,3, gendim=4, coeff=-1)

    @test dual(e(1,2;gendim=4)) == e(3,4, gendim=4)
    @test dual(e(1,3;gendim=4)) == e(2,4, gendim=4, coeff=-1)
    @test dual(e(1,4;gendim=4)) == e(2,3, gendim=4)

    @test dual(e(1,2,3,4, gendim=4))          == e(;gendim=4)

    @test dual(e(2,3,4, gendim=4))           == e(1;gendim=4)
    @test dual(e(1,3,4, gendim=4, coeff=-1))  == e(2;gendim=4)
    @test dual(e(1,2,4, gendim=4)          )  == e(3;gendim=4)
    @test dual(e(1,2,3, gendim=4, coeff=-1))  == e(4;gendim=4)

    @test dual(e(3,4, gendim=4)            )  == e(1,2;gendim=4)
    @test dual(e(2,4, gendim=4, coeff=-1)  )  == e(1,3;gendim=4)
    @test dual(e(2,3, gendim=4)            )  == e(1,4;gendim=4)
end

end#module


text:
for fconv in (:conv, :depthwiseconv), WT in (MaskedArray, MaskedSliceArray)
    local ∇conv_data, ∇conv_filter = Symbol.(:∇, fconv, [:_data, :_filter])
    conv_pullback, ∇conv_data_pullback = Symbol.([fconv, ∇conv_data], :_pullback)

    @eval begin
        NNlib.$fconv(x::AbstractArray{xT, N}, w::$WT{wT, N}; kwargs...) where {xT, wT, N} =
            NNlib.$fconv(x, freeze(w); kwargs...)
        NNlib.$fconv(x::AbstractArray{xT, N}, w::$WT{wT, N}, cdims::ConvDims;
                     kwargs...) where {xT, wT, N} = NNlib.$fconv(x, freeze(w), cdims; kwargs...)

        function ChainRulesCore.rrule(::typeof($fconv), x, w::$WT, cdims; kwargs...)
            y, pb = ChainRulesCore.rrule($fconv, x, freeze(w), cdims; kwargs...)
            function $conv_pullback(Δ)
                Δs = pb(Δ)
                return (Δs[1], Δs[2], w.mask .* Δs[3], Δs[4])
            end
            return y, $conv_pullback
        end

        function ChainRulesCore.rrule(::typeof($∇conv_data), x, w::$WT, cdims; kwargs...)
            y, pb = ChainRulesCore.rrule($∇conv_data, x, freeze(w), cdims; kwargs...)
            function $∇conv_data_pullback(Δ)
                Δs = pb(Δ)
                return (Δs[1], Δs[2], w.mask .* Δs[3], Δs[4])
            end
            return y, $∇conv_data_pullback
        end
    end
end


text:
#!/usr/bin/env julia

#=
--------------------------------------------------------------------------------------\n
Original Authors: BARC Project, Berkely MPC Laboratory -> https://github.com/MPC-Berkeley/barc
Modified by: Angshuman Goswami, Graduate Student, Clemson University
Date Create: 20/5/2016, Last Modified: 20/5/2016 \n
--------------------------------------------------------------------------------------\n
=#

using RobotOS
@rosimport barc.msg: ECU_raw, Ultrasound_xy, Z_DynBkMdl
@rosimport geometry_msgs.msg: Vector3
rostypegen()
using barc.msg
using geometry_msgs.msg
using JuMP
using Ipopt
#using AmplNLWriter
#using NLopt

# define model parameters
L_a     = 0.125         # distance from CoG to front axel
L_b     = 0.125         # distance from CoG to rear axel
dt      = 0.2              # time step of system

# preview horizon
N       =  12

# define targets [generic values]
#x_ref   = 10
#y_ref   = 10
#w       = sqrt(x_ref^2+y_ref^2)

# define decision variables
# states: position (x,y), yaw angle, and velocity
# inputs: acceleration, steering angle
println("Creating kinematic bicycle model ....")
mdl     = Model(solver = IpoptSolver(print_level=3))
#mdl     = Model(solver = BonminNLSolver())
#mdl     = Model(solver = NLoptSolver(algorithm=:LD_MMA))
@variable( mdl, x[1:(N+1)] )
@variable( mdl, y[1:(N+1)] )
@variable( mdl, psi[1:(N+1)] )
@variable( mdl, v[1:(N+1)] )
@variable( mdl, r[1:N] )
@variable( mdl, obsx1[1:N+1] )
@variable( mdl, obsy1[1:N+1] )
#@variable( mdl, obsx2[1:N+1] )
#@variable( mdl, obsy2[1:N+1] )
@variable( mdl, a[1:N] )
@variable( mdl, targetx[1:N+1] )
@variable( mdl, targety[1:N+1] )
@variable( mdl, epsi1[1:N+1] )
#@variable( mdl, epsi2[1:N+1] )
@variable( mdl, w1[1:N+1] )
@variable( mdl, w2[1:N+1] )

# define objective function
@NLobjective(mdl, Min, 1*(x[N+1] - targetx[1])^2 + 1*(y[N+1] - targety[1])^2 + 0.1*((r[N])^2 + (r[N-1])^2 + (r[N-2])^2 + (r[N-3])^2 + (r[N-4])^2 + (r[N-5])^2 + (r[N-6])^2 + (r[N-7])^2 + (r[N-8])^2 + (r[N-9])^2 + (r[N-10])^2 + (r[N-11])^2)  + 10^4*((epsi1[N+1]) + (epsi1[N]) + (epsi1[N-1]) + (epsi1[N-2]) + (epsi1[N-3]) + (epsi1[N-4]) + (epsi1[N-5]) + (epsi1[N-6]) + (epsi1[N-7]) + (epsi1[N-8]) + (epsi1[N-9]) + (epsi1[N-10]) + (epsi1[N-11])))

# define constraints
# define system dynamics
# Reference: R.Rajamani, Vehicle Dynamics and Control, set. Mechanical Engineering Series,
#               Spring, 2011, page 26
@NLparameter(mdl, x0                   == 0);              @NLconstraint(mdl, x[1]                  == x0);
@NLparameter(mdl, y0                   == 0);              @NLconstraint(mdl, y[1]                  == y0);
@NLparameter(mdl, psi0                 == 0);              @NLconstraint(mdl, psi[1]                == psi0 );
@NLparameter(mdl, v0                   == 0);              @NLconstraint(mdl, v[1]                  == v0);
@NLparameter(mdl, obsx10               == 20);             @NLconstraint(mdl, obsx1[1]              == obsx10);
@NLparameter(mdl, obsy10               == 0);              @NLconstraint(mdl, obsy1[1]              == obsy10);
#@NLparameter(mdl, obsx20               == 20);             @NLconstraint(mdl, obsx2[1]              == obsx20);
#@NLparameter(mdl, obsy20               == 0);              @NLconstraint(mdl, obsy2[1]              == obsy20);
@NLparameter(mdl, targetx0             == 10);             @NLconstraint(mdl, targetx[1]            == targetx0);
@NLparameter(mdl, targety0             == 10);             @NLconstraint(mdl, targety[1]            == targety0);
@NLparameter(mdl, w10                  == 10);             @NLconstraint(mdl, w1[1]                 == w10);
@NLparameter(mdl, w20                  == 10);             @NLconstraint(mdl, w2[1]                 == w20);

@NLexpression(mdl, obs1[i=1:N+1]       ,            (obsx1[i]-x[i])^2      + (obsy1[i]-y[i])^2)
#@defNLExpr(mdl, obs2[i=1:N+1]       ,            (obsx2[i]-x[i])^2      + (obsy2[i]-y[i])^2)

for i in 1:N
    @NLconstraint(mdl, x[i+1]                                 == x[i]      + dt*(v[i]*cos( psi[i] )  ) )
    @NLconstraint(mdl, y[i+1]                                 == y[i]      + dt*(v[i]*sin( psi[i] )  ) )
    @NLconstraint(mdl, psi[i+1]                               == psi[i]    + dt*(r[i])  )
    @NLconstraint(mdl, v[i+1]                                 == v[i]      + dt*(a[i])  )
    @NLconstraint(mdl, obsx1[i+1]                             == obsx1[i])
    @NLconstraint(mdl, obsy1[i+1]                             == obsy1[i])
    #@NLconstraint(mdl, obsx2[i+1]                             == obsx2[i])
    #@NLconstraint(mdl, obsy2[i+1]                             == obsy2[i])
    @NLconstraint(mdl, w1[i+1]                                == w1[i])
    @NLconstraint(mdl, w2[i+1]                                == w2[i])

    @NLconstraint(mdl, targetx[i+1]                           == targetx[i])
    @NLconstraint(mdl, targety[i+1]                           == targety[i])
    @NLconstraint(mdl, epsi1[i+1]                             == epsi1[i] + dt*(epsi1[i]))
    #@NLconstraint(mdl, epsi2[i+1]                             == epsi2[i])

    @NLconstraint(mdl,  0              <= a[i]                      <=1)
    @NLconstraint(mdl,  0              <= v[i]                      <=1.5)
    @NLconstraint(mdl, -2.7            <= r[i]                      <=2.7)

    @NLconstraint(mdl, epsi1[i] >=0)
    #@NLconstraint(mdl, epsi2[i] >=0)

    @NLconstraint(mdl, obs1[i] + epsi1[i] >=3)
    #@NLconstraint(mdl, obs2[i] + epsi2[i] >=4)

end

# status update
println("initial solve ...")
status = solve(mdl)
println("finished initial solve!")

function state_estimate_callback(msg::Z_DynBkMdl)
    # update mpc
    setvalue(x0,     msg.x)
    setvalue(y0,     msg.y)
    setvalue(psi0,   msg.psi)

end

function velocity_callback(msg::Vector3)
    # update mpc
    setvalue(v0,     msg.x)

end

function ultrasound1_xy_callback(msg::Ultrasound_xy)
    # update obstacles
    setvalue(obsx10             ,    msg.frontx)
    setvalue(obsy10             ,    msg.fronty)
    #setvalue(obsx20             ,    msg.frontleftx)
    #setvalue(obsy20             ,    msg.frontylefty)
    setvalue(w10                ,    msg.frontrightx)
    #setvalue(w20                ,    msg.frontrighty)

end

function target_callback(msg::Vector3)
    # update target
    setvalue(targetx0             ,    msg.x)
    setvalue(targety0             ,    msg.y)

end

function main()
    # initiate node, set up publisher / subscriber topics
    init_node("mpc1")
    pub = Publisher("ecu_1", ECU_raw, queue_size=10)
    pub1= Publisher("debug1_MPC", Z_DynBkMdl, queue_size=10)
    s1  = Subscriber("z_vhcl_1", Z_DynBkMdl, state_estimate_callback, queue_size=10)
    s2  = Subscriber("v_sim_1", Vector3, velocity_callback, queue_size=10)
    s3  = Subscriber("ultrasound_xy", Ultrasound_xy, ultrasound1_xy_callback, queue_size=10)
    s4  = Subscriber("target", Vector3, target_callback, queue_size=10)
    loop_rate = Rate(5)
    counter = 0
    while ! is_shutdown()
        # run mpc, publish command
	tic()
        status = solve(mdl)
        toc()

        # get optimal solutions
        a_opt   = getvalue(a[1])

        r_opt   = getvalue(r[1])

        if( status == :Infeasible)
           a_opt = 0
        end

        psi_opt = getvalue(psi[2])

        obsx1_opt   = getvalue(obsx1[1])
        obsy1_opt   = getvalue(obsy1[1])
        targetx_opt = getvalue(targetx[1])
        targety_opt = getvalue(targety[1])
        x_opt       = getvalue(epsi1[1])
        y_opt       = getvalue(epsi1[2])
        v_opt       = getvalue(epsi1[3])
        epsi1_opt   = getvalue(epsi1[4])
        # to ensure that the MPC runs multiple iteration before the car moves to get rid of intial infeasible results
        counter += 1
        if counter > 75
            # publish commands

              cmd = ECU_raw(a_opt, r_opt,psi_opt)
              publish(pub, cmd)

              #for debugging
              debug1 = Z_DynBkMdl(x_opt, y_opt, v_opt, obsx1_opt, obsy1_opt,epsi1_opt)
              publish(pub1, debug1)
        end
        rossleep(loop_rate)
    end
end

if ! isinteractive()
    main()
end

max_stars_repo_name: agoswam/Simulator

max_issues_repo_name: agoswam/Simulator

max_forks_repo_name: agoswam/Simulator


text:
#testbug_GB100
include("GB100.jl")
using DelimitedFiles
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV, Random, DataFrames, Gadfly
using MLBase

rng = Random.seed!(2019)

gbdt = GB100.GBDT(;
  sampling_rate         = 0.75,                # Sampling rate
  learning_rate         = 0.1,                 # Learning rate
  num_iterations        = 100,                # NuGBer of iterations
  tree_options          = Dict(
      :nsubfeatures     => 0,
      :min_samples_leaf => 1,
      :max_depth        => 5))
gbl                     = GB100.GBLearner(gbdt, :regression)

x, y                    = GB100.test_data(500; p = 400, sigma = sqrt(1))
x_test, y_test          = GB100.test_data(500; p = 400, sigma = sqrt(1))
mu                      = GB100.fried(x)
mu_test                 = GB100.fried(x_test)

# initialize the hypers
hypers = GB100.Hypers(x, y, 1, 1, 0.8)
hypers.σ = sqrt(1)
hypers.τ = 3 / 10 * var(y) / gbdt.num_iterations
hypers.α = 1
hypers.λ = 2.0  # hypers.λ = 1.0 or hypers.λ = nothing
hypers.s =  ones(size(x, 2)) / size(x, 2)
hypers

# initialize opt
opts = GB100.Opts()
opts.num_burn = 5

rmse2 = GB100.TBD_predict_and_rmse!(
    x,
    y,
    x_test,
    y_test,
    mu_test,
    gbdt,
    hypers,
    opts,
    rng)

hypers
sum(hypers.varcount) # 2765(burn 460)
sum(hypers.varcount[1:5]) # 1218(burn 460)
sum(hypers.s[1:5])

hypers.s
sum(hypers.s[1:5])

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
## set𝓖
function set𝓖!(aps::Vector{T},s::Symbol) where T<:AbstractElement
    𝓖 = QuadratureRule[s]
    nₑ = length(aps)
    nᵢ = nₑ*length(𝓖)
    data = Dict(:w=>zeros(nᵢ))
    data[:ξ] = zeros(nᵢ)
    isrk = false
    D = length(𝓖[1])
    if s∈(:SegRK2,:SegRK3,:SegRK4,:SegRK5,:TriRK3,:TriRK6,:TriRK13,:TetRK14,:TetRK27)
        isrk = true
        data[:wᵇ] = zeros(nᵢ)
        D > 3 ? data[:η] = zeros(nᵢ) : nothing
        D > 4 ? data[:γ] = zeros(nᵢ) : nothing
    else
        D > 2 ? data[:η] = zeros(nᵢ) : nothing
        D > 3 ? data[:γ] = zeros(nᵢ) : nothing
    end
    if T<:ReproducingKernel{SNode}
        index = zeros(Int,nᵢ)
        nₜ = 0
        n = 0
        𝝭 = Dict{Symbol,Vector{Float64}}()
        for ap in aps
            empty!(ap.𝓖)
            for ξ in 𝓖
                n += 1
                index[n] = nₜ
                nₜ += length(ap.𝓒)
                push!(ap.𝓖,SNode(n,data,index,𝝭))
                set𝓖(n,data,ξ,isrk)
            end
        end
    else
        n = 0
        for ap in aps
            empty!(ap.𝓖)
            for ξ in 𝓖
                n += 1
                push!(ap.𝓖,Node(n,data))
                set𝓖(n,data,ξ,isrk)
            end
        end
    end
end

function set𝓖(i::Int,data::Dict{Symbol,Vector{Float64}},ξ::NTuple{2,Float64},isrk::Bool=false)
    data[:w][i] = ξ[1]
    data[:ξ][i] = ξ[2]
end
function set𝓖(i::Int,data::Dict{Symbol,Vector{Float64}},ξ::NTuple{3,Float64},isrk::Bool=false)
    data[:w][i] = ξ[1]
    data[:ξ][i] = ξ[2]
    isrk ? data[:wᵇ][i] = ξ[3] : data[:η][i] = ξ[3]
end
function set𝓖(i::Int,data::Dict{Symbol,Vector{Float64}},ξ::NTuple{4,Float64},isrk::Bool=false)
    data[:w][i] = ξ[1]
    data[:ξ][i] = ξ[2]
    data[:η][i] = ξ[3]
    isrk ? data[:wᵇ][i] = ξ[4] : data[:γ][i] = ξ[4]
end
function set𝓖(i::Int,data::Dict{Symbol,Vector{Float64}},ξ::NTuple{5,Float64},isrk::Bool=false)
    data[:w][i] = ξ[1]
    data[:ξ][i] = ξ[2]
    data[:η][i] = ξ[3]
    data[:γ][i] = ξ[4]
    data[:wᵇ][i] = ξ[5]
end

## get𝓖
function get𝓖(a::T,b::S) where {T<:AbstractElement{:Seg2},S<:AbstractElement{:Poi1}}
    i = findfirst(x->x.id==b.𝓒[1].id, a.𝓒)
    if i ≠ nothing && i ≤ 2
        for ξ in b.𝓖
            i == 1 ? (ξ.ξ = -1.0;ξ.n₁ = -1.0) : (ξ.ξ = 1.0;ξ.n₁ = 1.0)
        end
        return b.𝓖
    else
        return nothing
    end
end

function get𝓖(a::T,b::S) where {T<:AbstractElement{:Tri3},S<:AbstractElement{:Poi1}}
    i = findfirst(x->x.id==b.𝓒[1].id, a.𝓒)
    if i ≠ nothing
        x₁ = a.𝓒[1].x
        y₁ = a.𝓒[1].y
        x₂ = a.𝓒[2].x
        y₂ = a.𝓒[2].y
        x₃ = a.𝓒[3].x
        y₃ = a.𝓒[3].y
        n₁₁ = y₃-y₂;n₂₁ = y₁-y₃;n₃₁ = y₂-y₁
        n₁₂ = x₂-x₃;n₂₂ = x₃-x₁;n₃₂ = x₁-x₂
        s₁₁ = -n₁₂;s₂₁ = -n₂₂;s₃₁ = -n₃₂
        s₁₂ =  n₁₁;s₂₂ =  n₂₁;s₃₂ =  n₃₁
        𝐿₁² = n₁₁^2+n₁₂^2
        𝐿₂² = n₂₁^2+n₂₂^2
        𝐿₃² = n₃₁^2+n₃₂^2
        if haskey(b.𝓖[1].data,:η)
            ξ₀ = b.𝓖[1]
            n = ξ₀.id
            data = ξ₀.data
            if isa(ξ₀,SNode)
                index = ξ₀.index
                𝝭 = ξ₀.𝝭
                push!(index,0)
                for (s,v) in data
                    push!(v,0.0)
                end
                ξ = SNode(n+1,data,index,𝝭)
            else
                ξ = Node(n+1,data)
            end
        else
            ξ = b.𝓖[1]
        end
        if i == 1
            ξ.ξ = 1.0
            ξ.η = 0.0
            ξ.Δn₁s₁ = n₂₁*s₂₁/𝐿₂² - n₃₁*s₃₁/𝐿₃²
            ξ.Δn₁s₂n₂s₁ = n₂₁*s₂₂/𝐿₂² + n₂₂*s₂₁/𝐿₂² - n₃₁*s₃₂/𝐿₃² - n₃₂*s₃₁/𝐿₃²
            ξ.Δn₂s₂ = n₂₂*s₂₂/𝐿₂² - n₃₂*s₃₂/𝐿₃²
        elseif i == 2
            ξ.ξ = 0.0
            ξ.η = 1.0
            ξ.Δn₁s₁ = n₃₁*s₃₁/𝐿₃² - n₁₁*s₁₁/𝐿₁²
            ξ.Δn₁s₂n₂s₁ = n₃₁*s₃₂/𝐿₃² + n₃₂*s₃₁/𝐿₃² - n₁₁*s₁₂/𝐿₁² - n₁₂*s₁₁/𝐿₁²
            ξ.Δn₂s₂ = n₃₂*s₃₂/𝐿₃² - n₁₂*s₁₂/𝐿₁²
        else
            ξ.ξ = 0.0
            ξ.η = 0.0
            ξ.Δn₁s₁ = n₁₁*s₁₁/𝐿₁² - n₂₁*s₂₁/𝐿₂²
            ξ.Δn₁s₂n₂s₁ = n₁₁*s₁₂/𝐿₁² + n₁₂*s₁₁/𝐿₁² - n₂₁*s₂₂/𝐿₂² - n₂₂*s₂₁/𝐿₂²
            ξ.Δn₂s₂ = n₁₂*s₁₂/𝐿₁² - n₂₂*s₂₂/𝐿₂²
        end
        return [ξ]
    else
        return nothing
    end
end

function get𝓖(a::T,b::S) where {T<:AbstractElement{:Tri3},S<:AbstractElement{:Seg2}}
    i = findfirst(x->x.id==b.𝓒[1].id, a.𝓒)
    j = findfirst(x->x.id==b.𝓒[2].id, a.𝓒)
    if i ≠ nothing && j ≠ nothing && i ≤ 3 && j ≤ 3
        𝐿 = get𝐿(b)
        x₁ = a.𝓒[1].x
        y₁ = a.𝓒[1].y
        x₂ = a.𝓒[2].x
        y₂ = a.𝓒[2].y
        x₃ = a.𝓒[3].x
        y₃ = a.𝓒[3].y
        for ξ in b.𝓖
            if i == 1
                ξ.ξ = (1.0-ξ.ξ)/2.0
                ξ.η = 1.0-ξ.ξ
                ξ.n₁ = (y₂-y₁)/𝐿
                ξ.n₂ = (x₁-x₂)/𝐿
                ξ.s₁ = (x₂-x₁)/𝐿
                ξ.s₂ = (y₂-y₁)/𝐿
            elseif i == 2
                ξ.η = (1.0-ξ.ξ)/2.0
                ξ.ξ = 0.0
                ξ.n₁ = (y₃-y₂)/𝐿
                ξ.n₂ = (x₂-x₃)/𝐿
                ξ.s₁ = (x₃-x₂)/𝐿
                ξ.s₂ = (y₃-y₂)/𝐿
            else
                ξ.ξ = (1.0+ξ.ξ)/2.0
                ξ.η = 0.0
                ξ.n₁ = (y₁-y₃)/𝐿
                ξ.n₂ = (x₃-x₁)/𝐿
                ξ.s₁ = (x₁-x₃)/𝐿
                ξ.s₂ = (y₁-y₃)/𝐿
            end
            ξ.w *= 0.5
            ξ.𝑤 = ξ.w*𝐿
        end
        return b.𝓖
    else
        return nothing
    end
end

## coordinate convertion
@inline getξ(a::T,b::T,ξ::N) where {T<:AbstractElement{:Tri3},N<:AbstractNode} = ξ
function getξ(a::T,b::S,ξ::N) where {T<:AbstractElement{:Tri3},S<:AbstractElement{:Seg2},N<:AbstractNode}
    i = findfirst(x->x.id==b.𝓒[1].id, a.𝓒)
    if i == 1
        return ((1.0-ξ.ξ)/2.0,(1.0+ξ.ξ)/2.0,0.0)
    elseif i == 2
        return ((1.0-ξ.ξ)/2.0,0.0,0.0)
    else
        return ((1.0+ξ.ξ)/2.0,0.0,0.0)
    end
end

## Quadrature Points
const QuadratureRule = Dict(
:PoiGI1 => ((1.0,1.0),),
:Seg100 => ((1.0,-1.0+0.02*i) for i in 0:100),
:SegGI1 => ((2.0,0.0),),
:SegGI2 =>
(
    (1.0,-0.5773502691896257645091487805),
    (1.0, 0.5773502691896257645091487805)
),
:SegGI3 =>
(
    (0.555555555555555555555555555556,-0.774596669241483377035853079957),
    (0.88888888888888888888888888889,0.0),
    (0.555555555555555555555555555556, 0.774596669241483377035853079957)
),
:SegGI4 =>
(
    (0.347854845137453857373063949222,-0.861136311594052575223946488893),
    (0.652145154862546142626936050778,-0.339981043584856264802665759103),
    (0.652145154862546142626936050778, 0.339981043584856264802665759103),
    (0.347854845137453857373063949222, 0.861136311594052575223946488893)
),
:SegGI5 =>
(
    (0.23692688505618908751426404072,-0.906179845938663992797626878299),
    (0.47862867049936646804129151484,-0.5384693101056830910363144207),
    (0.568888888888888888888888888889, 0.0),
    (0.47862867049936646804129151484, 0.5384693101056830910363144207),
    (0.23692688505618908751426404072, 0.906179845938663992797626878299)
),
:SegGI6 =>
(
    (0.171324492379170345040296142173,-0.932469514203152027812301554494),
    (0.360761573048138607569833513838,-0.6612093864662645136613995950),
    (0.46791393457269104738987034399,-0.238619186083196908630501721681),
    (0.46791393457269104738987034399, 0.238619186083196908630501721681),
    (0.360761573048138607569833513838, 0.66120938646626451366139959502),
    (0.171324492379170345040296142173, 0.932469514203152027812301554494)
),
:SegGI7 =>
(
    (0.129484966168869693270611432679,-0.949107912342758524526189684048),
    (0.27970539148927666790146777142,-0.741531185599394439863864773281),
    (0.38183005050511894495036977549,-0.405845151377397166906606412077),
    (0.417959183673469387755102040816, 0.0),
    (0.381830050505118944950369775489, 0.405845151377397166906606412077),
    (0.279705391489276667901467771424, 0.741531185599394439863864773281),
    (0.129484966168869693270611432679, 0.949107912342758524526189684048)
),
:SegGI8 =>
(
    (0.10122853629037625915253135431,-0.96028985649753623168356086857),
    (0.22238103445337447054435599443,-0.796666477413626739591553936476),
    (0.313706645877887287337962201987,-0.525532409916328985817739049189),
    (0.36268378337836198296515044928,-0.18343464249564980493947614236),
    (0.362683783378361982965150449277, 0.18343464249564980493947614236),
    (0.31370664587788728733796220199, 0.525532409916328985817739049189),
    (0.222381034453374470544355994426, 0.796666477413626739591553936476),
    (0.10122853629037625915253135431, 0.96028985649753623168356086857)
),
:SegGI9 =>
(
    (0.0812743883615744119718921581105,-0.968160239507626089835576202904),
    (0.180648160694857404058472031243,-0.83603110732663579429942978807),
    (0.260610696402935462318742869419,-0.613371432700590397308702039342),
    (0.31234707704000284006863040658,-0.32425342340380892903853801464),
    (0.330239355001259763164525069287, 0.0),
    (0.31234707704000284006863040658, 0.32425342340380892903853801464),
    (0.260610696402935462318742869419, 0.613371432700590397308702039342),
    (0.180648160694857404058472031243, 0.83603110732663579429942978807),
    (0.081274388361574411971892158111, 0.968160239507626089835576202904)
),
:SegGI10 =>
(
    (0.066671344308688137593568809893,-0.973906528517171720077964012085),
    (0.149451349150580593145776339658,-0.865063366688984510732096688424),
    (0.219086362515982043995534934228,-0.679409568299024406234327365115),
    (0.26926671930999635509122692157,-0.433395394129247190799265943166),
    (0.295524224714752870173892994651,-0.14887433898163121088482600113),
    (0.295524224714752870173892994651, 0.14887433898163121088482600113),
    (0.26926671930999635509122692157, 0.433395394129247190799265943166),
    (0.219086362515982043995534934228, 0.679409568299024406234327365115),
    (0.149451349150580593145776339658, 0.865063366688984510732096688424),
    (0.066671344308688137593568809893, 0.973906528517171720077964012085)
),
:TriGI3 =>
(
    (1/3,2/3,1/6),
    (1/3,1/6,2/3),
    (1/3,1/6,1/6)
),
:TriGI4 =>
(
    (-0.562500000000000,0.333333333333333,0.333333333333333),
    (0.520833333333333,0.600000000000000,0.200000000000000),
    (0.520833333333333,0.200000000000000,0.600000000000000),
    (0.520833333333333,0.200000000000000,0.200000000000000)
),
:TriGI6 =>
(
    (0.223381589678011,0.108103018168070,0.445948490915965),
    (0.223381589678011,0.445948490915965,0.108103018168070),
    (0.223381589678011,0.445948490915965,0.445948490915965),
    (0.109951743655322,0.816847572980459,0.091576213509771),
    (0.109951743655322,0.091576213509771,0.816847572980459),
    (0.109951743655322,0.091576213509771,0.091576213509771)
),
:TriGI7 =>
(
    (0.125939180544800,0.101286507323500,0.101286507323500),
    (0.125939180544800,0.797426985353100,0.101286507323500),
    (0.125939180544800,0.101286507323500,0.797426985353100),
    (0.132394152788500,0.470142064105100,0.059715871789800),
    (0.132394152788500,0.470142064105100,0.470142064105100),
    (0.132394152788500,0.059715871789800,0.470142064105100),
    (0.225000000000000,0.333333333333300,0.333333333333300)
),
:TriGI12 =>
(
    (0.116786275726379,0.501426509658179,0.249286745170910),
    (0.116786275726379,0.249286745170910,0.501426509658179),
    (0.116786275726379,0.249286745170910,0.249286745170910),
    (0.050844906370207,0.873821971016996,0.063089014491502),
    (0.050844906370207,0.063089014491502,0.873821971016996),
    (0.050844906370207,0.063089014491502,0.063089014491502),
    (0.082851075618374,0.053145049844817,0.310352451033784),
    (0.082851075618374,0.053145049844817,0.636502499121399),
    (0.082851075618374,0.310352451033784,0.053145049844817),
    (0.082851075618374,0.310352451033784,0.636502499121399),
    (0.082851075618374,0.636502499121399,0.053145049844817),
    (0.082851075618374,0.636502499121399,0.310352451033784)
),
:TriGI13 =>
(
    (0.053347235608800,0.065130102902200,0.065130102902200),
    (0.053347235608800,0.869739794195600,0.065130102902200),
    (0.053347235608800,0.065130102902200,0.869739794195600),
    (0.077113760890300,0.312865496004900,0.048690315425300),
    (0.077113760890300,0.638444188569800,0.312865496004900),
    (0.077113760890300,0.048690315425300,0.638444188569800),
    (0.077113760890300,0.638444188569800,0.048690315425300),
    (0.077113760890300,0.312865496004900,0.638444188569800),
    (0.077113760890300,0.048690315425300,0.312865496004900),
    (0.175615257433200,0.260345966079000,0.260345966079000),
    (0.175615257433200,0.479308067841900,0.260345966079000),
    (0.175615257433200,0.260345966079000,0.479308067841900),
    (-0.14957004446770,0.333333333333300,0.333333333333300)
),
:TriGI16 =>
(
    (0.144315607677787,0.333333333333333,0.333333333333333),
    (0.095091634267285,0.081414823414554,0.459292588292723),
    (0.095091634267285,0.459292588292723,0.081414823414554),
    (0.095091634267285,0.459292588292723,0.459292588292723),
    (0.103217370534718,0.658861384496480,0.170569307751760),
    (0.103217370534718,0.170569307751760,0.658861384496480),
    (0.103217370534718,0.170569307751760,0.170569307751760),
    (0.032458497623198,0.898905543365938,0.050547228317031),
    (0.032458497623198,0.050547228317031,0.898905543365938),
    (0.032458497623198,0.050547228317031,0.050547228317031),
    (0.027230314174435,0.008394777409958,0.263112829634638),
    (0.027230314174435,0.008394777409958,0.728492392955404),
    (0.027230314174435,0.263112829634638,0.008394777409958),
    (0.027230314174435,0.263112829634638,0.728492392955404),
    (0.027230314174435,0.728492392955404,0.008394777409958),
    (0.027230314174435,0.728492392955404,0.263112829634638)
),
:TriGI25 =>
(
    (0.090817990382754,0.333333333333333,0.333333333333333),
    (0.036725957756467,0.028844733232685,0.485577633383657),
    (0.036725957756467,0.485577633383657,0.028844733232685),
    (0.036725957756467,0.485577633383657,0.485577633383657),
    (0.045321059435528,0.781036849029926,0.109481575485037),
    (0.045321059435528,0.109481575485037,0.781036849029926),
    (0.045321059435528,0.109481575485037,0.109481575485037),
    (0.072757916845420,0.141707219414880,0.307939838764121),
    (0.072757916845420,0.141707219414880,0.550352941820999),
    (0.072757916845420,0.307939838764121,0.141707219414880),
    (0.072757916845420,0.307939838764121,0.550352941820999),
    (0.072757916845420,0.550352941820999,0.141707219414880),
    (0.072757916845420,0.550352941820999,0.307939838764121),
    (0.028327242531057,0.025003534762686,0.246672560639903),
    (0.028327242531057,0.025003534762686,0.728323904597411),
    (0.028327242531057,0.246672560639903,0.025003534762686),
    (0.028327242531057,0.246672560639903,0.728323904597411),
    (0.028327242531057,0.728323904597411,0.025003534762686),
    (0.028327242531057,0.728323904597411,0.246672560639903),
    (0.009421666963733,0.009540815400299,0.066803251012200),
    (0.009421666963733,0.009540815400299,0.923655933587500),
    (0.009421666963733,0.066803251012200,0.009540815400299),
    (0.009421666963733,0.066803251012200,0.923655933587500),
    (0.009421666963733,0.923655933587500,0.009540815400299),
    (0.009421666963733,0.923655933587500,0.066803251012200)
),
:QuadGI1 =>
(
    (2.0,0.0,0.0),
),
:QuadGI4 =>
(
    (1.0,-0.5773502691896258,-0.5773502691896258),
    (1.0, 0.5773502691896258,-0.5773502691896258),
    (1.0, 0.5773502691896258, 0.5773502691896258),
    (1.0,-0.5773502691896258, 0.5773502691896258)
),
:SegRK2 =>
(
    (1.0,-1.0,1.0),
    (1.0, 1.0,1.0)
),
:SegRK3 =>
(
    (1/3,-1.0,1.0),
    (4/3, 0.0,0.0),
    (1/3, 1.0,1.0)
),
:SegRK4 =>
(
    (1/6,-1.0,    1.0),
    (5/6,-0.2^0.5,0.0),
    (5/6, 0.2^0.5,0.0),
    (1/6, 1.0,    1.0)
),
:SegRK5 =>
(
    (1/10, -1.0,      1.0),
    (49/90,-(3/7)^0.5,0.0),
    (32/45, 0.0,      0.0),
    (49/90, (3/7)^0.5,0.0),
    (1/10,  1.0,      1.0)
),
:TriRK3 =>
(
    (1/3,1.0,0.0,1/2),
    (1/3,0.0,1.0,1/2),
    (1/3,0.0,0.0,1/2)
),
:TriRK6 =>
(
    (0.0,1.0,0.0,1/6),
    (0.0,0.0,1.0,1/6),
    (0.0,0.0,0.0,1/6),
    (1/3,0.0,0.5,2/3),
    (1/3,0.5,0.0,2/3),
    (1/3,0.5,0.5,2/3)
),
:TriRK13 =>
(
    (-0.0277777777777778,1.0000000000000000,0.0000000000000000, 1/20),
    (-0.0277777777777778,0.0000000000000000,1.0000000000000000, 1/20),
    (-0.0277777777777778,0.0000000000000000,0.0000000000000000, 1/20),
    ( 0.0296296296296297,0.0000000000000000,0.5000000000000000,16/45),
    ( 0.0296296296296297,0.5000000000000000,0.0000000000000000,16/45),
    ( 0.0296296296296297,0.5000000000000000,0.5000000000000000,16/45),
    ( 0.0907407407407407,0.0000000000000000,0.8273268353539885,49/180),
    ( 0.0907407407407407,0.0000000000000000,0.1726731646460116,49/180),
    ( 0.0907407407407407,0.1726731646460116,0.0000000000000000,49/180),
    ( 0.0907407407407407,0.8273268353539885,0.0000000000000000,49/180),
    ( 0.0907407407407407,0.8273268353539885,0.1726731646460116,49/180),
    ( 0.0907407407407407,0.1726731646460116,0.8273268353539885,49/180),
    ( 0.4500000000000000,0.3333333333333333,0.3333333333333333,  0.0)
)
)

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
#!/usr/bin/env julia

# Different implementations of gradient computation in Fourier space, with
# performance comparisons.

# Some sample benchmark results (on Julia 1.4):
#
#     Transforms: (RFFT, FFT, FFT)
#     Input type: Float64
#     Global dimensions: (64, 32, 64)  ->  (33, 32, 64)
#     MPI topology: 2D decomposition (2×1 processes)
#
#     gradient_global_view!...            184.853 μs (0 allocations: 0 bytes)
#     gradient_global_view_explicit!...   124.993 μs (0 allocations: 0 bytes)
#     gradient_local!...                  146.369 μs (0 allocations: 0 bytes)
#     gradient_local_parent!...           145.743 μs (0 allocations: 0 bytes)
#     gradient_local_linear!...           145.679 μs (0 allocations: 0 bytes)
#     gradient_local_linear_explicit!...  145.543 μs (0 allocations: 0 bytes)
#
# This was obtained when running julia with the default optimisation level -O2.
#
# If one uses -O1 instead, `gradient_global_view_explicit!` becomes much
# slower than the rest, indicating that there's a lot of compiler optimisations
# going on specifically in that function.

# See also https://jipolanco.github.io/PencilFFTs.jl/dev/examples/gradient/.

using BenchmarkTools
using MPI
using PencilFFTs

using AbstractFFTs: fftfreq, rfftfreq
using Printf: @printf
using Random: randn!

const PA = PencilFFTs.PencilArrays

const INPUT_DIMS = (64, 32, 64)

const DEV_NULL = @static Sys.iswindows() ? "nul" : "/dev/null"

function generate_wavenumbers_r2c(dims::Dims{3})
    box_size = (2π, 2π, 2π)  # Lx, Ly, Lz
    sample_rate = 2π .* dims ./ box_size

    # In our case (Lx = 2π and Nx even), this gives kx = [0, 1, 2, ..., Nx/2].
    kx = rfftfreq(dims[1], sample_rate[1])

    # In our case (Ly = 2π and Ny even), this gives
    # ky = [0, 1, 2, ..., Ny/2-1, -Ny/2, -Ny/2+1, ..., -1] (and similarly for kz).
    ky = fftfreq(dims[2], sample_rate[2])
    kz = fftfreq(dims[3], sample_rate[3])

    (kx, ky, kz)
end

# Compute and return ∇θ in Fourier space, using global views.
function gradient_global_view!(∇θ_hat::NTuple{3,PencilArray},
                               θ_hat::PencilArray, kvec_global)
    # Generate OffsetArrays that take global indices.
    θ_glob = global_view(θ_hat)
    ∇θ_glob = map(global_view, ∇θ_hat)

    @inbounds for (n, I) in enumerate(CartesianIndices(θ_glob))
        i, j, k = Tuple(I)  # global indices

        # Wave number vector associated to current Cartesian index.
        kx = kvec_global[1][i]
        ky = kvec_global[2][j]
        kz = kvec_global[3][k]

        u = im * θ_glob[n]

        ∇θ_glob[1][n] = kx * u
        ∇θ_glob[2][n] = ky * u
        ∇θ_glob[3][n] = kz * u
    end

    ∇θ_hat
end

# This is the fastest implementation, and I'm not sure why!
function gradient_global_view_explicit!(∇θ_hat::NTuple{3,PencilArray},
                                        θ_hat::PencilArray, kvec_global)
    # Generate OffsetArrays that take global indices.
    θ_glob = global_view(θ_hat)
    ∇θ_glob = map(global_view, ∇θ_hat)

    rng = axes(θ_glob)  # (i1:i2, j1:j2, k1:k2)

    # Note: since the dimensions in Fourier space are permuted as (z, y, x), it
    # is faster to loop with `k` as the fastest index.
    @assert permutation(θ_hat) === Permutation(3, 2, 1)

    @inbounds for i in rng[1], j in rng[2], k in rng[3]
        # Wave number vector associated to current Cartesian index.
        kx = kvec_global[1][i]
        ky = kvec_global[2][j]
        kz = kvec_global[3][k]

        u = im * θ_glob[i, j, k]

        ∇θ_glob[1][i, j, k] = kx * u
        ∇θ_glob[2][i, j, k] = ky * u
        ∇θ_glob[3][i, j, k] = kz * u
    end

    ∇θ_hat
end

# Compute and return ∇θ in Fourier space, using local indices.
function gradient_local!(∇θ_hat::NTuple{3,PencilArray}, θ_hat::PencilArray,
                         kvec_local)
    @inbounds for (n, I) in enumerate(CartesianIndices(θ_hat))
        i, j, k = Tuple(I)  # local indices

        # Wave number vector associated to current Cartesian index.
        kx = kvec_local[1][i]
        ky = kvec_local[2][j]
        kz = kvec_local[3][k]

        u = im * θ_hat[n]

        ∇θ_hat[1][n] = kx * u
        ∇θ_hat[2][n] = ky * u
        ∇θ_hat[3][n] = kz * u
    end

    ∇θ_hat
end

# Compute and return ∇θ in Fourier space, using local indices on the raw data
# (which takes permuted indices).
function gradient_local_parent!(∇θ_hat::NTuple{3,PencilArray},
                                θ_hat::PencilArray, kvec_local)
    θ_p = parent(θ_hat) :: Array
    ∇θ_p = parent.(∇θ_hat)

    perm = permutation(θ_hat)

    @inbounds for (n, I) in enumerate(CartesianIndices(θ_p))
        # Unpermute indices to (i, j, k)
        J = perm \ I

        # Wave number vector associated to current Cartesian index.
        i, j, k = Tuple(J)  # local indices
        kx = kvec_local[1][i]
        ky = kvec_local[2][j]
        kz = kvec_local[3][k]

        u = im * θ_p[n]

        ∇θ_p[1][n] = kx * u
        ∇θ_p[2][n] = ky * u
        ∇θ_p[3][n] = kz * u
    end

    ∇θ_hat
end

# Similar to gradient_local!, but avoiding CartesianIndices (slightly faster).
function gradient_local_linear!(∇θ_hat::NTuple{3,PencilArray},
                                θ_hat::PencilArray, kvec_local)
    # We want to iterate over the arrays in memory order to maximise
    # performance. For this we need to take into account the permutation of
    # indices in the Fourier-transformed arrays. By default, the memory order in
    # Fourier space is (z, y, x) instead of (x, y, z), but this is never assumed
    # below. The wave numbers must be permuted accordingly.
    perm = permutation(θ_hat)  # e.g. Permutation(3, 2, 1)
    kvec_perm = perm * kvec_local  # e.g. (kz, ky, kx)

    # Create wave number iterator.
    kvec_iter = Iterators.product(kvec_perm...)

    @inbounds for (n, kvec_n) in enumerate(kvec_iter)
        # Apply inverse permutation to the current wave number vector.
        # Note that this permutation has zero cost, since perm is a
        # compile-time constant!
        # (This can be verified by comparing the performance of this function
        # with the "explicit" variant of `gradient_local_linear`, below.)
        κ = perm \ kvec_n  # = (kx, ky, kz)

        u = im * θ_hat[n]

        # Note that this is very easy to generalise to N dimensions...
        ∇θ_hat[1][n] = κ[1] * u
        ∇θ_hat[2][n] = κ[2] * u
        ∇θ_hat[3][n] = κ[3] * u
    end

    ∇θ_hat
end

# Less generic version of the above, assuming that the permutation is (3, 2, 1).
# It's basically the same but probably easier to understand.
function gradient_local_linear_explicit!(∇θ_hat::NTuple{3,PencilArray},
                                         θ_hat::PencilArray, kvec_local)
    @assert permutation(θ_hat) === Permutation(3, 2, 1)

    # Create wave number iterator in (kz, ky, kx) order, i.e. in the same order
    # as the array data.
    kvec_iter = Iterators.product(kvec_local[3], kvec_local[2], kvec_local[1])

    @inbounds for (n, kvec_n) in enumerate(kvec_iter)
        kz, ky, kx = kvec_n
        u = im * θ_hat[n]
        ∇θ_hat[1][n] = kx * u
        ∇θ_hat[2][n] = ky * u
        ∇θ_hat[3][n] = kz * u
    end

    ∇θ_hat
end

MPI.Init()

# Input data dimensions (Nx × Ny × Nz)
dims = INPUT_DIMS

kvec = generate_wavenumbers_r2c(dims)  # as tuple of Frequencies
kvec_collected = collect.(kvec)        # as tuple of Vector

# Apply a 3D real-to-complex (r2c) FFT.
transform = Transforms.RFFT()

# MPI topology information
comm = MPI.COMM_WORLD
Nproc = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

# Disable output on all but one process.
rank == 0 || redirect_stdout(open(DEV_NULL, "w"))

# Let MPI_Dims_create choose the decomposition.
proc_dims = let pdims = zeros(Int, 2)
    MPI.Dims_create!(Nproc, pdims)
    pdims[1], pdims[2]
end

# Create plan
plan = PencilFFTPlan(dims, transform, proc_dims, comm)
println(plan, "\n")

# Allocate data and initialise field
θ = allocate_input(plan)
randn!(θ)

# Perform distributed FFT
θ_hat = plan * θ

# Compute and compare gradients using different methods.
# Note that these return a tuple of 3 PencilArrays representing a vector
# field.
∇θ_hat_base = allocate_output(plan, Val(3))
∇θ_hat_other = similar.(∇θ_hat_base)

# Local wave numbers: (kx[i1:i2], ky[j1:j2], kz[k1:k2]).
kvec_local = let rng = range_local(θ_hat)  # = (i1:i2, j1:j2, k1:k2)
    getindex.(kvec, rng)
end

gradient_global_view!(∇θ_hat_base, θ_hat, kvec)

@printf "%-40s" "gradient_global_view!..."
@btime gradient_global_view!($∇θ_hat_other, $θ_hat, $kvec_collected)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

@printf "%-40s" "gradient_global_view! (lazy)..."
@btime gradient_global_view!($∇θ_hat_other, $θ_hat, $kvec)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

# Slow version: kvec as a tuple of Vector
@printf "%-40s" "gradient_global_view_explicit!..."
@btime gradient_global_view_explicit!($∇θ_hat_other, $θ_hat, $kvec_collected)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

# Fast version: kvec as a tuple of lazy vectors (Frequencies)
@printf "%-40s" "gradient_global_view_explicit! (lazy)..."
@btime gradient_global_view_explicit!($∇θ_hat_other, $θ_hat, $kvec)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

@printf "%-40s" "gradient_local!..."
@btime gradient_local!($∇θ_hat_other, $θ_hat, $kvec_local)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

@printf "%-40s" "gradient_local_parent!..."
@btime gradient_local_parent!($∇θ_hat_other, $θ_hat, $kvec_local)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

@printf "%-40s" "gradient_local_linear!..."
@btime gradient_local_linear!($∇θ_hat_other, $θ_hat, $kvec_local)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

@printf "%-40s" "gradient_local_linear_explicit!..."
@btime gradient_local_linear_explicit!($∇θ_hat_other, $θ_hat, $kvec_local)
@assert all(∇θ_hat_base .≈ ∇θ_hat_other)

# Get gradient in physical space.
∇θ = plan \ ∇θ_hat_base


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_SimpleStateEvents

Simulate SimpleStateEvents model.

Demonstrates/tests simple state events
"""
module Simulate_SimpleStateEvents

import ..ModiaMath
include(joinpath("models", "SimpleStateEvents.jl"))
import .SimpleStateEvents

model  = SimpleStateEvents.Model()
result = ModiaMath.simulate!(model, stopTime=10.0; log=true)
ModiaMath.plot(result, (:s, :v, :f), heading="sSimulate_SimpleStateEvents.jl")

end


text:
export assign, assignlattice

function assign(n,sigma)
    magicnumber = 8
    K = Int64(ceil(sqrt(n*sigma/magicnumber)))
    println("We use $K x $K cells.")
    return assign(n,sigma,K,K)
end

function assignlattice(n,sigma)
    magicnumber = 8
    K = Int64(ceil(sqrt(n*sigma/magicnumber)))
    println("We use $K x $K cells.")
    return assignlattice(n,sigma,K,K)
end

function assignlattice(n,sigma,nx,ny)
    x0,x1 = 0.0,1.0
    y0,y1 = 0.0,1.0
    m = Int64(floor(sqrt(n)))
    edgesx = range(x0,stop=x1,length=m+1)
    edgesy = range(y0,stop=y1,length=m+1)
    x = zeros(m*m)
    y = zeros(m*m)
    r = ones(m*m)
    count = 1
    for i=1:m
        for j=1:m
            x[count] = edgesx[i]
            y[count] = edgesy[j]
            r[count] = 2*sigma/(m*m)
            count += 1
        end
    end
    delta1 = edgesx[2]-edgesx[1]
    radii = r[1]
    delta2 = delta1-2*radii
    println("Lattice distance: $delta1")
    println("Boundary distance $delta2")
    println("Radii: $radii")
    return assign(x0,x1,nx,y0,y1,ny,x,y,r)
end

function assign(n,sigma,nx,ny)
    x0,x1 = 0.0,1.0
    y0,y1 = 0.0,1.0
    x = x0.+(x1-x0)*rand(n)
    y = y0.+(y1-y0)*rand(n)
    r = ones(n)*(sigma/n)/2;
    return assign(x0,x1,nx,y0,y1,ny,x,y,r)
end

function assign(x0,x1,nx,y0,y1,ny,x,y,r)
    T = typeof(x[1])
    collectionx = [zeros(T,0) for i=1:nx,j=1:ny]
    collectiony = [zeros(T,0) for i=1:nx,j=1:ny]
    collectionr = [zeros(T,0) for i=1:nx,j=1:ny]
    lx = x1-x0
    ly = y1-y0
    edgesx = range(x0,stop=x1,length=nx+1)
    edgesy = range(y0,stop=y1,length=ny+1)
    @assert length(x) == length(y) == length(r)
    no = length(x)
    for n=1:no
        i = findindex(x[n],edgesx)
        j = findindex(y[n],edgesy)

        # Obstacle is in this cell but we have to check the neighbors
        for nexti in [-1,0,1]
            for nextj in [-1,0,1]
                shiftedx = x[n]+nexti*r[n]
                shiftedy = y[n]+nextj*r[n]
                shiftedxperiodic = x0+mod(shiftedx-x0,lx)
                shiftedyperiodic = y0+mod(shiftedy-y0,ly)
                shiftedi = mmod(i+nexti,nx)
                shiftedj = mmod(j+nextj,ny)

                if (edgesx[shiftedi]<=shiftedxperiodic<edgesx[shiftedi+1]) & (edgesy[shiftedj]<=shiftedyperiodic<edgesy[shiftedj+1])
                    deltax = ifelse(shiftedx<x0,lx,ifelse(shiftedx>x1,-lx,0))
                    deltay = ifelse(shiftedy<y0,ly,ifelse(shiftedy>y1,-ly,0))

                    append!(collectionx[shiftedi,shiftedj],x[n]+deltax)
                    append!(collectiony[shiftedi,shiftedj],y[n]+deltay)
                    append!(collectionr[shiftedi,shiftedj],r[n])
                end
            end
        end
    end

    Collection = Array{Obstacles,2}(undef,nx,ny)
    for i=1:nx
        for j=1:ny
            Collection[i,j] = Obstacles(a2sa(collectionx[i,j]),
                                        a2sa(collectiony[i,j]),
                                        a2sa(collectionr[i,j]),
                                        length(collectionx[i,j]))
        end
    end

    Meta = Metadata(nx,ny,edgesx,edgesy)

    return Collection, Meta
end


text:
function convert_arguments(P::PlotFunc, d::KernelDensity.UnivariateKDE)
    ptype = plottype(P, Lines) # choose the more concrete one
    to_plotspec(ptype, convert_arguments(ptype, d.x, d.density))
end

function convert_arguments(::Type{<:Poly}, d::KernelDensity.UnivariateKDE)
    points = Vector{Point2f0}(undef, length(d.x) + 2)
    points[1] = Point2f0(d.x[1], 0)
    points[2:end-1] .= Point2f0.(d.x, d.density)
    points[end] = Point2f0(d.x[end], 0)
    (points,)
end

function convert_arguments(P::PlotFunc, d::KernelDensity.BivariateKDE)
    ptype = plottype(P, Heatmap)
    to_plotspec(ptype, convert_arguments(ptype, d.x, d.y, d.density))
end

function searchrange(x, xlims)
    min, max = xlims
    i1 = searchsortedfirst(x, min)
    i2 = searchsortedlast(x, max)
    return i1:i2
end

function trim_density(k::KernelDensity.UnivariateKDE, xlims)
    range = searchrange(k.x, xlims)
    KernelDensity.UnivariateKDE(k.x[range], k.density[range])
end

function _density(x; trim = false)
    k = KernelDensity.kde(x)
    return trim ? trim_density(k, extrema_nan(x)) : k
end

"""
    density(values; npoints = 200, offset = 0.0, direction = :x)

Plot a kernel density estimate of `values`.
`npoints` controls the resolution of the estimate, the baseline can be
shifted with `offset` and the `direction` set to :x or :y.
`bandwidth` and `boundary` are determined automatically by default.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Density) do scene
    Theme(
        color = :gray85,
        strokecolor = :black,
        strokewidth = 1,
        strokearound = false,
        npoints = 200,
        offset = 0.0,
        direction = :x,
        boundary = automatic,
        bandwidth = automatic,
    )
end

function plot!(plot::Density{<:Tuple{<:AbstractVector}})
    x = plot[1]

    points = lift(x, plot.direction, plot.boundary, plot.offset,
        plot.npoints, plot.bandwidth) do x, dir, bound, offs, n, bw

        k = KernelDensity.kde(x;
            npoints = n,
            (bound === automatic ? NamedTuple() : (boundary = bound,))...,
            (bw === automatic ? NamedTuple() : (bandwidth = bw,))...,
        )

        ps = Vector{Point2f0}(undef, length(k.x) + 2)
        if dir === :x
            ps[1] = Point2f0(k.x[1], offs)
            ps[2:end-1] .= Point2f0.(k.x, k.density .+ offs)
            ps[end] = Point2f0(k.x[end], offs)
        elseif dir === :y
            ps[1] = Point2f0(offs, k.x[1])
            ps[2:end-1] .= Point2f0.(k.density .+ offs, k.x)
            ps[end] = Point2f0(offs, k.x[end])
        else
            error("Invalid direction $dir, only :x or :y allowed")
        end
        ps
    end

    linepoints = lift(points, plot.strokearound) do ps, sa
        if sa
            push!(copy(ps), ps[1])
        else
            ps[2:end-1]
        end
    end

    poly!(plot, points, color = plot.color, strokewidth = 0)
    lines!(plot, linepoints, color = plot.strokecolor, linewidth = plot.strokewidth)
end


text:
struct Expression{N,A<:AbstractOperator} <: AbstractExpression
	x::NTuple{N,Variable}
	L::A
	function Expression{N}(x::NTuple{N,Variable}, L::A) where {N,A<:AbstractOperator}

		# checks on L
		ndoms(L,1) > 1 && throw(ArgumentError(
        "cannot create expression with LinearOperator with `ndoms(L,1) > 1`"))

		#checks on x
		szL = size(L,2)
		szx = size.(x)
		check_sz = length(szx) == 1 ? szx[1] != szL : szx != szL
		check_sz && throw(ArgumentError(
        "Size of the operator domain $(size(L, 2)) must match size of the variable $(size.(x))"))

		dmL = domainType(L)
		dmx = eltype.(x)
		check_dm = length(dmx) == 1 ? dmx[1] != dmL : dmx != dmL
		check_dm && throw(ArgumentError(
	    "Type of the operator domain $(domainType(L)) must match type of the variable $(eltype.(x))"))

		new{N,A}(x,L)
	end
end

include("utils.jl")
include("multiplication.jl")
include("addition.jl")
include("abstractOperator_bind.jl")


text:
using Distributions
using LinearAlgebra: I # Identity matrix
using Random

using Minijyro

Random.seed!(42)

# Generate some data.
N = 100
D = 5
true_w = randn(D)
X = randn(N, D)
noise = 0.1 * randn(N)
y_obs = X * true_w + noise

@jyro function model(xs)
    D = size(xs)[2]
    w ~ MvNormal(zeros(D), I)
    y ~ MvNormal(xs * w, 0.1*I)
end

cond_model = condition(model, Dict(:y => y_obs))
samples, stats = nuts(cond_model, (X,), 1000)

@show abs.(true_w - mean(samples))


text:
using Test
using Altro
using BenchmarkTools
using TrajectoryOptimization

TEST_TIME = false

@testset "Benchmark Problems" begin
    include("benchmark_problems.jl")
    include(joinpath(@__DIR__,"..","examples","quickstart.jl"))
end

@testset "Solvers" begin
    include("constructors.jl")
    include("solve_tests.jl")
end

@testset "Solver Options" begin
    include("solver_opts.jl")
    include("solver_stats.jl")
end


text:
module Test11

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp
using Test

@testset "Test11" begin
    a = re"[a-z]"
    a.when = :le
    a = re.rep1(a)
    a.actions[:exit] = [:one]
    b = re"[a-z][a-z0-9]*"
    b.actions[:exit] = [:two]

    machine = Automa.compile(re.cat(a | b, '\n'))
    actions = Dict(
        :one => :(push!(logger, :one)),
        :two => :(push!(logger, :two)),
        :le  => :(p ≤ n))

    ctx = Automa.CodeGenContext(generator=:table)
    @test_throws ErrorException Automa.generate_exec_code(ctx, machine, actions)

    for generator in (:inline, :goto), checkbounds in (true, false), clean in (true, false)
        ctx = Automa.CodeGenContext(generator=generator, checkbounds=checkbounds, clean=clean)
        validate = @eval function (data, n)
            logger = Symbol[]
            $(Automa.generate_init_code(ctx, machine))
            p_end = p_eof = sizeof(data)
            $(Automa.generate_exec_code(ctx, machine, actions))
            return logger, cs == 0 ? :ok : cs < 0 ? :error : :incomplete
        end
        @test validate(b"a\n", 0) == ([:two], :ok)
        @test validate(b"a\n", 1) == ([:one, :two], :ok)
        @test validate(b"a1\n", 1) == ([:two], :ok)
        @test validate(b"aa\n", 1) == ([:two], :ok)
        @test validate(b"aa1\n", 1) == ([:two], :ok)
        @test validate(b"aa\n", 2) == ([:one, :two], :ok)
        @test validate(b"aa1\n", 2) == ([:two], :ok)
        @test validate(b"1\n", 1) == ([], :error)
    end
end

end


text:
module FreeRay

# using Julia package
using Parsers, DataFrames, Plots, SpecialFunctions

# Model utilities
include("utils.jl")

# Core modules
include("geometry.jl")
include("boundary.jl")
include("atmosphere.jl")
include("options.jl")
include("impcoeff.jl")
include("impedance.jl")
include("windprofile.jl")

include("environment.jl")
include("plotray.jl")
include("PlotShd.jl")
include("plottlr.jl")

end


text:
DEFAULT_VK_PHYSICAL_DEVICE_FEATURES = Dict(
    "robustBufferAccess"                         => 0,
    "fullDrawIndexUint32"                        => 0,
    "imageCubeArray"                             => 0,
    "independentBlend"                           => 0,
    "geometryShader"                             => 0,
    "tessellationShader"                         => 0,
    "sampleRateShading"                          => 0,
    "dualSrcBlend"                               => 0,
    "logicOp"                                    => 0,
    "multiDrawIndirect"                          => 0,
    "drawIndirectFirstInstance"                  => 0,
    "depthClamp"                                 => 0,
    "depthBiasClamp"                             => 0,
    "fillModeNonSolid"                           => 0,
    "depthBounds"                                => 0,
    "wideLines"                                  => 0,
    "largePoints"                                => 0,
    "alphaToOne"                                 => 0,
    "multiViewport"                              => 0,
    "samplerAnisotropy"                          => 0,
    "textureCompressionETC2"                     => 0,
    "textureCompressionASTC_LDR"                 => 0,
    "textureCompressionBC"                       => 0,
    "occlusionQueryPrecise"                      => 0,
    "pipelineStatisticsQuery"                    => 0,
    "vertexPipelineStoresAndAtomics"             => 0,
    "fragmentStoresAndAtomics"                   => 0,
    "shaderTessellationAndGeometryPointSize"     => 0,
    "shaderImageGatherExtended"                  => 0,
    "shaderStorageImageExtendedFormats"          => 0,
    "shaderStorageImageMultisample"              => 0,
    "shaderStorageImageReadWithoutFormat"        => 0,
    "shaderStorageImageWriteWithoutFormat"       => 0,
    "shaderUniformBufferArrayDynamicIndexing"    => 0,
    "shaderSampledImageArrayDynamicIndexing"     => 0,
    "shaderStorageBufferArrayDynamicIndexing"    => 0,
    "shaderStorageImageArrayDynamicIndexing"     => 0,
    "shaderClipDistance"                         => 0,
    "shaderCullDistance"                         => 0,
    "shaderFloat64"                              => 0,
    "shaderInt64"                                => 0,
    "shaderInt16"                                => 0,
    "shaderResourceResidency"                    => 0,
    "shaderResourceMinLod"                       => 0,
    "sparseBinding"                              => 0,
    "sparseResidencyBuffer"                      => 0,
    "sparseResidencyImage2D"                     => 0,
    "sparseResidencyImage3D"                     => 0,
    "sparseResidency2Samples"                    => 0,
    "sparseResidency4Samples"                    => 0,
    "sparseResidency8Samples"                    => 0,
    "sparseResidency16Samples"                   => 0,
    "sparseResidencyAliased"                     => 0,
    "variableMultisampleRate"                    => 0,
    "inheritedQueries"                           => 0,
    )


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualMachineScaleSetRollingUpgradesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Cancels the current virtual machine scale set rolling upgrade.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: OperationStatusResponse
"""
function virtualMachineScaleSetRollingUpgradesCancel(_api::VirtualMachineScaleSetRollingUpgradesApi, resourceGroupName::String, vmScaleSetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", OperationStatusResponse, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/rollingUpgrades/cancel", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the status of the latest virtual machine scale set rolling upgrade.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: RollingUpgradeStatusInfo
"""
function virtualMachineScaleSetRollingUpgradesGetLatest(_api::VirtualMachineScaleSetRollingUpgradesApi, resourceGroupName::String, vmScaleSetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", RollingUpgradeStatusInfo, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/rollingUpgrades/latest", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Starts a rolling upgrade to move all virtual machine scale set instances to the latest available Platform Image OS version. Instances which are already running the latest available OS version are not affected.
Param: resourceGroupName::String (required)
Param: vmScaleSetName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: OperationStatusResponse
"""
function virtualMachineScaleSetRollingUpgradesStartOSUpgrade(_api::VirtualMachineScaleSetRollingUpgradesApi, resourceGroupName::String, vmScaleSetName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", OperationStatusResponse, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/virtualMachineScaleSets/{vmScaleSetName}/osRollingUpgrade", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "vmScaleSetName", vmScaleSetName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualMachineScaleSetRollingUpgradesCancel, virtualMachineScaleSetRollingUpgradesGetLatest, virtualMachineScaleSetRollingUpgradesStartOSUpgrade


text:
@enum NodeType CALL CALLUNIVAR VARIABLE VALUE PARAMETER SUBEXPRESSION LOGIC COMPARISON EXTRA

export CALL, CALLUNIVAR, VARIABLE, VALUE, PARAMETER, SUBEXPRESSION, LOGIC, COMPARISON

immutable NodeData
    nodetype::NodeType
    index::Int
    parent::Int
end

export NodeData

# Only need to store parents, then transpose adjacency matrix!

# for CALL, index is into list of operators
# for CALLUNIVAR, index is into list of univariate operators
# for VARIABLE, index is variable index
# for VALUE, index is into list of constants
# for PARAMETER, index is into list of parameters
# for SUBEXPRESSION, index is into list of subexpressions
# for LOGIC, index is into list of logical operators (inputs and outputs are boolean)
# for COMPARISON, index is into lost of comparison operators
# for EXTRA, index is extension specific

const operators = [:+,:-,:*,:^,:/,:ifelse,:max,:min]
const USER_OPERATOR_ID_START = length(operators) + 1

const operator_to_id = Dict{Symbol,Int}()
for i in 1:length(operators)
    operator_to_id[operators[i]] = i
end
export operator_to_id, operators

const univariate_operators = Symbol[:+,:-,:abs]
const univariate_operator_to_id = Dict{Symbol,Int}(:+ => 1, :- => 2, :abs => 3)
const univariate_operator_deriv = Any[:(one(x)),:(-one(x)),:(ifelse(x >= 0, one(x),-one(x)))]
export univariate_operator_to_id, univariate_operators

for (op, deriv) in Calculus.symbolic_derivatives_1arg()
    push!(univariate_operators, op)
    push!(univariate_operator_deriv,deriv)
    univariate_operator_to_id[op] = length(univariate_operators)
end
const USER_UNIVAR_OPERATOR_ID_START = length(univariate_operators) + 1

const logic_operators = [:&&,:||]
const logic_operator_to_id = Dict{Symbol,Int}()
for i in 1:length(logic_operators)
    logic_operator_to_id[logic_operators[i]] = i
end
export logic_operator_to_id, logic_operators

const comparison_operators = [:<=,:(==),:>=,:<,:>]
const comparison_operator_to_id = Dict{Symbol,Int}()
for i in 1:length(comparison_operators)
    comparison_operator_to_id[comparison_operators[i]] = i
end
export comparison_operator_to_id, comparison_operators

# user-provided operators

const user_operator_map = Dict{Int,MathProgBase.AbstractNLPEvaluator}()

# we use the MathProgBase NLPEvaluator interface, where the
# operator takes the place of the objective function.
# users should implement eval_f and eval_grad_f for now.
# we will eventually support hessians too
function register_multivariate_operator(s::Symbol,f::MathProgBase.AbstractNLPEvaluator)
    !haskey(operator_to_id, s) || error("Operator $s has already been defined")
    id = length(operators)+1
    push!(operators,s)
    operator_to_id[s] = id
    user_operator_map[id] = f
end

export register_multivariate_operator

const user_univariate_operator_f = Dict{Int,Any}()
const user_univariate_operator_fprime = Dict{Int,Any}()
const user_univariate_operator_fprimeprime = Dict{Int,Any}()

# for univariate operators, just take in functions to evaluate
# zeroth, first, and second order derivatives
function register_univariate_operator(s::Symbol,f,fprime,fprimeprime)
    !haskey(univariate_operator_to_id, s) || error("Operator $s has already been defined")
    id = length(univariate_operators)+1
    push!(univariate_operators,s)
    univariate_operator_to_id[s] = id
    user_univariate_operator_f[id] = f
    user_univariate_operator_fprime[id] = fprime
    user_univariate_operator_fprimeprime[id] = fprimeprime
end

export register_univariate_operator

function has_user_multivariate_operators(nd::Vector{NodeData})
    for k in 1:length(nd)
        nod = nd[k]
        if nod.nodetype == CALL && nod.index >= USER_OPERATOR_ID_START
            return true
        end
    end
    return false
end


text:
# consts

function GAZCost( start::Int64, finish::Int64)
  #puzzle 1
  #return abs( start- finish )

  #puzzle2
  cost = 0
  for i=1:abs(start-finish)
    cost += i
  end
  return cost
end

# methods
function computeGAZ( entry::Vector{Int64}, pivot::Int64)

  totgaz = 0
  n = size(entry)[1]
  if ( pivot > n)
    return Inf;
  elseif ( pivot < 1 )
    return Inf;
  end

  for i=1:n
    totgaz += GAZCost( entry[i] , pivot )
  end
  return totgaz
end

function findGAZ( entry::Vector{Int64}, index::Int64)

  current = computeGAZ(entry, index)
  upward  = computeGAZ(entry, index+1)
  downward= computeGAZ(entry, index-1)

  if ( (current<upward) && (current<downward) )
    return current
  end

  next_index =  (upward > downward) ? (index - 1) : (index + 1)
  return findGAZ( entry, next_index)

end

function main()

  println("#######################")
  println("### AoC2021 - Day 7 ###")
  println("#######################")

  println(">> @arg1 : Input file")
  println()

  if ( size(ARGS)[1] != 1 )
    println("FATAL ERROR : Invalid number of args. Input file is required as @arg1.")
    quit()
  end

  fd = open(ARGS[1], "r")
  data = split(read(fd, String), ",")

  crabs = [parse(Int,x) for x in data]
  n = size(crabs)[1]
  if (n <= 1)
    println("INVALID INPUT : Not enough data for meaningful answer.")
    quit()
  end

  sort!(crabs)

  # find candidate median halfway between min & max
  min = crabs[1]
  max = crabs[n]
  median = ( min + max ) / 2
  println(" Min : $min , Max : $max , Median : $median")

  # total gaz
  minGaz = findGAZ( crabs, trunc(Int, median) )

  println("MIN GAZ : $minGaz")

  # close file
  close(fd)

end

@time main()


text:
module TestSearch3
include("../../../examples/search3.jl")
end  # module