text:
# This file is part of IntegerTriangles.jl.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module TrianglesTraitCard

using Nemo, TrianglesBase, TrianglesExamples, TrianglesUtils

export TraitCard, TransTraits, Traits, AssociatedTriangles, AllTraitCards

"""
Build the trait card of one or many triangles with Anums.
```
julia> TraitCard(LaguerreTriangle)
Triangle:  A021009 [1, 1, 1, 2, 4, 1, 6, 18, 9, 1]
Reverse:   A021010 [1, 1, 1, 1, 4, 2, 1, 9, 18, 6]
Inverse:   A021009 [1, -1, 1, 2, -4, 1, -6, 18, -9, 1]
RevInv:    A021010 [1, 1, -1, 1, -4, 2, 1, -9, 18, -6]
Diagonal:  A084950 [1, 1, 2, 1, 6, 4, 24, 18, 1, 120]
Sum:       A002720 [1, 2, 7, 34, 209, 1546, 13327, 130922, 1441729, 17572114]
EvenSum:   A331325 [1, 1, 3, 15, 97, 745, 6571, 65359, 723969, 8842257]
OddSum:    A331326 [0, 1, 4, 19, 112, 801, 6756, 65563, 717760, 8729857]
AltSum:    A009940 [1, 0, -1, -4, -15, -56, -185, -204, 6209, 112400]
DiagSum:   C001040 [1, 1, 3, 10, 43, 225, 1393, 9976, 81201, 740785]
Central:   A295383 [1, 4, 72, 2400, 117600, 7620480, 614718720, 59364264960]
LeftSide:  A000142 [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880]
RightSide: A000012 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
PosHalf:   A025167 [1, 3, 17, 139, 1473, 19091, 291793, 5129307, 101817089]
NegHalf:   A025166 [1, -1, 1, 7, -127, 1711, -23231, 334391, -5144063, 84149983]
TransUnos: A002720 [1, 2, 7, 34, 209, 1546, 13327, 130922, 1441729, 17572114]
TransSqrs: A105219 [0, 1, 8, 63, 544, 5225, 55656, 653023, 8379008, 116780049]
TransNat0: A103194 [0, 1, 6, 39, 292, 2505, 24306, 263431, 3154824, 41368977]
TransNat1: C000262 [1, 3, 13, 73, 501, 4051, 37633, 394353, 4596553, 58941091]
PolyVal2:  A087912 [1, 3, 14, 86, 648, 5752, 58576, 671568, 8546432, 119401856]
PolyVal3:  A277382 [1, 4, 23, 168, 1473, 14988, 173007, 2228544, 31636449]
```
"""
const ModuleTrianglesTraitCard = ""

const SEARCHONLINE = false # true

function TraitCard(T::ℚTri, name)
    @warn "WRONG TYPE! Will not process!"
end

"""
(SIGNATURES)
"""
function Traits(T::ℤTri, name, N=10, an=false)
    an = an && ! oeis_notinstalled()
    sep = " "

    println("\n=================")
    println(name)
    println()
    len = max(N, min(N - 1, length(T)))
    len = len < 10 ? len : 72

    S = Flat(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "Triangle", sep, SeqToFixlenString(S, len))

    R = Reverse(T); S = Flat(R)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "Reverse", sep, SeqToFixlenString(S, len))

    I = Inverse(T)
    if I != []

        S = Flat(I);
        anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
        println(anum, sep, name, sep, "Inverse", sep, SeqToFixlenString(S, len))

        RI = Reverse(I); S = Flat(RI);
        anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
        println(anum, sep, name, sep, "RevInv", sep, SeqToFixlenString(S, len))
    end

    IR = Inverse(R)
    if IR != []

        S = Flat(IR);
        anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
        println(anum, sep, name, sep, "InvRev", sep, SeqToFixlenString(S, len))
    end

    S = DiagonalTriangle(T); S = Flat(S)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "DiagTri", sep, SeqToFixlenString(S, len))

    S = ConvolutionTriangle(T); S = Flat(S)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "ConvTri", sep, SeqToFixlenString(S, len))

    S = PolyTriangle(T); S = Flat(S)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "PolyTri", sep, SeqToFixlenString(S, len))

    S = sum(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "Sum", sep, SeqToFixlenString(S, len))

    S = EvenSum(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "EvenSum", sep, SeqToFixlenString(S, len))

    S = OddSum(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "OddSum", sep, SeqToFixlenString(S, len))

    S = AltSum(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "AltSum", sep, SeqToFixlenString(S, len))

    S = DiagSum(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "DiagSum", sep, SeqToFixlenString(S, len))

    S = Middle(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "Middle", sep, SeqToFixlenString(S, len))

    S = Central(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "Central", sep, SeqToFixlenString(S, len))

    S = LeftSide(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "LeftSide", sep, SeqToFixlenString(S, len))

    S = RightSide(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "RightSide", sep, SeqToFixlenString(S, len))

    # --------
    P = Polynomial(T)
    # --------

    S = PosHalf(P)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "PosHalf", sep, SeqToFixlenString(S, len))

    S = NegHalf(P)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "NegHalf", sep, SeqToFixlenString(S, len))

    S = BinConv(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "BinConv", sep, SeqToFixlenString(S, len))

    S = IBinConv(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "IBinConv", sep, SeqToFixlenString(S, len))

    S = Trans(T, [ZZ(n^2) for n = 0:19])
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "TransSqrs", sep, SeqToFixlenString(S, len))

    S = Trans(T, [ZZ(n) for n = 0:19])
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "TransNat0", sep, SeqToFixlenString(S, len))

    S = Trans(T, [ZZ(n) for n = 1:19])
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "TransNat1", sep, SeqToFixlenString(S, len))

    S = PolyVal2(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "PolyVal2", sep, SeqToFixlenString(S, len))

    S = PolyVal3(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "PolyVal3", sep, SeqToFixlenString(S, len))

    S = PolyValn(T)
    anum = an ? GetSeqnum(S, SEARCHONLINE) : ""
    println(anum, sep, name, sep, "PolyValn", sep, SeqToFixlenString(S, len))
end

"""
Print the standard traits generated by linear transformation.
```
julia> TransTraits(BinomialTriangle)
BinomialTriangle
TransSqrs: A001788 [0, 1, 6, 24, 80, 240, 672, 1792, 4608, 11520]
TransNat0: A001787 [0, 1, 4, 12, 32, 80, 192, 448, 1024, 2304]
TransNat1: A001792 [1, 3, 8, 20, 48, 112, 256, 576, 1280, 2816]
```
"""
function TransTraits(triangle::Function)
    dim = 19
    T = triangle(dim)
    println("$triangle ")
    len = min(10, min(dim - 1, length(T)))

    S = Trans(T, [ZZ(1) for n = 0:dim]); anum = GetSeqnum(S, false)
    print("TransUnos: ", anum, " "); S[1:len] |> Println
    S = Trans(T, [ZZ((-1)^n) for n = 0:dim]); anum = GetSeqnum(S, false)
    print("TransAlts: ", anum, " "); S[1:len] |> Println
    S = Trans(T, [ZZ(n^2) for n = 0:dim]);    anum = GetSeqnum(S, false)
    print("TransSqrs: ", anum, " "); S[1:len] |> Println
    S = Trans(T, [ZZ(n) for n = 0:dim]);      anum = GetSeqnum(S, false)
    print("TransNat0: ", anum, " "); S[1:len] |> Println
    S = Trans(T, [ZZ(n) for n = 1:dim]);      anum = GetSeqnum(S, false)
    print("TransNat1: ", anum, " "); S[1:len] |> Println
    println()
end

"""
(SIGNATURES)
"""
function AssociatedTriangles(T::ℤTri, N)

    len = max(N, min(N - 1, length(T)))

    println("\nTriangle:\n")
    for row in T[1:len] Println(row) end
    println()
    print("Flat triangle:  ");
    S = Flat(T);
    S[1:len] |> Println

    println("\nInverse triangle:\n")
    S = InverseTriangle(T)
    if S != []
        for row in S[1:len] Println(row) end
        println()
        S = Flat(S); print("Inverse:  "); S[1:len] |> Println
    end

    println("\nDiagonal triangle:\n")
    S = DiagonalTriangle(T)
    Println.(S[1:len])
    println()
    S = Flat(S); print("Diagonal:  "); S[1:len] |> Println

    println("\nPolynomial values, array:\n")
    Println.(PolyArray(T))

    println("\nPolynomial values, triangle:\n")
    S = PolyTriangle(T)
    Println.(S[1:len])
    println()
    S = Flat(S); print("PolyVal:  "); S[1:len] |> Println
end

"""
Prints a list of traits of the triangle with Anums.
"""
function TraitCard(triangle::Function)
    dim = 19
    T = triangle(dim)
    name = "$triangle"
    c = findfirst("Triangle", name)
    Traits(T, name[1:c[1]-1], dim, true)
end

"""
Prints a list of traits of the triangle without Anums
and, if assoc=true, also a list of associated triangles.
"""
function TraitCard(triangle::Function, dim, assoc=false)
    T = triangle(dim)
    name = "$triangle"
    c = findfirst("Triangle", name)
    Traits(T, name[1:c[1]-1], dim)
    AssociatedTriangles(T, dim)
end

function AllTraitCards()
    for T in TRIANGLES
        TraitCard(T)
    end
end

# START-TEST-########################################################
using TrianglesExamples

function test()
    dim = 12
    TraitCard(BinomialTriangle, dim)
    TraitCard(LaguerreTriangle, dim)
    TraitCard(LahTriangle, dim)
    TraitCard(CatalanTriangle, dim)
    TraitCard(MotzkinTriangle, dim)
    TraitCard(NarayanaTriangle, dim)
    TraitCard(SchroederBigTriangle, dim)
    TraitCard(JacobsthalTriangle, dim)
    TraitCard(FibonacciTriangle, dim)
    TraitCard(EulerianTriangle, dim)
    TraitCard(UniTriangle, dim)
end

function demo()
    # TransTraits(BinomialTriangle)
    # Function composition
    #EulerianClassic = Reverse ∘ EulerianTriangle
    #TraitCard(EulerianClassic)
end

function perf()
    #AllTraitCards()
end

function main()
    #test()
    demo()
    perf()
end

main()

end # module


text:
using PlutoSliderServer
ENV["JULIA_DEBUG"] = PlutoSliderServer

test_dir = tempname(cleanup=false)
cp(@__DIR__, test_dir)

try
    # open the folder on macos:
    run(`open $(test_dir)`)
catch end

port = 2345

# cdn = nothing
cdn = "http://127.0.0.1:$(port)/pluto_asset/"

PlutoSliderServer.run_directory(test_dir;
    static_export=true, run_server=true,
    SliderServer_serve_static_export_folder=true, SliderServer_port=port, SliderServer_host="127.0.0.1", Export_baked_state=false, Export_slider_server_url="http://127.0.0.1:$(port)/", Export_pluto_cdn_root=cdn)

ext: jl


text:
#TODO: accessing parents of attributes
#TODO: implement remaining xpath functions
#TODO: parenthesized expressions
#TODO: xmlns namespace parsing
#TODO: correct ordering of output
#TODO: $QName string interpolation
#TODO: " and '

# XPath Spec: http://www.w3.org/TR/xpath/

import Base: |

const xpath_axes = Dict(
    "ancestor" => :ancestor,
    "ancestor-or-self" => :ancestor_or_self,
    "attribute" => :attribute,
    "child" => :child,
    "descendant" => :descendant,
    "descendant-or-self" => :descendant_or_self,
    "following" => :following,
    "following-sibling" => :following_sibling,
#    "namespace" => :namespace,
    "parent" => :parent,
    "preceding" => :preceding,
    "preceding-sibling" => :preceding_sibling,
    "self" => :self)

const xpath_types = Dict(
    "comment" => (:comment,AbstractString),
    "text" => (:text,AbstractString),
#    "processing-instruction" => (:processing_instruction, ??),
    "node" => (:node,Any))

const xpath_functions = Dict( # (name, min args, max args)
    #node-set
    "last" => (:last,0,0,Int),
    "position" => (:position,0,0,Int),
    "count" => (:count,1,1,Int),
    "local-name" => (:local_name,0,1,AbstractString),
    #"namespace-uri" => (:namespace_uri,0,1,AbstractString),
    "name" => (:name,0,1,AbstractString),

    #string
    "string" => (:string_fn,0,1,AbstractString),
    "concat" => (:concat,2,typemax(Int),AbstractString),
    "starts-with" => (:startswith,2,2,Bool),
    "contains" => (:contains,2,2,Bool),
    "substring-before" => (:substring_before,2,2,AbstractString),
    "substring-after" => (:substring_after,2,2,AbstractString),
    "substring" => (:substring,2,3,AbstractString),
    "string-length" => (:string_length,0,1,Int),
    "normalize-space" => (:normalize_space,0,1,AbstractString),
    "translate" => (:translate,3,3,AbstractString),

    #boolean
    "boolean" => (:bool,1,1,Bool),
    "not" => (:not,1,1,Bool),
    "true" => (:true_,0,0,Bool),
    "false" => (:false_,0,0,Bool),
    #"lang" => (:lang,1,1,Bool),

    #number
    "number" => (:number_fn,0,1,Number),
    "sum" => (:sum,1,1,Float64),
    "floor" => (:floor,1,1,Int),
    "ceiling" => (:ceiling,1,1,Int),
    "round" => (:round,1,1,Float64))

macro xpath_str(xpath)
    xp, returntype = xpath_parse(xpath, true)
    :( XPath{$(typeof(xpath)), $(returntype)}($(xp)) )
end

function consume_whitespace(xpath, k)
    #consume leading space
    res = iterate(xpath, k)
    while res !== nothing
        c, k2 = res
        if !isspace(c)
            break
        end
        k = k2
        res = iterate(xpath, k)
    end
    k
end

const xpath_separators = Set(['+','(',')','[',']','<','>','!','=','|','/','*',','])

function xpath_parse(xpath::T, ismacro=false) where T<:AbstractString
    k = firstindex(xpath)
    k, parsed, returntype, has_last_fn = xpath_parse_expr(xpath, k, 0, ismacro)
    if k ≤ ncodeunits(xpath) # assume that xpath_parse_expr uses proper character indexing
        error("failed to parse to the end of the xpath (stopped at $k)")
    end
    return parsed, returntype
end

#function xpath_parse_filters{T<:AbstractString}(xpath::T)
#    k = consume_whitespace(xpath, start(xpath))
#    c,k = next(xpath,k)
#    if c != '['
#        error("expected first character of xpath filter to be a '['")
#    end
#    k, parsed, returntype, has_last_fn = xpath_parse_expr(xpath, k, 0, false)
#    if done(xpath,k)
#        error("unexpected end to xpath filter, expected a ']' (stopped at $k)")
#    end
#    c,k = next(xpath,k)
#    if !done(xpath,consume_whitespace(xpath,k)) || c != ']'
#        error("failed to parse to the end of the xpath filter (stopped at $k)")
#    end
#    return parsed, returntype
#end

macro xpath_parse(arg1, arg2)
    :(
        if $(esc(:ismacro))
            a2 = $(esc(arg2))
            if !isa(a2,Expr) && !isa(a2,AbstractString)
                a2 = Expr(:quote,a2)
            end
            $(esc(:parsed)) = Expr(:call, :push!, $(esc(:parsed)), Expr(:tuple,Expr(:quote,$(esc(arg1))),a2))
        else
            push!($(esc(:parsed))::Vector{SymbolAny}, ($(esc(arg1)),$(esc(arg2))))
        end
    )
end
macro xpath_fn(arg1, arg2)
    :(
        if $(esc(:ismacro))
            a2 = $(esc(arg2))
            if !isa(a2,Expr) && !isa(a2,AbstractString)
                a2 = Expr(:quote,a2)
            end
            Expr(:tuple,Expr(:quote,$(esc(arg1))),a2)
        else
            ($(esc(arg1)),$(esc(arg2)))
        end
    )
end

function xpath_parse(xpath::T, k, ismacro) where T<:AbstractString
    if ismacro
        parsed = :(SymbolAny[])
    else
        parsed = SymbolAny[]
    end
    k = consume_whitespace(xpath, k)
    res = iterate(xpath, k)
    if res === nothing
        error("empty xpath expressions is not valid")
    end
    # 1. Consume root node
    c, k2 = res
    if c == '/'
        @xpath_parse :root :node
        k = k2
    end
    returntype::DataType = ETree
    first::Bool = true
    while iterate(xpath, k) !== nothing
        # i..j has text, k is current character
        havename::Bool = false
        axis::Symbol = :child
        colonpos::Int = 0
        doublecolon::Bool = false
        dot::Bool = false
        parens::Bool = false
        name::T = ""
        c, k2 = iterate(xpath, k)
        i = k
        j = 0
        if c == '/'
            @xpath_parse :descendant_or_self :node
            returntype = Any
            i = k = k2 #advance to next
        end
        # 2. Consume node name
        while iterate(xpath, k) !== nothing
            c, k2 = iterate(xpath,k)
            if c == ':'
                # 2a. Consume axis name
                if !havename && j == 0
                    error("unexpected : at $k $i:$j")
                end
                if colonpos != 0
                    if !havename
                        name = xpath[i:j]
                    end
                    if doublecolon
                        error("unexpected :: at $k")
                    end
                    havename = false
                    axis_ = get(xpath_axes, name, nothing)
                    if axis_ === nothing
                        error("unknown axis $name")
                    end
                    axis = axis_::Symbol
                    colonpos = 0
                    doublecolon = true
                    i = k2
                    j = 0
                else # colon == 0
                    colonpos = k
                end #if
            else # c != ":"
                if colonpos != 0
                    j = colonpos
                    colonpos = 0
                end #if
                # 2b. Consume node name
                if j == 0 && c == '*'
                    havename = true
                    name = "*"
                    i = k = k2
                    break
                elseif isspace(c) || in(c, xpath_separators)
                    if j != 0
                        @assert(!havename)
                        havename = true
                        name = xpath[i:j]
                        j = 0
                    end
                    if c == '('
                        k2 = consume_whitespace(xpath, k2)
                        if iterate(xpath,k2) === nothing
                            error("unexpected end to xpath after (")
                        end
                        c, k3 = iterate(xpath,k2)
                        if c != ')'
                            error("unexpected character before ) in nodetype() expression at $k2")
                        end
                        k = k3
                        parens = true
                        break
                    elseif !isspace(c)
                        break
                    end #if
                    i = k2
                elseif havename # && !isspace && !separator
                    break
                elseif c == '-' && j == 0
                    error("TODO: -negation")
                else # text character
                    j = k
                end #if
            end #if
            k = k2
        end # if
        if !havename
            if j!=0
                havename = true
                name = xpath[i:j]
            elseif first != true
                if iterate(xpath,k) === nothing
                    error("xpath should not end with a /")
                end
                error("expected name before $c at $k")
            else
                break
            end
        elseif j!=0
            @assert(false)
        end
        first = false
        if parens
            nodetype = get(xpath_types, name, nothing)
            if nodetype === nothing
                error("unknown node type or function $name at $k")
            end
            @xpath_parse axis nodetype[1]::Symbol
            returntype = nodetype[2]::DataType
        elseif name[1] == '.'
            if doublecolon
                error("xml names may not begin with a . (at $k)")
            elseif length(name) == 2 && name[2] == '.'
                @xpath_parse :parent :element
                returntype = ETree
            elseif length(name) == 1
                @xpath_parse :self :node
            else
                error("xml names may not begin with a . (at $k)")
            end
        elseif name[1] == '@' || axis == :attribute
            if axis != :attribute
                k2 = consume_whitespace(name, 2)
                name = name[k2:end]
            end
            if name == "*"
                @xpath_parse :attribute nothing
            else
                @xpath_parse :attribute name
            end
            returntype = AbstractString
        elseif name[1] == '$'
            @xpath_parse axis :element
            @xpath_parse :name Expr(:call, :string, esc(Symbol(name[2:end])))
        else
            @xpath_parse axis :element
            if name != "*"
                @xpath_parse :name name
            end
            returntype = ETree
        end #if
        while iterate(xpath,k) !== nothing
            c, k2 = iterate(xpath,k)
            if isspace(c)
                k = k2
                continue
            elseif c == '/'
                if iterate(xpath,k) === nothing
                    error("xpath should not end with a /")
                #elseif returntype !== ETree # this is a valid XPath
                #    error("xpath has an unexpected / at $k -- previous selector does not return a Node")
                end
                k = k2
                break
            elseif c == '|'
                k, parsed2, rt2 = xpath_parse(xpath, k2, ismacro)
                if rt2 !== returntype
                    returntype = Any
                    #error("xpath return types on either side of | don't appear to match")
                end
                if ismacro
                    parsed = :( push!(SymbolAny[], (:(|), ($(parsed), $(parsed2)))) )
                else
                    parsed = push!(SymbolAny[], (:(|), (parsed, parsed2)))
                end
                return k, parsed, returntype
            elseif c == '['
                i = k
                k = k2
                k, filter, rettype, has_last_fn = xpath_parse_expr(xpath, k, 0, ismacro)
                if has_last_fn
                    @xpath_parse :filter_with_last filter
                else
                    @xpath_parse :filter filter
                end
                k = consume_whitespace(xpath, k)
                if iterate(xpath, k) === nothing
                    error("unmatched ] at $i")
                end
                c, k2 = iterate(xpath, k)
                if (c != ']')
                    error("expected matching ] at $k for [ at $i, found $c")
                end
                k = k2
                if iterate(xpath, k) !== nothing
                    c, k2 = iterate(xpath, k)
                end
            else
                return k, parsed, returntype #hope something else can parse it
            end #if
        end #while
    end # while
    return k, parsed, returntype
end # function

function xpath_parse_expr(xpath::T, k, precedence::Int, ismacro) where T<:AbstractString
    i = k = consume_whitespace(xpath, k)
    j = 0
    prevtokenspecial = true
    while iterate(xpath, k) !== nothing
        c, k2 = iterate(xpath, k)
        if prevtokenspecial && c == '*'
            nothing
        elseif c == '@' || c == ':'
            prevtokenspecial = true
            k = k2
            continue
        elseif c == '"' || c == '\''
            c2::Char = 0
            escape = false
            while c2 != c && !escape
                if ismacro && c2 == '\\' && !escape
                    escape = true
                else
                    escape = false
                end
                j = k
                k = k2
                if iterate(xpath, k) === nothing
                    error("unterminated string literal $c at $k")
                end
                c2, k2 = iterate(xpath, k)
            end
            k = k2
            break
        elseif isspace(c) || in(c, xpath_separators)
            if c == '/'
                j = k
            end
            break
        end
        prevtokenspecial = false
        j = k
        k = k2
    end
    if j == 0
        error("expected expression at $k")
    end
    k = consume_whitespace(xpath, k)
    if iterate(xpath, k) === nothing
        c = '\0'
        k2 = k
    else
        c, k2 = iterate(xpath, k)
    end
    has_fn_last::Bool = false
    #if ismacro
    #    fn::Expr
    #else
    #    fn::SymbolAny
    #end
    if '0' <= xpath[i] <= '9' || xpath[i] == '-'
        # parse token as a number
        num = parse(Float64,xpath[i:j])
        fn = @xpath_fn :number num
        returntype = Number
    elseif xpath[i] == '"' || xpath[i] == '\''
        if ismacro
            str = PipeBuffer()
            sexpr = Expr(:call, :string)
            escape = false
            var = parenvar = false
            substr_k = iterate(xpath, i)[2]
            j = iterate(xpath, j)[2]
            while substr_k != j
                c, substr_k = iterate(xpath,substr_k)
                if var == true
                    if bytesavailable(str) == 0
                        if !parenvar && c == '('
                            parenvar = true
                            continue
                        end
                        if !isletter(c) && c!="_"
                            error("invalid interpolation syntax at $substr_k")
                        end
                        write(str,c)
                        continue
                    elseif !isletter(c) && !isnumeric(c) && c!='_' && c!='!'
                        push!(sexpr.args, Expr(:call,:string,esc(Symbol(String(take!(str))))))
                        var = false
                        if parenvar
                            if c != ')' # we aren't interested in writing a general purpose string parser
                                error("invalid interpolation syntax at $substr_k")
                            end
                            continue
                        end
                    else
                        write(str,c)
                        continue
                    end
                end
                if c == '\\'
                    if escape
                        write(str,c)
                    else
                        escape = true
                    end
                else
                    escape = false
                    if c == '$'
                        var = true
                        parenvar = false
                        bytesavailable(str) != 0 && push!(sexpr.args, String(take!(str)))
                    else
                        write(str,c)
                    end
                end
            end
            if var == true
                (bytesavailable(str) != 0 && !parenvar) || error("invalid interpolation syntax at $j")
                push!(sexpr.args, Expr(:call,:string,esc(Symbol(String(take!(str))))))
            else
                bytesavailable(str) != 0 && push!(sexpr.args, String(take!(str)))
            end
            if length(sexpr.args) == 1
                sexpr = ""
            elseif length(sexpr.args) == 2
                sexpr = sexpr.args[2]
            end
        else
            sexpr = xpath[iterate(xpath,i)[2]:j]
        end
        fn = @xpath_fn :string sexpr
        returntype = AbstractString
    else
        if c == '('
            name = xpath[i:j]
            k, fn_, returntype, has_fn_last = consume_function(xpath, k2, name, ismacro)
        else
            fn_ = nothing
        end
        if fn_ === nothing
            k, fn_, returntype = xpath_parse(xpath, i, ismacro)
            if returntype == Any
                fn_ = @xpath_fn :xpath_any fn_
            elseif returntype == ETree
                fn_ = @xpath_fn :xpath fn_
            elseif returntype == AbstractString
                if !ismacro && length(fn_) == 1 && fn_[1][1]::Symbol == :attribute
                    fn_ = fn_[1]
                else
                    fn_ = @xpath_fn :xpath_str fn_
                end
            else
                @assert(false)
            end
            returntype = Vector{returntype}
        end
        if ismacro
            fn = fn_::Expr
        else
            fn = fn_::SymbolAny
        end
    end
    k = consume_whitespace(xpath, k)
    while iterate(xpath,k) !== nothing
        c1,k1 = iterate(xpath,k)
        if c1 == ']' || c1 == ')' || c1 == ','
            break
        end
        if iterate(xpath,k2) === nothing
            error("unexpected end to xpath")
        end
        c2,k2 = iterate(xpath,k1)
        i = k #backup k
        if c1 == 'o' && c2 == 'r' # lowest precedence (0)
            if iterate(xpath,k2) === nothing
                error("unexpected end to xpath")
            end
            c3,k3 = iterate(xpath,k2)
            if !isspace(c3)
                error("expected a space after operator at $k")
            end
            op_precedence = 0
            op = :or
            k = k3
            returntype = Bool

        elseif c1 == 'a' && c2 == 'n'
            if iterate(xpath,k2) === nothing
                error("unexpected end to xpath")
            end
            c3,k3 = iterate(xpath,k2)
            if c3 != 'd'
                error("invalid operator $c at $k")
            end
            if iterate(xpath,k3) === nothing
                error("unexpected end to xpath")
            end
            c3,k2 = iterate(xpath,k3)
            if !isspace(c3)
                error("expected a space after operator at $k")
            end
            op_precedence = 1
            op = :and
            k = k3
            returntype = Bool

        elseif c1 == '='
            op_precedence = 2
            op = :(=)
            k = k1
            returntype = Bool
        elseif c1 == '!' && c2 == '='
            op_precedence = 2
            op = :(!=)
            k = k2
            returntype = Bool

        elseif c1 == '>'
            op_precedence = 3
            if c2 == '='
                op = :(>=)
                k = k2
            else
                op = :(>)
                k = k1
            end
            returntype = Bool
        elseif c1 == '<'
            op_precedence = 3
            if c2 == '='
                op = :(<=)
                k = k2
            else
                op = :(<)
                k = k1
            end
            returntype = Bool

        elseif c1 == '+'
            op_precedence = 4
            op = :(+)
            k = k1
            returntype = Number
        elseif c1 == '-'
            op_precedence = 4
            op = :(-)
            k = k1
            returntype = Number

        else # highest precedence (5)
            if iterate(xpath,k2) === nothing
                error("unexpected end to xpath")
            end
            c3,k3 = iterate(xpath,k2)
            if iterate(xpath,k3) === nothing
                error("unexpected end to xpath")
            end
            op_precedence = 5
            if c1 == 'd' && c2 == 'i' && c3 == 'v'
                op = :div
            elseif c1 == 'm' && c2 == 'o' && c3 == 'd'
                op = :mod
            else
                error("invalid operator $c1 at $k")
            end
            c4,k = iterate(xpath,k3)
            if !isspace(c4)
                error("expected a space after operator at $k")
            end
            returntype = Number
        end
        if precedence > op_precedence
            k = i #restore k
            break
        end
        k, fn2, rt2, has_fn_last2 = xpath_parse_expr(xpath, k, op_precedence+1, ismacro)
        k = consume_whitespace(xpath, k)

        if ismacro
            fn = Expr(:tuple,:(:binop),Expr(:tuple,Expr(:quote,op), fn, fn2))
        else
            fn = (:binop, (op, fn, fn2))
        end
        has_fn_last |= has_fn_last2
    end
    return k, fn, returntype, has_fn_last
end

function consume_function(xpath, k, name, ismacro)
    #consume a function call
    k = consume_whitespace(xpath, k)
    if iterate(xpath,k) === nothing
        error("unexpected end to xpath after (")
    end
    fntype = get(xpath_functions, name, nothing)
    if fntype === nothing
        return k, nothing, Nothing, false
    end
    minargs = fntype[2]::Int
    maxargs = fntype[3]::Int
    fnreturntype = fntype[4]::DataType
    if ismacro
        args = :(SymbolAny[])
    else
        args = SymbolAny[]
    end

    c, k2 = iterate(xpath,k)
    if c == ','
        error("unexpected , in functions args at $k")
    end
    has_fn_last::Bool = (fntype[1] == :last)
    len_args = 0
    while c != ')'
        k, arg, returntype, has_fn_last2 = xpath_parse_expr(xpath, k, 0, ismacro)
        if ismacro
            args = Expr(:call, :push!, args, arg)
        else
            push!(args, arg)
        end
        len_args += 1
        has_fn_last |= has_fn_last2
        k = consume_whitespace(xpath, k)
        if iterate(xpath,k) === nothing
            error("unexpected end to xpath after (")
        end
        c, k2 = iterate(xpath, k)
        if c != ',' && c != ')'
            error("unexpected character $c at $k")
        end
        k = k2
    end
    if !(minargs <= len_args <= maxargs)
        error("incorrect number of arguments for function $name (found $(length(args)))")
    end
    if ismacro
        fn = Expr(:tuple, Expr(:quote, fntype[1]::Symbol), args)
    else
        fn = (fntype[1]::Symbol, args)
    end
    return k2, fn, fnreturntype, has_fn_last
end

isroot(pd::ETree) = (pd.parent == pd)

struct XPath{T<:AbstractString,
                returntype <: Union{Vector{ETree},
                      Vector{AbstractString},
                      Vector{Any},
                      Bool,
                      Number,
                      Int,
                      AbstractString,
                      Any}}
    # an XPath filter is a series of XPath segments implemented as
    # (:cmd, data) pairs. For example,
    # "//A/..//*[2]" should be parsed as:
    # [(:root,:element), (:descendant_or_self,:node), (:child,:element), (:name,"A")),
    #  (:parent,:element), (:descendant_or_self,:node), (:child,:element), (:filter,(:number,2))]
    # All data strings are expected to be of type T
    filter::SymbolAny
end

mutable struct XPath_Collector
    nodes::Vector{ETree}
    filter::Any
    index::Int
    function XPath_Collector()
        new(ETree[], nothing, 0)
    end
end

xpath(filter::T) where {T<:AbstractString} = (xp = xpath_parse(filter); XPath{T,xp[2]}(xp[1]))

function xpath(pd::Vector, xp::XPath{T,Vector{returntype}}) where {T,returntype}
    output = returntype[]
    for ele in pd
        add = xpath_expr(ele, xp, xp.filter, 1, -1, Vector{returntype})::Vector{returntype}
        output = append!(output, setdiff(add, output))
    end
    return output::Vector{returntype}
end
function xpath(pd::Vector, xp::XPath{T,returntype}) where {T,returntype}
    output = returntype[]
    for ele in pd
        push!(output, xpath_expr(pd, xp, xp.filter, 1, -1, returntype)::returntype)
    end
    return output
end
xpath(pd, xp::XPath{T,returntype}) where {T,returntype} = xpath_expr(pd, xp, xp.filter, 1, -1, returntype)::returntype
xpath(pd, filter::T) where {T<:AbstractString} = xpath(pd, xpath(filter))

function xpath_combined_checked(pd1::XPath, pd2::XPath)
    a1 = pd1.filter[1]
    a2 = pd2.filter[1]
    if a1 == a2
        filt = a1
    else
        if !(( a1 == :xpath || a1 == :xpath_str || a1 == :xpath_any) &&
            (a2 == :xpath || a2 == :xpath_str || a2 == :xpath_any))
            error("can't combine xpath expressions that don't return XML nodes")
        end
        filt = :xpath_any
    end
    xp = SymbolAny[]
    push!(xp, (:(|),(pd1.filter[2], pd2.filter[2])))
    return (filt, xp)
end
|(pd1::XPath{T}, pd2::XPath{S}) where {T,S} =
    XPath{Union{T,S},Any}( xpath_combined_checked(pd1,pd2) )
|(pd1::XPath{T,ret1}, pd2::XPath{S,ret2}) where {T,S,ret1<:Vector,ret2<:Vector} =
    XPath{Union{T,S},Vector{Any}}( xpath_combined_checked(pd1,pd2) )
|(pd1::XPath{T,ret}, pd2::XPath{S,ret}) where {T,S,ret<:Vector} =
    XPath{Union{T,S},ret}( xpath_combined_checked(pd1,pd2) )
|(pd1::XPath{T,ret}, pd2::XPath{S,ret}) where {T,S,ret} =
    XPath{Union{T,S},ret}( xpath_combined_checked(pd1,pd2) )
#*{T,S,ret}(pd::XPath{T,ret}, filters::S) = XPath{Union{T,S},ret}( ??? )

xpath_boolean(a::Bool) = a
xpath_boolean(a::Int) = a != 0
xpath_boolean(a::Float64) = a != 0 && !isnan(a)
xpath_boolean(a::AbstractString) = !isempty(a)
xpath_boolean(a::Vector) = !isempty(a)
xpath_boolean(a::ETree) = true

xpath_number(a::Bool) = a ? 1 : 0
xpath_number(a::Int) = a
xpath_number(a::Float64) = a
xpath_number(a::AbstractString) = try parse(Float64,a) catch ex NaN end
xpath_number(a::Vector) = xpath_number(xpath_string(a))
xpath_number(a::ETree) = xpath_number(xpath_string(a))

xpath_string(a::Bool) = string(a)
xpath_string(a::Int) = string(a)
function xpath_string(a::Float64)
    if a == 0
        return "0"
    elseif isinf(a)
        return (a < 0 ? "-Infinity" : "Infinity")
    elseif isinteger(a)
        return string(Int(a))
    else
        return string(a)
    end
end
xpath_string(a::AbstractString) = a
xpath_string(a::Vector) = length(a) == 0 ? "" : xpath_string(a[1])
xpath_string(a::ETree) = string_value(a)

function xpath_normalize(s::AbstractString)
    normal = IOBuffer()
    space = false
    first = false
    for c in s
        if isspace(c)
            if !space && first
                space = true
            end
        else
            if space
                write(normal,' ')
                space = false
            end
            if !first
                first = true
            end
            write(normal,c)
        end
    end
    String(take!(normal))
end

function xpath_translate(a::AbstractString,b::AbstractString,c::AbstractString)
    b = collect(b)
    c = collect(c)
    tr = IOBuffer()
    for ch in a
        i = findfirst(b,ch)
        if i == 0
            write(tr, ch)
        elseif i <= length(c)
            write(tr, c[i])
        end
    end
    String(take!(tr))
end

function xpath_expr(pd, xp::XPath{T}, filter::Tuple{Symbol,Any}, position::Int, last::Int, output_hint::DataType) where T<:AbstractString
    op = filter[1]::Symbol
    args = filter[2]
    if op == :attribute
        if !isa(pd, ETree)
            return AbstractString[]
        elseif isa(args, Nothing)
            return pd.attr
        else
            attr = get(pd.attr, args::T, nothing)
            if attr === nothing
                return AbstractString[]
            else
                return AbstractString[attr]
            end
        end
    elseif op == :number
        return args::Number
    elseif op == :string
        return args::AbstractString
    elseif op == :position
        @assert(position > 0)
        return position
    elseif op == :last
        @assert(last >= 0)
        return last
    elseif op == :count
        result = xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Vector)::Vector
        return length(result)
    elseif op == :not
        return !(xpath_boolean(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Bool))::Bool)
    elseif op == :true_
        return true
    elseif op == :false_
        return false
    elseif op == :bool
        return xpath_boolean(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Bool))::Bool
    elseif op == :binop
        op = args[1]::Symbol
        if op == :and
            a = xpath_boolean(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Bool))::Bool
            if a
                return xpath_boolean(xpath_expr(pd, xp, args[3]::SymbolAny, position, last, Bool))::Bool
            end
            return false
        elseif op == :or
            a = xpath_boolean(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Bool))::Bool
            if a
                return true
            end
            return xpath_boolean(xpath_expr(pd, xp, args[3]::SymbolAny, position, last, Bool))::Bool
        end
        a = xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Any)
        b = xpath_expr(pd, xp, args[3]::SymbolAny, position, last, Any)
        if op == :(+)
            return xpath_number(a) + xpath_number(b)
        elseif op == :(-)
            return xpath_number(a) - xpath_number(b)
        elseif op == :div
            return xpath_number(a) / xpath_number(b)
        elseif op == :mod
            return xpath_number(a) % xpath_number(b)
        else
            if !isa(a,Vector)
                a = (a,)
            end
            if !isa(b,Vector)
                b = (b,)
            end
            for a = a
                for b = b
                    if isa(a, ETree)
                        if isa(b, ETree)
                            #nothing
                        elseif isa(b, Int) || isa(b, Float64)
                            a = xpath_number(a)
                        elseif isa(b, Bool)
                            a = xpath_boolean(a)
                        elseif isa(b, AbstractString)
                            a = xpath_string(a)
                        else
                            @assert(false)
                        end
                    elseif isa(b, ETree)
                        if isa(a, Int) || isa(a, Float64)
                            b = xpath_number(b)
                        elseif isa(a, Bool)
                            b = xpath_boolean(b)
                        elseif isa(a, AbstractString)
                            b = xpath_string(b)
                        else
                            @assert(false)
                        end
                    end #if
                    if op == :(=) || op == :(!=)
                        if isa(a, Bool) || isa(b, Bool)
                            a = xpath_boolean(a)
                            b = xpath_boolean(b)
                        elseif isa(a, Int) || isa(b, Int)
                            a = xpath_number(a)
                            b = xpath_number(b)
                        else
                            a = xpath_string(a)
                            b = xpath_string(b)
                        end
                        if op == :(=)
                            if a == b
                                return true
                            end
                        else
                            if a != b
                                return true
                            end
                        end
                    else # op != :(=) && op != :(!=)
                        a = xpath_number(a)
                        b = xpath_number(b)
                        if op == :(>)
                            if a > b
                                return true
                            end
                        elseif op == :(>=)
                            if a >= b
                                return true
                            end
                        elseif op == :(<)
                            if a < b
                                return true
                            end
                        elseif op == :(<=)
                            if a <= b
                                return true
                            end
                        else
                            @assert(false)
                        end
                    end #if
                end #for b
            end #for a
            return false
        end #if
    elseif op == :xpath
        if output_hint == Bool
            return xpath(pd, :node, xp, args::Vector{SymbolAny}, 1, Int[], 1, XPath_Collector(), Bool)::Bool
        elseif output_hint == Vector{ETree} || output_hint == Vector || output_hint == Any
            out = ETree[]
            xpath(pd, :node, xp, args::Vector{SymbolAny}, 1, Int[], 1, XPath_Collector(), out)
            return out
        else
            @assert(false, "unexpected output hint $output_hint")
        end
    elseif op == :xpath_str
        if output_hint == Bool
            return xpath(pd, :node, xp, args::Vector{SymbolAny}, 1, Int[], 1, XPath_Collector(), Bool)::Bool
        elseif output_hint == Vector{AbstractString} || output_hint == Vector || output_hint == Any
            out = AbstractString[]
            xpath(pd, :node, xp, args::Vector{SymbolAny}, 1, Int[], 1, XPath_Collector(), out)
            return out
        else
            @assert(false)
        end
    elseif op == :xpath_any
        if output_hint == Bool
            return xpath(pd, :node, xp, args::Vector{SymbolAny}, 1, Int[], 1, XPath_Collector(), Bool)::Bool
        elseif output_hint == Vector{Any} || output_hint == Vector || output_hint == Any
            out = Any[]
            xpath(pd, :node, xp, args::Vector{SymbolAny}, 1, Int[], 1, XPath_Collector(), out)
            return out
        else
            @assert(false, "unexpected output hint $output_hint")
        end
    elseif op == :string_fn
        if length(args) == 0
            a = xpath_string(pd)::AbstractString
        else
            a = xpath_string(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::AbstractString
        end
        return a
    elseif op == :contains
        a = xpath_string(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::AbstractString
        b = xpath_string(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Any))::AbstractString
        return !(isempty(something(findfirst(b, a), 0:-1)))::Bool
    elseif op == :startswith
        a = xpath_string(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::AbstractString
        b = xpath_string(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Any))::AbstractString
        return startswith(a,b)::Bool
    elseif op == :name
        if !isempty(args)
            a = xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Vector)::Vector
            if isempty(a)
                return ""
            else
                return xpath_string(a[1])::AbstractString
            end
        else
            return xpath_string(pd)::AbstractString
        end
    elseif op == :concat
        a = IOBuffer()
        for arg = args
            write(a, xpath_string(xpath_expr(pd, xp, arg::SymbolAny, position, last, Any))::AbstractString)
        end
        return String(take!(a))::AbstractString
    elseif op == :substring_before
        a = xpath_string(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::AbstractString
        b = xpath_string(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Any))::AbstractString
        i = first(something(findfirst(b, a), 0:-1))
        if i < 1
            return ""
        else
            return a[1:prevind(a,i)]::AbstractString
        end
    elseif op == :substring_after
        a = xpath_string(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::AbstractString
        b = xpath_string(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Any))::AbstractString
        i = last(something(findfirst(b, a), 0:-1))
        if i < 1
            return ""
        else
            return a[nextind(a,i):end]::AbstractString
        end
    elseif op == :substring
        a = xpath_string(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::AbstractString
        b = parse(Int, xpath_number(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Any)))::Int
        if length(args) > 2
            c = parse(Int, xpath_number(xpath_expr(pd, xp, args[3]::SymbolAny, position, last, Any)))::Int
            return a[chr2ind(a,b):chr2ind(a,c)]
        else
            return a[chr2ind(a,b):end]
        end
    elseif op == :string_length
        if isempty(args)
            a = xpath_string(pd)::AbstractString
        else
            a = xpath_string((xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any)))::AbstractString
        end
        return length(a)::Int
    elseif op == :normalize_space
        if isempty(args)
            a = xpath_string(pd)::AbstractString
        else
            a = xpath_string((xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any)))::AbstractString
        end
        return xpath_normalize(a)::AbstractString
    elseif op == :translate
        a = xpath_string(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::AbstractString
        b = xpath_string(xpath_expr(pd, xp, args[2]::SymbolAny, position, last, Any))::AbstractString
        c = xpath_string(xpath_expr(pd, xp, args[3]::SymbolAny, position, last, Any))::AbstractString
        return xpath_translate(a,b,c)::AbstractString
    elseif op == :number_fn
        if isempty(args)
            return xpath_number(pd)::Number
        else
            return xpath_number(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::Number
        end
    elseif op == :sum
        a = 0.0
        for n in xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Vector)::Vector
            a += xpath_number(n)
        end
        return a::Float64
    elseif op == :floor
        a = ifloor(xpath_number(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::Number)
        return a::Int
    elseif op == :ceiling
        a = iceil(xpath_number(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::Number)
        return a::Int
    elseif op == :round
        a = xpath_number(xpath_expr(pd, xp, args[1]::SymbolAny, position, last, Any))::Number
        if -0.5 <= a < -0.0 || a === -0.0
            return -0.0
        else
            return floor(a+0.5)
        end
    else
        error("invalid or unimplmented op $op")
    end
    @assert(false)
end

function xpath_output(pd::ETree, output)
    if isa(output,Vector)
        if !in(pd, output)
            push!(output, pd)
        end
    end
    xpath_boolean(pd)::Bool
end
function xpath_output(string::AbstractString, output)
    if isa(output,Vector)
        if !in(string, output)
            push!(output, string)
        end
    end
    xpath_boolean(string)::Bool
end
function xpath_output(strings::Vector{AbstractString}, output)
    if isa(output,Vector)
        for string in strings
            if !in(string, output)
                push!(output, string)
            end
        end
    end
    xpath_boolean(strings)::Bool
end

macro xpath(node)
    esc( :( iscounted |= xpath($(node), name::Symbol, xp, filter, index, position, position_index, collector, output) ) )
end
macro xpath_descendant(node)
    esc( :( iscounted |= xpath_descendant($(node), name::Symbol, xp, filter, index, position, position_index, collector, output) ) )
end
function xpath(pd, nodetype_filter::Symbol, xp::XPath{T}, filter::Vector{SymbolAny}, index::Int, position::Vector{Int}, position_index::Int, collector::XPath_Collector, output) where T<:AbstractString
    #return value is whether the node is "true" for the input to a boolean function
    #implements axes: child, descendant, parent, ancestor, self, root, descendant-or-self, ancestor-or-self
    #implements filters: name
    if nodetype_filter === :element
        if !isa(pd, ETree)
            return false
        end
    elseif nodetype_filter === :comment
        error("xpath comments are not currently saved as part of the xml")
    elseif nodetype_filter === :text
        if !isa(pd, AbstractString)
            return false
        end
    else
        @assert(nodetype_filter === :node)
    end
    if index > length(filter)
        return xpath_output(pd, output)
    end
    axis = filter[index][1]::Symbol
    name = filter[index][2]
    index += 1
    iscounted::Bool = false

    # FILTERS
    if axis == :filter
        s = length(position)+1
        if s <= position_index
            resize!(position, position_index)
            for i = s:position_index-1
                position[s] = -1
            end
            position[end] = 1
        else
            position[position_index] += 1
        end
        p = position[position_index]
        bool = xpath_expr(pd, xp, name::SymbolAny, p, -1, Bool)
        if isa(bool, Int)
            iscounted = bool::Int == p
        elseif isa(bool, Float64)
            iscounted = bool::Float64 == p
        elseif isa(bool, Vector{ETree})
            iscounted = length(bool::Vector{ETree}) != 0
        else
            iscounted = xpath_boolean(bool)::Bool
        end
        if iscounted
            position_index += 1
            iscounted = false
            name = :node
            @xpath pd
        end

    elseif axis == :filter_with_last
        if collector.filter === nothing
            @assert(collector.index == 0)
            collector.nodes = ETree[]
            collector.filter = name::SymbolAny
            collector.index = index
        else
            @assert(collector.filter === name)
            @assert(collector.index === index)
        end
        push!(collector.nodes, pd)
        iscounted = false

    elseif axis == :attribute
        attrs = xpath_expr(pd, xp, filter[index-1]::SymbolAny, -1, -1, Vector{AbstractString})::Vector{AbstractString}
        name = :node
        @xpath attrs

    elseif axis == :name
        if name::T == pd.name
            name = :node
            @xpath pd
        else
            iscounted = false
        end

    elseif axis == :(|)
        filter1 = name[1]::Vector{SymbolAny}
        filter2 = name[2]::Vector{SymbolAny}
        index = 1
        name = :node
        filter = filter1
        @xpath pd
        filter = filter2
        @xpath pd

    # AXES
    else
        if !isa(pd, ETree)
            return false
        end
        if axis == :root
            root = pd
            while !isroot(root)
                root = root.parent
            end
            @xpath root
        elseif axis == :parent
            @xpath pd.parent
        elseif axis == :ancestor
            parent = pd
            while !isroot(parent)
                parent = parent.parent
                @xpath parent
            end
        elseif axis == :ancestor_or_self
            parent = pd
            @xpath parent
            while !isroot(parent)
                parent = parent.parent
                @xpath parent
            end
        elseif axis == :self
            @xpath pd
        elseif axis == :child
            if name::Symbol == :node
                for (key,attr) in pd.attr
                    @xpath attr
                end
            end
            for child in pd.elements
                if isa(child, ETree)
                    @xpath child::ETree
                elseif isa(child, AbstractString)
                    @xpath child::AbstractString
                else
                    @assert(false)
                end
            end
        elseif axis == :descendant
            @xpath_descendant pd
        elseif axis == :descendant_or_self
            @xpath pd
            @xpath_descendant pd

        #TODO: more axes
        #elseif axis == :attribute
        #elseif axis == :following
        #elseif axis == :following-sibling
        #elseif axis == :preceding
        #elseif axis == :preceding-sibling

        #TODO: axis in xpath_types

        # ERROR - NO MATCH
        else
            error("encountered unsupported axis $axis")
        end
        while collector.filter !== nothing
            nodes = collector.nodes
            collector_filter = collector.filter::SymbolAny
            index = collector.index
            name = :node
            collector.filter = nothing
            collector.index = 0
            last = length(nodes)
            count = 1
            for pd = nodes
                if xpath_boolean(xpath_expr(pd, xp, collector_filter, count, last, Bool))
                    @xpath pd
                end
                count += 1
            end
            empty!(nodes)
        end
        for i = position_index:length(position)
            if position[i] > 0
                position[i] = 0
            end
        end
    end
    return iscounted
end

function xpath_descendant(pd::ETree, name::Symbol, xp::XPath, filter::Vector{SymbolAny}, index::Int, position::Vector{Int}, position_index::Int, collector::XPath_Collector, output)
    iscounted = false
    for child in pd.elements
        if name::Symbol == :node
            for (key,attr) in pd.attr
                @xpath attr
            end
        end
        if isa(child, ETree)
            @xpath child::ETree
            @xpath_descendant child::ETree
        elseif isa(child, AbstractString)
            @xpath child::AbstractString
        else
            @assert(false)
        end
    end
    iscounted
end

getindex(pd::ETree,x::AbstractString) = xpath(pd,x)
getindex(pd::ETree,x::XPath) = xpath(pd,x)
getindex(pd::Vector{ETree},x::AbstractString) = xpath(pd, x)
getindex(pd::Vector{ETree},x::XPath) = xpath(pd, x)


text:
export Pass

# subtypes are expected to have a 'ref::API.LLVMPassRef' field
abstract type Pass end

Base.unsafe_convert(::Type{API.LLVMPassRef}, pass::Pass) = pass.ref

#
# Module passes
#

export ModulePass

if VERSION >= v"1.6.0-DEV.809"

function module_pass_callback(ptr, data)
    mod = Module(convert(API.LLVMModuleRef, ptr))
    runner_box = Base.unsafe_pointer_to_objref(data)
    runner = runner_box[]
    changed = runner(mod)::Core.Bool
    convert(Bool, changed)
end

@checked struct ModulePass <: Pass
    ref::API.LLVMPassRef
    root::Any

    function ModulePass(name::String, runner::Core.Function)
        callback = @cfunction(module_pass_callback, API.LLVMBool, (Ptr{Cvoid}, Ptr{Cvoid}))
        runner_box = Ref(runner)

        ref = API.LLVMCreateModulePass2(name, callback, runner_box)
        refcheck(ModulePass, ref)

        return new(ref, runner_box)
    end
end

else

@checked struct ModulePass <: Pass
    ref::API.LLVMPassRef
    root::Any

    function ModulePass(name::String, runner::Core.Function)
        function callback(ptr::Ptr{Cvoid})::Core.Bool
            mod = Module(convert(API.LLVMModuleRef, ptr))
            return runner(mod)::Core.Bool
        end

        ref = API.LLVMCreateModulePass(name, callback)
        refcheck(ModulePass, ref)

        return new(ref, callback)
    end
end

end

#
# Function passes
#

export FunctionPass

if VERSION >= v"1.6.0-DEV.809"

function function_pass_callback(ptr, data)
    fn = Function(convert(API.LLVMValueRef, ptr))
    runner_box = Base.unsafe_pointer_to_objref(data)
    runner = runner_box[]
    changed = runner(fn)::Core.Bool
    convert(Bool, changed)
end

@checked struct FunctionPass <: Pass
    ref::API.LLVMPassRef
    root::Any

    function FunctionPass(name::String, runner::Core.Function)
        callback = @cfunction(function_pass_callback, API.LLVMBool, (Ptr{Cvoid}, Ptr{Cvoid}))
        runner_box = Ref(runner)

        ref = API.LLVMCreateFunctionPass2(name, callback, runner_box)
        refcheck(FunctionPass, ref)

        return new(ref, runner_box)
    end
end

else

@checked struct FunctionPass <: Pass
    ref::API.LLVMPassRef
    root::Any

    function FunctionPass(name::String, runner::Core.Function)
        function callback(ptr::Ptr{Cvoid})::Core.Bool
            fn = Function(convert(API.LLVMValueRef, ptr))
            return runner(fn)::Core.Bool
        end

        ref = API.LLVMCreateFunctionPass(name, callback)
        refcheck(FunctionPass, ref)

        return new(ref, callback)
    end
end

end

#
# Basic block passes
#

export BasicBlockPass

@checked struct BasicBlockPass <: Pass
    ref::API.LLVMPassRef
    root::Any

    function BasicBlockPass(name::String, runner::Core.Function)
        VERSION >= v"1.4.0-DEV.589" && error("BasicBlockPass functionality has been removed from Julia and LLVM")

        function callback(ptr::Ptr{Cvoid})::Core.Bool
            bb = BasicBlock(convert(API.LLVMBasicBlockRef, ptr))
            return runner(bb)::Core.Bool
        end

        ref = API.LLVMCreateBasicBlockPass(name, callback)
        refcheck(BasicBlockPass, ref)

        return new(ref, callback)
    end
end


text:
function get_lr(l, r, acc_as)
    if l == 1
        return acc_as[r]
    else
        return acc_as[r] - acc_as[l-1]
    end
end

function dp(l, r, as, acc_as, memo)
    if memo[l, r] > 0
        return memo[l, r]
    elseif l == r
        return memo[l, r] = 0
    else
        res = minimum([dp(l, i, as, acc_as, memo) +
            dp(i+1, r, as, acc_as, memo) +
            get_lr(l, r, acc_as) for i in l:r-1])
        return memo[l, r] = res
    end
end

function main()
    N = parse(Int, readline())
    as = parse.(Int, split(readline()))
    acc_as = cumsum(as)

    memo = zeros(Int, N, N)
    println(dp(1, N, as, acc_as, memo))
end

main()


text:
export add_Hopf!
export Hopf_variables
export replace_Hopf_variable!

using HarmonicBalance: is_rearranged, rearrange_standard
using HarmonicBalance: LinearResponse.get_implicit_Jacobian

"Add a Hopf bifurcation to the system."
function add_Hopf!(eom::DifferentialEquation, Δω::Num)
    for var in get_variables(eom), ω in eom.harmonics[var]
        add_harmonic!(eom, var, ω + Δω)
        add_harmonic!(eom, var, ω - Δω)
    end
end

"Successively adds n Hopf bifurcations at the same frequency Δω"
add_Hopf!(eom::DifferentialEquation; frequency::Num, multiplicity::Int) = [add_Hopf!(eom, frequency) for k in 1:multiplicity]

"""
    Hopf_variables(eom, freq)
"""
Hopf_variables(eom::HarmonicEquation, freq::Num) = filter(x -> any(isequal.(freq, get_all_terms(x.ω))), eom.variables)

"""
    Obtain the Jacobian of `eom` with a free (Hopf) variable `free_var`.
    `free_var` marks the variable which is free due to U(1) symmetry. Its entry in the Jacobian matrix must be discarded
    (the discarded eigenvalue is 0, corresponding to a free phase.)
"""
function _Hopf_Jacobian(eom::HarmonicEquation, free_var::HarmonicVariable)
    eom_Jac = rearrange_standard(eom)
    free_idx = findall(x -> isequal(x, free_var), eom_Jac.variables)
    deleteat!(eom_Jac.equations, free_idx)
    deleteat!(eom_Jac.variables, free_idx)

    # the free variable can be fixed to zero, this is also done in the corresponding HarmonicEquation later
    eom_Jac = substitute_all(eom_Jac, free_var => 0)
    J = get_implicit_Jacobian(eom_Jac)
end

"""
    Construct a `Problem` in the case where U(1) symmetry is present
    due to having added a Hopf variable with frequency `Hopf_ω`.
"""
function _Hopf_Problem(eom::HarmonicEquation, Hopf_ω::Num)

    isempty(Hopf_variables(eom, Hopf_ω)) ? error("No Hopf variables found!") : nothing
    !any(isequal.(eom.parameters, Hopf_ω)) ? error(Hopf_ω, " is not a parameter of the harmonic equation!") : nothing

    free_var = Hopf_variables(eom, Hopf_ω)[end] # eliminate one of the Cartesian variables, it does not matter which

    # get the Hopf Jacobian before altering anything - this is the usual Jacobian but the entries corresponding
    # to the free variable are removed
    J = _Hopf_Jacobian(eom, free_var)

    new_symbol = HarmonicBalance.declare_variable(string(Hopf_ω), first(get_independent_variables(eom)))
    rules = Dict(Hopf_ω => new_symbol, free_var.symbol => Num(0))
    eom.equations = expand_derivatives.(substitute_all(eom.equations, rules))
    eom.parameters = setdiff(eom.parameters, [Hopf_ω]) # Hopf_ω is now NOT a parameter anymore

    free_var.type = "Hopf"
    free_var.ω = Num(0)
    free_var.symbol = new_symbol
    free_var.name = var_name(new_symbol)

    # define Problem as usual but with the Hopf Jacobian (always computed implicitly)
    p = Problem(eom; Jacobian=J)
    return p
end


text:
#=

Define all the scientific constants needed in this module here

=#

const c₂                 = 1.4387769
const cMassMol           = 1.66053873e-27
const cSqrtLn2divSqrtPi  = 0.469718639319144059835
const cLn2               = 0.6931471805599
const cSqrtLn2           = 0.8325546111577
const cSqrt2Ln2          = 1.1774100225
const cc_                = 2.99792458e8
const cBolts_            = 1.3806503e-23
const p_ref              = 1013.25  # reference pressure [hPa]
const t_ref              = 296.0    # reference temperature [K]
const nm_per_m           = 1.0e7


text:
using Images
using Javis
using LaTeXStrings
using ReferenceTests
using Test
using VideoIO

@testset "Unit" begin
    include("unit.jl")
end
@testset "SVG LaTeX tests" begin
    include("svg.jl")
end
@testset "Animations" begin
    include("animations.jl")
end


text:
using GeometricAlgebra
using GeometricAlgebra: StaticArrays.SVector

@testset "constructors" begin
	@test iszero(zero(Blade{(1,1,1),2,0b101,Float64}))
	@test iszero(zero(Blade{(1,1,1)}(big(42), 0b000)))
	@test iszero(zero(Multivector{(1,1,1),1,Vector{Int}}))
	@test iszero(zero(MixedMultivector{(1,1,1),Vector{Int}}))

	@test first(basis((1,1,1))) == Blade{(1,1,1)}(1, 0b001)
	@test  last(basis((1,1,1))) == Blade{(1,1,1)}(1, 0b100)
end

@testset "equality" begin
	a = Blade{(1,1)}(42, 0b01)
	@test a == Blade{(1,1)}(42, 0b01)
	@test a == Blade{(1,1)}(42.0, 0b01)
	@test a == Blade{(1,1)}(BigFloat(42), 0b01)
	@test a != Blade{(1,1)}(-1, 0b01)
	@test a != Blade{(1,1)}(42, 0b10)
	@test a != Blade{(1,1)}(42, 0b11)
	@test a != Blade{(1,1,1)}(42, 0b01)

	a = Multivector{(1,1,1),1}([3, 5, 7])
	@test a == Multivector{(1,1,1),1}(Float16[3, 5, 7])
	@test a == Multivector{(1,1,1),1}(Number[3, 5, 7 + 0im])
	@test a != Multivector{(1,1,1),1}([3, -1, 7])
	@test a != Multivector{(1,1,1),2}([3, 5, 7])
	@test a != Multivector{(-1,1,1),1}([3, 5, 7])

	a = MixedMultivector{(1,1)}([1, 10, 20, 100])
	@test a == MixedMultivector{(1,1)}(BigInt[1, 10, 20, 100])
	@test a != MixedMultivector{(1,1)}([1, 10, -1, 100])
	@test a != MixedMultivector{(-1,1)}([1, 10, 20, 100])

	# equal if both zero even if grade differs
	@test Blade{(1,1)}(0, 0b01) == Blade{(1,1)}(0, 0b10)
	@test Multivector{(1,1),1}([0, 0]) == Multivector{(1,1),0}([0])
	@test Blade{(1,1)}(0, 0b01) != Blade{(1,1)}(1, 0b10)

	@test Blade{(1,1)}(0, 0b01) == 0
	@test Blade{(1,1)}(1, 0b00) == 1
	@test Blade{(1,1)}(2, 0b00) != 1
	@test Blade{(1,1)}(1, 0b01) != 1
end

# TODO: test getcomponent and Base.getindex

@testset "best_type()" begin
	x = Blade{(1,1),1,0b01,Int}(42)

	@testset "given $label" for (label, a) ∈ ["intances" => x, "types" => typeof(x)]

		@test best_type(Blade, a) === Blade{(1,1),k,bits,Int} where {k,bits}
		@test best_type(Multivector, a) === Multivector{(1,1),k,Vector{Int}} where k
		@test best_type(MixedMultivector, a) === MixedMultivector{(1,1),Vector{Int}}

		@test best_type(Blade, a; bits=Val(0b11)) === Blade{(1,1),2,0b11,Int}
		@test best_type(Blade, a; promote_eltype_with=BigInt) === Blade{(1,1),k,bits,BigInt} where {k,bits}
		@test best_type(Blade, a; bits=Val(0b00), promote_eltype_with=Bool) === Blade{(1,1),0,0b00,Int}
		@test best_type(Multivector, a; grade=Val(2)) === Multivector{(1,1),2,Vector{Int}}
		@test best_type(Multivector, a; promote_eltype_with=Complex{Int}) === Multivector{(1,1),k,Vector{Complex{Int}}} where k
		@test best_type(Multivector, a; grade=Val(1), promote_eltype_with=Float32) === Multivector{(1,1),1,Vector{Float32}}
		@test best_type(MixedMultivector, a; promote_eltype_with=Float32) === MixedMultivector{(1,1),Vector{Float32}}

		# single-argument versions preserve parameters
		@test best_type(x) == (x isa Type ? x : typeof(x))
		@test best_type(a; promote_eltype_with=BigFloat, bits=Val(0b0)) === Blade{(1,1),0,0b0,BigFloat}
	end

	y = Multivector{(1,1),1,Vector{Float64}}(0:1)
	@test best_type(y) == typeof(y)
	@test best_type(Multivector, x, y) == Multivector{(1,1),k,Vector{Float64}} where k

	@test_throws Exception best_type(MixedMultivector, MixedMultivector{(1,1,1),Vector{Int}}, Multivector{(1,1),0,Vector{Int}})
end

@testset "conversion" begin
	scalar_types = [
		Int,
		Int128,
		BigInt,
		Float64,
		BigFloat,
		Complex{Int},
		Complex{Float64},
		Rational{Int8},
	]
	@testset "eltype conversions" begin
		blades = Blade{(1,)}.(one.(scalar_types), 0b1)
		@testset "$(eltype(a)) -> $(eltype(b))" for a ∈ blades, b ∈ blades
			@test convert(typeof(b), a) == b
		end
	end

	x, y = basis((1, 1))
	@test float(x) isa Blade{(1, 1),k,bits,<:AbstractFloat} where {k,bits}
	@test big(x) isa Blade{(1, 1),k,bits,BigInt} where {k,bits}
	@test complex(x) isa Blade{(1, 1),k,bits,Complex{Int}} where {k,bits}

	fns = float, big, complex, real
	objs = x*y, x + y, x*y - 1
	@testset "$(GeometricAlgebra.multivectortype(a))" for a ∈ objs
		@testset "$fn()" for fn ∈ fns
			@test typeof(fn(a)) == fn(typeof(a))
			@test typeof(fn(complex(a))) == fn(typeof(complex(a)))
		end
		@test real(a + 2im*a) == a
		@test imag(2a + im*a) == a
	end

	@test Multivector(x)::Multivector == x
	@test MixedMultivector(x)::MixedMultivector == x
	@test MixedMultivector(x + 1) == x + 1
end

@testset "promotion" begin
	x, y, z = basis((1, 1, 1))
	@test ==(promote(2x, x + x)...)
	@test ==(promote(x, x + 0)...)
	@test ==(promote(x + y, x + y + 0)...)

	@test ==(promote(1, x^2)...)
	@test ==(promote(2, (x + y)^2)...)
	@test ==(promote(1, (x + 0)^2)...)

	# objects are promoted when stored in vectors
	@test first([1, x]) isa Blade
	@test first([2, x + y]) isa Multivector
	@test first([3, z + 1]) isa MixedMultivector
	@test sum([1, x, x + y, z + 1]) == 2 + 2x + y + z

	@test prod([x, x*y, y*x*z]) == x*z
end

@testset "isapprox()" begin
	x, y, z = basis((1, 1, 1))
	@test x^2 ≈ 1
	@test sin(1f0)x ≈ sin(1)x
	@test (x + y)^2 ≈ 2
	@test (1 + x*y)^4 ≈ -4

	@test ((x + y)/√2)^Int(1e10) ≈ 1 rtol=1e-4

	# TODO: test zero-equalities
end

@testset "storagetype" begin
	sig = (1, 1, 1)
	b = basis(sig)
	B = convert.(Multivector{sig,1,SVector{k,Int} where k}, b)

	@test b[1]b[2] isa Blade{sig,2}
	@test B[1]B[2] isa MixedMultivector{sig,<:SVector}

	# we want operations with blades to preserve the storage type
	isstatic(a) = GeometricAlgebra.unwrap(GeometricAlgebra.storagetype(a)) <: SVector
	@test isstatic(B[1] + b[2])
	@test isstatic(B[1] - b[2])
	@test isstatic(B[1]b[2])
	@test isstatic(B[1]/b[2])
	@test isstatic((B[1] + 1)∧b[2])
	@test isstatic((B[1] + 1)*b[2])
end


text:
# # GARTEUR SM-AG19 Testbed: Harmonic Vibration Analysis

# Source code: [`garteur_hva_tut.jl`](garteur_hva_tut.jl)

# ## Description

# This virtual test application is based on the test article
# used by the GARTEUR Structures & Materials Action Group 19
# which organized a Round Robin exercise where 12 European laboratories
# tested a single structure between 1995 and 1997. The benchmark structure
# was a laboratory structure built to simulate the dynamic behaviour
# of an aeroplane. The structure was initially built for a benchmark
# study on experimental modal analysis conducted by the
# Structures and Materials Action Group (SM-AG19) of the Group
# for Aeronautical Research and Technology in EURope (GARTEUR).
# The test-bed was designed and manufactured by ONERA, France.

# ![](IMAC97photo.png)

# ### References

# [1] Ground Vibration Test Techniques, compiled by A Gravelle, GARTEUR
# Structures & Materials Action Group 19 Technical report TP-115, 1999.
# [2] Etienne Balmes, Jan R. Wright, GARTEUR GROUP ON GROUND VIBRATION TESTING |
# RESULTS FROM THE TEST OF A SINGLE STRUCTURE BY 12 LABORATORIES IN EUROPE,
# Proceedings of DETC'97, 1997 ASME Design Engineering Technical Conferences,
# September 14-17, 1997, Sacramento, California.
# [3] 3M(TM) Viscoelastic Damping Polymer 112 Series,  Technical Data, May 2017.

# ## Goals

# - Show how to construct model from multiple connected beams.
# - Demonstrate the use of massless connectors.
# - Demonstrate the use of point masses.
# - Demonstrate the use of grounded springs.
# - Illustrate verification of the solution of the free vibration problem.

##
# ## Geometry of the testbed airplane.

# It was a rather simple structure which was reasonably dynamically
# representative of a simple airplane structure. It was composed of several beams
# simulating a fuselage with wings and a tail. Wing tip drums allowed to adjust
# bending and torsion frequencies similarly to airplane ones, with some very
# close modal frequencies.

# ![](garteur-geom.png)

# The script included below defines the geometry of the structure, the
# cross-sectional properties, the connectivity, and the location of the nodes.
include("garteur_geometry_tut.jl")

# The geometry is visualized in the tutorial
# [garteur_geometry_vis_tut](garteur_geometry_vis_tut.jl).

##
# ## Material

# Material properties can be now used to create a material: isotropic elasticity
# model of the `FinEtoolsDeforLinear` package is instantiated.
using FinEtoolsDeforLinear
# The material of the structure is aluminum.
# The elastic modulus:
E = 70000.0 * phun("MPa")
nu = 0.31;
# The mass density:
rho = 2700 * phun("kg/m^3")
alu = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

# Material for the massless connectors has the mass density set to zero;
# otherwise it has the same properties as the aluminum material  of the
# structure.
massless = MatDeforElastIso(DeforModelRed3D, 0.0, alu.E, alu.nu, 0.0)

# This simple function returns material based on the label of the beam elements.
getmaterial(labl) = begin
    if labl >= 7
        return massless
    end
    return alu
end

##
# ## Fields

# Now we start constructing the discrete finite element model.
# We begin by constructing the requisite fields, geometry and displacement.
# These are the so-called "configuration variables", all initialized to 0.
# This is that geometry field.
geom0 = NodalField(fens.xyz)
# This is the displacement field, three unknown displacements per node.
u0 = NodalField(zeros(size(fens.xyz, 1), 3))
# This is the rotation field, three unknown rotations per node are represented
# with a rotation matrix, in total nine numbers. The utility function
# `initial_Rfield`
using FinEtoolsFlexBeams.RotUtilModule: initial_Rfield
Rfield0 = initial_Rfield(fens)

# Finally, this is the displacement and rotation field for incremental changes,
# incremental displacements and incremental rotations. In total, 6 unknowns per
# node. Note that the the incremental displacements are in general complex.
dchi = NodalField(0.0im .* zeros(size(fens.xyz, 1), 6))

# There are no support conditions.
applyebc!(dchi)
# The  the number of free(unknown) degrees of freedom is equal to the total
# number of degrees of freedom in the system.
numberdofs!(dchi);

##
# ## Identify support points and locations of sensors

# Suspension points are at these nodes:
suspln = selectnode(fens; box = initbox!(Float64[], vec([0.0*L 0.0*L 0.805*L])), inflate = tolerance)
susprn = selectnode(fens; box = initbox!(Float64[], vec([0.0*L -0.0*L 0.805*L])), inflate = tolerance)
suspbn = selectnode(fens; box = initbox!(Float64[], vec([-2.0*L 0.0*L 0.0*L])), inflate = tolerance)

# Find out at which nodes the sensors are:
sensors = Dict()
sensors = let
    # The sensors at the tip of the left and right wing drum
    sensor112n = selectnode(fens; box = initbox!(Float64[], vec([+1.8*L 9.8*L 0.96*L])), inflate = tolerance)
    sensors[112] = sensor112n
    sensor12n = selectnode(fens; box = initbox!(Float64[], vec([+1.8*L -9.8*L .96*L])), inflate = tolerance)
    sensors[12] = sensor12n
    sensor111n = selectnode(fens; box = initbox!(Float64[], vec([-1.8*L 9.8*L 0.96*L])), inflate = tolerance)
    sensors[111] = sensor111n
    sensor11n = selectnode(fens; box = initbox!(Float64[], vec([-1.8*L -9.8*L .96*L])), inflate = tolerance)
    sensors[11] = sensor11n
        # The joint between the horizontal and vertical tail parts
    sensor202n = selectnode(fens; box = initbox!(Float64[], vec([-8*L 0 3.8*L])), inflate = tolerance)
    sensors[202] = sensor202n
    sensors
end

##
# ## Assemble the global discrete system

# For disambiguation we will refer to the stiffness and mass functions by
# qualifying them with the corotational-beam module, `FEMMCorotBeamModule`.
using FinEtoolsFlexBeams.FEMMCorotBeamModule
CB = FEMMCorotBeamModule

# Note that we have an array of finite element sets. We compute the matrices for
# each set separately and accumulate them into the final overall matrix. Thus
# we can construct the stiffness and mass matrix as follows.
using  SparseArrays

# Loop over all the finite element sets and add up their contributions, the
# stiffness in the mass matrix for each.
Kf, Kd, M = let
    Kf = spzeros(dchi.nfreedofs, dchi.nfreedofs)
    Kd = spzeros(dchi.nfreedofs, dchi.nfreedofs)
    M = spzeros(dchi.nfreedofs, dchi.nfreedofs)
    for fes in fesa
        labl  = fes.label[1]
        femm = CB.FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), getmaterial(labl));
        if labl == 7 # connectors representing the damping layer
            Kd += CB.stiffness(femm, geom0, u0, Rfield0, dchi);
        else
            Kf += CB.stiffness(femm, geom0, u0, Rfield0, dchi);
        end
        M += CB.mass(femm, geom0, u0, Rfield0, dchi);
    end
    Kf, Kd, M
end

##
# ## Additional concentrated masses.

using LinearAlgebra

using FinEtoolsFlexBeams.FEMMPointMassModule
PM = FEMMPointMassModule

# There is a sensor on the tail.
femmcm1 =  PM.FEMMPointMass(IntegDomain(FESetP1(reshape(sensors[202], 1, 1)), PointRule()), FFltMat(2*L*L/5*L/5*2*rho*LinearAlgebra.I(3)));

# These are the forward/interior locations on the wing drums where the
# compensation masses are attached.
mass1n = selectnode(fens; box = initbox!(Float64[], vec([1.8*L 9.2*L .96*L])), inflate = tolerance)
mass2n = selectnode(fens; box = initbox!(Float64[], vec([1.8*L -9.2*L .96*L])), inflate = tolerance)
femmcm2 =  PM.FEMMPointMass(IntegDomain(FESetP1(reshape([mass1n; mass2n;], 2, 1)), PointRule()), FFltMat(0.2*phun("kg")*LinearAlgebra.I(3)));

Mp = PM.mass(femmcm1, geom0, u0, Rfield0, dchi) + PM.mass(femmcm2, geom0, u0, Rfield0, dchi);

##
# ## Bungee supports

# This is the assumed stifffness of the bungee cords (each one separately).
bungeecoefficient = 4000*phun("N/m");

using LinearAlgebra

using FinEtoolsFlexBeams.FEMMPointGroundedSpringModule
BS = FEMMPointGroundedSpringModule

# There are three suspension points at the top of the fuselage. We assume that these bungee supports exert only reaction in the vertical direction.
femmbs =  BS.FEMMPointGroundedSpring(IntegDomain(FESetP1(reshape([suspln; susprn; suspbn;], 3, 1)), PointRule()),
FFltMat([bungeecoefficient*[0;0;1]*[0;0;1]' 0*LinearAlgebra.I(3); 0*LinearAlgebra.I(3) 0*LinearAlgebra.I(3)]));

Kb = BS.stiffness(femmbs, geom0, u0, Rfield0, dchi)

##
# ## Damping

# The overall damping ratios measured in the physical experiments were
# approximately 1%. Damping levels thought appropriate for obtaining realistic
# damping were obtained through the use of a viscoelastic layer glued to the
# wing beam with an aluminum constraining plate on top.

# ### Viscoelastic damping layer

# The properties of the material of the viscoelastic layer were described in the
# technical specification [3]. Since the properties are frequency dependent, we
# take as a representative value the numbers obtained for 20 Hz.

# The viscoelastic used was the 3M acrylic viscoelastic polymer ISD 112 in the
# form of a tape of 76 mm width with thickness of 50 microm.  This material is
# particularly well suited for the testbed operating range of 5-50 Hz and at 20
# degrees C where the loss factor is near its peak of 0.4 [3]. These are
# representative quantities taken at 20 Hz.

eta = 0.5
omega_5_50 = 2*pi*20

# The connectors between the wing beam and the constraining plate are taken as
# representative of the stiffness of the constraining layer. The viscoelastic
# damping is then taken as proportional to this stiffness.
Cd = (eta / omega_5_50) .* Kd

# ### Aluminum structure

# We assume the loss factor of the fuselage, wing, and tail to be 0.01. This can
# be used to derive the damping model for these parts of the aircraft structure
# in the form of Rayleigh damping.
zeta= 0.01/2; # damping ratio
zeta_1= zeta; # damping ratio for mode 1
zeta_2= zeta; # damping ratio for mode 2
# Stiffness and mass proportional damping parameters.
omega_1=2*pi*6; # Guess
omega_2=2*pi*30; # Guess
rdmass = (2*1.0/(omega_2^2-omega_1^2)).*[omega_2^2 -omega_1^2]*[zeta_1*omega_1; zeta_2*omega_2];
rdstiffness = (2*1.0/(omega_2^2-omega_1^2)).*[-1 1]*[zeta_1*omega_1; zeta_2*omega_2];

# The airframe damping matrix is a mixture of the mass matrix and the stiffness
# matrix.
Cf = rdmass .* M + rdstiffness .* Kf;

# These are the system matrices. The stiffness consists of the contribution of
# the airframe, the connectors between the wing components, and the bungees of
# the suspension.
Kt = Kf + Kd + Kb
# The mass matrix is the contribution of the structure and the attached point
# masses (compensation masses).
Mt = M + Mp
# The damping matrix is composed of the damping of the overall structure, and of
# the damping layer (the connectors between the wing components).
Ct = Cf + Cd

##
# ## Loading

# Here we assume that the stinger  was attached at the location of the sensor
# 12. The force is vertical (in the Z direction) and the magnitude of the force
# is arbitrary.
forceat = 12
fmagn= 1.0;
loadbdry = FESetP1(reshape(sensors[forceat], 1, 1))
lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
# The force is applied in the vertical direction, and we assume it is positive
# upwards.
fi = ForceIntensity(FFlt[0, 0.0, -fmagn, 0, 0, 0]);

# The force loading is now integrated over the "volume" of the integration
# domain.
F = CB.distribloads(lfemm, geom0, dchi, fi, 3);

##
# ## Solve the Harmonic vibration problem

# The frequency sweep will start at `fromf` and continue through the frequency
# `tof`.  The frequencies will be logarithmically distributed throughout this
# range.
fromf = 3.0
tof = 70.0
nf = 150
frequencies = logspace(log10(fromf), log10(tof), nf);

# Only a single number per sensor will be collected, the Z direction
# displacement.
receptance12 = fill(0.0im, nf)
receptance112 = fill(0.0im, nf)

# Now Loop over
for   fi in 1:length(frequencies)
    f =  frequencies[fi];
    om = 2*pi*f;
    # Solve the system of complex  equations of balance:
    U = (-om^2*Mt + 1im*om*Ct + Kt) \ F
    # Distribute the vector of the solution:
    scattersysvec!(dchi, U)
    # Now sample the solution at the locations of the two sensors, 12 and 112:
    v = fill(0.0im, 6)
    gathervalues_asvec!(dchi, v, sensors[12])
    p_d = v[1:3];
    receptance12[fi] = p_d[3]/fmagn;
    gathervalues_asvec!(dchi, v, sensors[112])
    p_d = v[1:3];
    receptance112[fi] = p_d[3]/fmagn;
end

# The receptances were calculated above. The mobility and accelerance may be
# obtained in a postprocessing step.
oms = (2*pi) .* frequencies;
mobility12 = receptance12 .* (-1im*oms);
mobility112 = receptance112 .* (-1im*oms);
accelerance12 = receptance12 .* (-oms.^2);
accelerance112 = receptance112 .* (-oms.^2);

results = Dict()
results[12] = Dict("receptance"=>receptance12, "mobility"=>mobility12, "accelerance"=>accelerance12)
results[112] = Dict("receptance"=>receptance112, "mobility"=>mobility112, "accelerance"=>accelerance112)

##
# ## Present the results graphically

using Gnuplot

# Plot the amplitude of the response curves. We output two curves.
# The first for the driving-point FRF:
quantity = "accelerance"; units = "m/s^2/N"
outputat = 12
y = abs.(results[outputat][quantity]) / phun(units)
@gp  "set terminal windows 0 "  :-
@gp  :- frequencies y " lw 2 lc rgb 'blue' with lines title 'output at $(outputat)' "  :-

# The second for the cross transfer:
outputat = 112
y = abs.(results[outputat][quantity]) / phun(units)
@gp  :- frequencies y " lw 2 lc rgb 'red' with lines title 'output at $(outputat)' "  :-
@gp  :- "set logscale y" :-
@gp  :- "set xlabel 'Frequency [Hz]'" :-
@gp  :- "set ylabel 'abs(H) [$(units)]'" :-
@gp  :- "set title 'Force at $(forceat), $(quantity)'"

# Plot the phase shift of the response curves. Again we output two curves,
# the first for the driving-point FRF:
outputat = 12
y = atan.(imag(results[outputat][quantity]), real(results[outputat][quantity]))/pi*180
@gp  "set terminal windows 1 "  :-
@gp  :- frequencies y " lw 2 lc rgb 'blue' with lines title 'output at $(outputat)' "  :-
# The second for the cross transfer:
outputat = 112
y = atan.(imag(results[outputat][quantity]), real(results[outputat][quantity]))/pi*180
@gp  :- frequencies y " lw 2 lc rgb 'red' with lines title 'output at $(outputat)' "  :-
@gp  :- "set xlabel 'Frequency [Hz]'" :-
@gp  :- "set ylabel 'Phase shift [deg]'" :-
@gp  :- "set title 'Force at $(forceat), $(quantity)'"

true


text:
module OmniSci

using Thrift, LibGEOS

import Thrift: process, meta, distribute, ThriftMeta, ThriftMetaAttribs
import Base: show, convert
import DataFrames: DataFrame, rename!
import Random: randstring
import Dates: DateTime, Date, Time, TimeType, datetime2unix, Second, unix2datetime, hour, minute, second
import Sockets: connect
import LibGEOS: readgeom, writegeom
import GeoInterface:  Point, LineString, Polygon, MultiPolygon, AbstractPoint, AbstractLineString, AbstractPolygon, AbstractMultiPolygon
import DecFP: Dec32, Dec64, Dec128
import Decimals: Decimal
import Tables: eachcolumn, rows, schema, columns

#### types and enums ####
export TDatum, TColumn, TExecuteMode

#### functions ####
export connect, disconnect
export get_status, get_hardware_info
export get_tables_meta, get_table_details
export get_users, get_databases
export sql_execute, sql_execute_df, sql_execute_gdf
export render_vega
export get_dashboards, get_dashboard_grantees
export load_table_binary_columnar, load_table_binary_arrow, load_table
export create_table, import_table
export get_roles, get_all_roles_for_user
export DataFrame

#Thrift interface code, mostly generated from Thrift.jl
include("common/common_types.jl")
include("completion_hints/completion_hints_types.jl")
include("serialized_result_set/serialized_result_set_types.jl")
include("extension_functions/extension_functions_types.jl")
include("mapd/mapd_types.jl")
include("mapd/MapD.jl")

#Package code
include("enums.jl") #hand-maintained code for more convenient enums
include("constructors.jl")
include("ipc.jl")
include("client.jl")
include("dataframe_show.jl")

end # module OmniSci


text:
using Test
using LinearAlgebra
using SpecialFunctions
using VoronoiDelaunay

using RootsAndPoles
import RootsAndPoles: IndexablePoint2D

const RP = RootsAndPoles

function approxmatch(A::AbstractArray, B::AbstractArray)
    length(A) == length(B) || return false

    # This could be less than O(n²) but it's not worth the bookkeeping effort
    @inbounds for i in eachindex(A)
        amatch = false
        @inbounds for j in eachindex(B)
            if A[i] ≈ B[j]
                amatch = true
                break
            end
        end
        if !amatch
            return false
        end
    end
    return true
end

@testset "RootsAndPoles" begin
    include("RootsAndPoles.jl")

    @testset "Simple Rational Function" begin include("SimpleRationalFunction.jl") end
    @testset "Complex Modes" begin include("ComplexModes.jl") end
    @testset "Lossy Multilayered Waveguide" begin include("LossyMultilayeredWaveguide.jl") end
    @testset "Graphene Transmission Line" begin include("GrapheneTransmissionLine.jl") end
    @testset "Default" begin include("DefaultFunction.jl") end
end


text:
#
# Date created: 2022-03-08
# Author: aradclif
#
#
############################################################################################
#### Experiments with programmatic generation

prog = "for i = 1:5
i + 1
end"
ex = Meta.parse(prog)

prog2 = "for i = 1:5, j = 1:4
i + j + 1
end"
ex2 = Meta.parse(prog2)

prog3 = "for i = 1:5 i + 1 end"
ex3 = Meta.parse(prog3)

prog4 = "for i = 1:5, j = 1:4 i + j + 1 end"
ex4 = Meta.parse(prog4)

Meta.show_sexpr(ex) == Meta.show_sexpr(ex3)
Meta.show_sexpr(ex2) == Meta.show_sexpr(ex4)

A = [1 2 3; 4 5 6]
B = zeros(Int, 1, 3)

prog5 = "for i = axes(A, 1)
i + 1
end"
ex5 = Meta.parse(prog5)

prog6 = "for i = axes(A, 1), j = axes(A, 2)
i + j + 1
end"
ex6 = Meta.parse(prog6)

prog7 = "for i = eachindex(A)
i + 1
end"
ex7 = Meta.parse(prog7)

#### Literal reconstructions
# building ex programmatically
_ex = Expr(:for)
e1 = Expr(:(=), :i)
e12 = Expr(:call, :(:), 1, 5)
push!(e1.args, e12)
push!(_ex.args, e1)
e2 = Expr(:block)
# e21 = :(#= none:2 =#)
# push!(e2.args, e21)
e22 = Expr(:call, :+, :i, 1)
push!(e2.args, e22)
push!(_ex.args, e2)

# building ex2 programmatically
_ex = Expr(:for)
e1 = Expr(:block)
e11 = Expr(:(=), :i)
e112 = Expr(:call, :(:), 1, 5)
push!(e11.args, e112)
push!(e1.args, e11)
e12 = Expr(:(=), :j)
e122 = Expr(:call, :(:), 1, 4)
push!(e12.args, e122)
push!(e1.args, e12)
push!(_ex.args, e1)
e2 = Expr(:block)
# e21 = :(#= none:2 =#)
# push!(e2.args, e21)
e22 = Expr(:call, :+, :i, :j, 1)
push!(e2.args, e22)
push!(_ex.args, e2)

# building ex, revised to conform with @nloops
loops = Expr(:for)
ex = Expr(:(=), Symbol(:i_, 1))
ax = Expr(:call, :axes, :A, 1)
push!(ex.args, ax)
push!(loops.args, ex)
body = Expr(:block)
e = Expr(:call, :+, Symbol(:i_, 1), 1)
# e = Expr(:call, :println, Expr(:call, :+, Symbol(:i_, 1), 1))
# e = Expr(:macrocall, :@show, Expr(:call, :+, Symbol(:i_, 1), 1))
push!(body.args, e)
push!(loops.args, body)

# building ex2, revised to conform with @nloops
N = ndims(A)
loops = Expr(:for)
block = Expr(:block)
for d = N:-1:1
    ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
    push!(block.args, ex)
end
push!(loops.args, block)
body = Expr(:block)
# create e somewhere...
# a couple examples:
# e = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
# e = Expr(:call, :+, Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...), 1)
a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
e = Expr(:(=), a, Expr(:call, :+, a, 1))
push!(body.args, e)
push!(loops.args, body)

function loopgen(N::Int) #A::AbstractArray{T, N} where {T, N}
    if N == 1
        loops = Expr(:for)
        ex = Expr(:(=), Symbol(:i_, 1), Expr(:call, :axes, :A, 1))
        push!(loops.args, ex)
    else
        loops = Expr(:for)
        block = Expr(:block)
        for d = N:-1:1
            ex = Expr(:(=), Symbol(:i_, d), Expr(:call, :axes, :A, d))
            push!(block.args, ex)
        end
        push!(loops.args, block)
    end
    return loops
end
loopgen(A::AbstractArray{T, N}) where {T, N} = loopgen(N)

ls = loopgen(A)
push!(ls.args, body)

# Another experiment
h(x, y) = x * x + y
body = Expr(:block)
a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
e = Expr(:(=), a, Expr(:call, :h, a, 1))
push!(body.args, e)
ls = loopgen(A)
push!(ls.args, body)

function reducebody1(f, A::AbstractArray{T, N}) where {T, N}
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    e = Expr(:(=), a, Expr(:call, Symbol(f), a, 1))
    push!(body.args, e)
    body
end

ls = loopgen(A)
body = reducebody1(h, A)
push!(ls.args, body)

g(x, y) = x * (x + y)
ls = loopgen(A)
body = reducebody1(g, A)
push!(ls.args, body)

# Body expressions
ls = loopgen(A)
body = Expr(:block)
a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
b = Expr(:ref, :B, ntuple(d -> length(axes(B, d)) == 1 ? 1 : Symbol(:i_, d), N)...)
e = Expr(:(=), b, Expr(:call, :+, b, a))
# e = Expr(:(=), b, Expr(:call, :h, b, a))
push!(body.args, e)
push!(ls.args, body)

function reducebody(f, B::AbstractArray{T, N}, A::AbstractArray{T, N}) where {T, N}
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> length(axes(B, d)) == 1 ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    body
end

ls = loopgen(A)
body = reducebody(+, B, A)
push!(ls.args, body)
# Anonymous
y = (x, y) -> x * x
ls = loopgen(A)
body = reducebody(y, B, A)
push!(ls.args, body)

function reduce_quote(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = similar(A, Base.promote_op(f, T), Dᴮ)
    ls = loopgen(A)
    body = reducebody(f, B, A)
    push!(ls.args, body)
    ls
end

function reducebody(f, N, D)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    # type-stable? slower, actually
    # a = Expr(:ref, :A)
    # b = Expr(:ref, :B)
    # append!(b.args, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N))
        # for d = 1:N
    #     push!(a.args, Symbol(:i_, d))
    #     push!(b.args, D[d] == 1 ? 1 : Symbol(:i_, d))
    # end
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    body
end
function mapreducebody(f, op, N, D)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> D[d] == 1 ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(op), b, Expr(:call, Symbol(f), a)))
    push!(body.args, e)
    body
end
reducebody(|, 2, (1,3))
mapreducebody(abs2, +, 2, (1,3))
#### An attempt at reduce quote for use in @generated
function reduce_quote(F, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    push!(ls.args, body)
    # return quote
    #     $ls
    #     return B
    # end
    ls
end
function mapreduce_quote(F, OP, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    op = OP.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(op), b, Expr(:call, Symbol(f), a)))
    push!(body.args, e)
    push!(ls.args, body)
    # return quote
    #     $ls
    #     return B
    # end
    ls
end
Db = size(B)
D_b = ntuple(d -> StaticInt(Db[d]), N)
reduce_quote(typeof(+), 2, typeof(D_b))
m2reduce_quote(typeof(+), 2, typeof(D_b))

function mreduce(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = zeros(Base.promote_op(f, T), Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _mreduce!(f, B, A, Dᴮ′)
    B
end

@generated function _mreduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, T, N, D}
    reduce_quote(F, N, D)
end
function mapmreduce(f, op, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = zeros(Base.promote_op(f, T), Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _mapmreduce!(f, op, B, A, Dᴮ′)
    B
end

@generated function _mapmreduce!(f::F, op::OP, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, OP, T, N, D}
    mapreduce_quote(F, OP, N, D)
end

A = rand(1000, 1000, 1000);
@benchmark mreduce(+, A, (1,))
@benchmark reduce(+, A, dims=(3,))
@benchmark m2reduce(+, A, (3,))
@benchmark m2treduce(+, A, (3,))
@benchmark mapmreduce(abs2, +, A, (3,))
@benchmark mapreduce(abs2, +, A, dims=(3,))
@benchmark mapm2reduce(abs2, +, A, (3,))
@benchmark mapm2treduce(abs2, +, A, (3,))
mreduce(|, A, (1,))
reduce(|, A, dims=(1,), init=0)
m2reduce(+, A, (1,))
m2treduce(+, A, (1,))
mapmreduce(abs2, +, A, (1,))
mapreduce(abs2, +, A, dims=(1,))
mapm2reduce(abs2, +, A, (1,))
mapm2treduce(abs2, +, A, (1,))

function m2reduce_quote(F, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return B
    end
end
function m2reduce(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = zeros(Base.promote_op(f, T), Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _m2reduce!(f, B, A, Dᴮ′)
    B
end
@generated function _m2reduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, T, N, D}
    m2reduce_quote(F, N, D)
end

function m2treduce_quote(F, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(f), b, a))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return B
    end
end
function m2treduce(f, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = zeros(Base.promote_op(f, T), Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _m2treduce!(f, B, A, Dᴮ′)
    B
end
@generated function _m2treduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, T, N, D}
    m2treduce_quote(F, N, D)
end

function mapm2reduce_quote(F, OP, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    op = OP.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(op), b, Expr(:call, Symbol(f), a)))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return B
    end
end
function mapm2reduce(f, op, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = zeros(Base.promote_op(f, T), Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _mapm2reduce!(f, op, B, A, Dᴮ′)
    B
end

@generated function _mapm2reduce!(f::F, op::OP, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, OP, T, N, D}
    mapm2reduce_quote(F, OP, N, D)
end

function mapm2treduce_quote(F, OP, N::Int, D)
    ls = loopgen(N)
    body = Expr(:block)
    params = D.parameters
    f = F.instance
    op = OP.instance
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), b, Expr(:call, Symbol(op), b, Expr(:call, Symbol(f), a)))
    push!(body.args, e)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return B
    end
end
function mapm2treduce(f, op, A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    Dᴬ = size(A)
    Dᴮ = ntuple(d -> d ∈ dims ? 1 : Dᴬ[d], N)
    B = zeros(Base.promote_op(f, T), Dᴮ)
    Dᴮ′ = ntuple(d -> StaticInt(Dᴮ[d]), N)
    _mapm2treduce!(f, op, B, A, Dᴮ′)
    B
end

@generated function _mapm2treduce!(f::F, op::OP, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, OP, T, N, D}
    mapm2treduce_quote(F, OP, N, D)
end

################################################################
# 2022-03-09
# Tests
n = 4
A = rand(rand(1:100, n)...);
for pre ∈ (:lv, :lvt), post ∈ (:sum, :prod, :maximum, :minimum)
    φ = Symbol(pre, post)
    @timev @eval $φ(A, dims=(2,4))
    A = rand(rand(1:100, n)...);
    @timev @eval $φ(A, dims=(2,4))
    @assert @eval $φ(A, dims=(2,4)) ≈ $post(A, dims=(2,4))
end
################################################################


text:
using Base.Test
using Meganet

npixel   = 20
nex      = 12
nchannel = 3

for TYPE=[Float64,Float32]
    L        = getAffineScalingLayer(TYPE,[npixel,nchannel,nex])

    @testset "affineScalingLayer $TYPE" begin
        testAbstractMeganetElement(L)
    end
end