text:
using EscapeAnalysis, InteractiveUtils, Test, JET

mutable struct MutableSome{T}
    value::T
end
mutable struct MutableCondition
    cond::Bool
end

@testset "EscapeAnalysis" begin

@testset "basics" begin
    let # simplest
        result = analyze_escapes((Any,)) do a # return to caller
            return nothing
        end
        @test has_return_escape(result.state.arguments[2])
    end

    let # return
        result = analyze_escapes((Any,)) do a
            return a
        end
        @test has_return_escape(result.state.arguments[1]) # self
        @test has_return_escape(result.state.arguments[2]) # argument
    end

    let # global store
        result = analyze_escapes((Any,)) do a
            global aa = a
            return nothing
        end
        @test has_global_escape(result.state.arguments[2])
    end

    let # global load
        result = analyze_escapes() do
            global gr
            return gr
        end
        i = findfirst(has_return_escape, result.state.ssavalues)
        @assert !isnothing(i)
        has_global_escape(result.state.ssavalues[i])

        result = analyze_escapes((Any,)) do a
            global gr
            (gr::MutableSome{Any}).x = a
        end
        @test has_global_escape(result.state.arguments[2])
    end

    # https://github.com/aviatesk/EscapeAnalysis.jl/pull/16
    let # don't propagate escape information for bitypes
        result = analyze_escapes((Nothing,)) do a
            global bb = a
        end
        @test !(has_global_escape(result.state.arguments[2]))
    end
end

@testset "control flows" begin
    let # branching
        result = analyze_escapes((Any,Bool,)) do a, c
            if c
                return nothing # a doesn't escape in this branch
            else
                return a # a escapes to a caller
            end
        end
        @test has_return_escape(result.state.arguments[2])
    end

    let # π node
        result = analyze_escapes((Any,)) do a
            if isa(a, Regex)
                identity(a) # compiler will introduce π node here
                return a    # return escape !
            else
                return nothing
            end
        end
        @assert any(@nospecialize(x)->isa(x, Core.PiNode), result.ir.stmts.inst)
        @test has_return_escape(result.state.arguments[2])
    end

    let # loop
        result = analyze_escapes((Int,)) do n
            c = MutableCondition(false)
            while n > 0
                rand(Bool) && return c
            end
            nothing
        end
        i = findfirst(==(MutableCondition), result.ir.stmts.type)
        @assert !isnothing(i)
        @test has_return_escape(result.state.ssavalues[i])
    end

    let # exception
        result = analyze_escapes((Any,)) do a
            try
                nothing
            catch err
                return a # return escape
            end
        end
        @test has_return_escape(result.state.arguments[2])
    end
end

let # more complex
    result = analyze_escapes((Bool,)) do c
        x = Vector{MutableCondition}() # return escape
        y = MutableCondition(c) # return escape
        if c
            push!(x, y)
            return nothing
        else
            return x # return escape
        end
    end

    i = findfirst(==(Vector{MutableCondition}), result.ir.stmts.type)
    @assert !isnothing(i)
    @test has_return_escape(result.state.ssavalues[i])
    i = findfirst(==(MutableCondition), result.ir.stmts.type)
    @assert !isnothing(i)
    @test has_return_escape(result.state.ssavalues[i])
end

let # simple allocation
    result = analyze_escapes((Bool,)) do c
        mm = MutableCondition(c) # just allocated, never escapes
        return mm.cond ? nothing : 1
    end

    i = findfirst(==(MutableCondition), result.ir.stmts.type) # allocation statement
    @assert !isnothing(i)
    @test has_no_escape(result.state.ssavalues[i])
end

@testset "inter-procedural" begin
    m = Module()

    # FIXME currently we can't prove the effect-freeness of `getfield(RefValue{String}, :x)`
    # because of this check https://github.com/JuliaLang/julia/blob/94b9d66b10e8e3ebdb268e4be5f7e1f43079ad4e/base/compiler/tfuncs.jl#L745
    # and thus it leads to the following two broken tests

    @eval m @noinline f_no_escape(x) = (broadcast(identity, x); nothing)
    let
        result = @eval m $analyze_escapes() do
            f_no_escape(Ref("Hi"))
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test_broken has_no_escape(result.state.ssavalues[i])
    end

    @eval m @noinline f_no_escape2(x) = broadcast(identity, x)
    let
        result = @eval m $analyze_escapes() do
            f_no_escape2(Ref("Hi"))
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test_broken has_no_escape(result.state.ssavalues[i])
    end

    @eval m @noinline f_global_escape(x) = (global xx = x) # obvious escape
    let
        result = @eval m $analyze_escapes() do
            f_global_escape(Ref("Hi"))
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_global_escape(result.state.ssavalues[i])
    end

    # if we can't determine the matching method statically, we should be conservative
    let
        result = @eval m $analyze_escapes((Ref{Any},)) do a
            may_exist(a)
        end
        @test has_all_escape(result.state.arguments[2])
    end
    let
        result = @eval m $analyze_escapes((Ref{Any},)) do a
            Base.@invokelatest f_no_escape(a)
        end
        @test has_all_escape(result.state.arguments[2])
    end

    # handling of simple union-split (just exploit the inliner's effort)
    @eval m begin
        @noinline unionsplit_noescape(_)      = string(nothing)
        @noinline unionsplit_noescape(a::Int) = a + 10
    end
    let
        T = Union{Int,Nothing}
        result = @eval m $analyze_escapes(($T,)) do value
            a = $MutableSome{$T}(value)
            unionsplit_noescape(a.value)
            return nothing
        end
        inds = findall(==(MutableSome{T}), result.ir.stmts.type) # find allocation statement
        @assert !isempty(inds)
        for i in inds
            @test has_no_escape(result.state.ssavalues[i])
        end
    end

    # appropriate conversion of inter-procedural context
    # https://github.com/aviatesk/EscapeAnalysis.jl/issues/7
    @eval m @noinline f_no_escape_simple(a) = Base.inferencebarrier(nothing)
    let
        result = @eval m $analyze_escapes() do
            aaa = Ref("foo") # shouldn't be "return escape"
            a = f_no_escape_simple(aaa)
            nothing
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_no_escape(result.state.ssavalues[i])
    end
    let
        result = @eval m $analyze_escapes() do
            aaa = Ref("foo") # still should be "return escape"
            a = f_no_escape_simple(aaa)
            return aaa
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_return_escape(result.state.ssavalues[i])
    end

    # should propagate escape information imposed on return value to the aliased call argument
    @eval m @noinline function f_return_escape(a)
        println("hi") # prevent inlining
        return a
    end
    let
        result = @eval m $analyze_escapes() do
            obj = Ref("foo")           # should be "return escape"
            ret = f_return_escape(obj)
            return ret                 # alias of `obj`
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_return_escape(result.state.ssavalues[i])
    end

    @eval m @noinline function f_no_return_escape(a)
        println("hi") # prevent inlining
        return "hi"
    end
    let
        result = @eval m $analyze_escapes() do
            obj = Ref("foo")              # better to not be "return escape"
            ret = f_no_return_escape(obj)
            return ret                    # must not alias to `obj`
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test !has_return_escape(result.state.ssavalues[i])
    end
end

@testset "builtins" begin
    let # throw
        r = analyze_escapes((Any,)) do a
            throw(a)
        end
        @test has_thrown_escape(r.state.arguments[2])
    end

    let # implicit throws
        r = analyze_escapes((Any,)) do a
            getfield(a, :may_not_field)
        end
        @test has_thrown_escape(r.state.arguments[2])

        r = analyze_escapes((Any,)) do a
            sizeof(a)
        end
        @test has_thrown_escape(r.state.arguments[2])
    end

    let # :===
        result = analyze_escapes((Bool, String)) do cond, s
            m = cond ? MutableSome(s) : nothing
            c = m === nothing
            return c
        end
        i = findfirst(==(MutableSome{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_no_escape(result.state.ssavalues[i])
    end

    let # sizeof
        ary = [0,1,2]
        result = @eval analyze_escapes() do
            ary = $(QuoteNode(ary))
            sizeof(ary)
        end
        i = findfirst(==(Core.Const(ary)), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_no_escape(result.state.ssavalues[i])
    end

    let # ifelse
        result = analyze_escapes((Bool,)) do c
            r = ifelse(c, Ref("yes"), Ref("no"))
            return r
        end
        inds = findall(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isempty(inds)
        for i in inds
            @test has_return_escape(result.state.ssavalues[i])
        end
    end
    let # ifelse (with constant condition)
        result = analyze_escapes() do
            r = ifelse(true, Ref("yes"), Ref(nothing))
            return r
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_return_escape(result.state.ssavalues[i])
        i = findfirst(==(Base.RefValue{Nothing}), result.ir.stmts.type) # find allocation statement
        @assert !isnothing(i)
        @test has_no_escape(result.state.ssavalues[i])
    end

    let # typeassert
        result = analyze_escapes((Any,)) do x
            y = x::String
            return y
        end
        r = findfirst(x->isa(x,Core.ReturnNode), result.ir.stmts.inst)
        @assert !isnothing(r)
        @test has_return_escape(result.state.arguments[2], r)
        @test !has_all_escape(result.state.arguments[2])
    end
end

@testset "Exprs" begin
    let
        result = analyze_escapes((String,)) do s
            m = MutableSome(s)
            GC.@preserve m begin
                return nothing
            end
        end
        i = findfirst(==(MutableSome{String}), result.ir.stmts.type) # find allocation statement
        @test !isnothing(i)
        @test has_no_escape(result.state.ssavalues[i])
    end

    let # :isdefined
        result = analyze_escapes((String, Bool, )) do a, b
            if b
                s = Ref(a)
            end
            return @isdefined(s)
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # find allocation statement
        @test !isnothing(i)
        @test has_no_escape(result.state.ssavalues[i])
    end

    let # :foreigncall
        result = analyze_escapes((Vector{String}, Int, )) do a, b
            return isassigned(a, b) # TODO: specialize isassigned
        end
        @test has_all_escape(result.state.arguments[2])
    end
end

@testset "field analysis" begin
    let
        result = analyze_escapes((String,)) do a # => ReturnEscape
            o = MutableSome(a) # no need to escape
            f = getfield(o, :value)
            return f
        end
        i = findfirst(==(MutableSome{String}), result.ir.stmts.type)
        r = findfirst(x->isa(x, Core.ReturnNode), result.ir.stmts.inst)
        @assert !isnothing(i) && !isnothing(r)
        @test has_return_escape(result.state.arguments[2], r)
        @test_broken !has_return_escape(result.state.ssavalues[i], r)
    end

    let
        result = analyze_escapes((String,)) do a # => ReturnEscape && GlobalEscape
            o = MutableSome(a) # no need to escape
            global oo = o
            f = getfield(o, :value)
            return f
        end
        i = findfirst(==(MutableSome{String}), result.ir.stmts.type)
        r = findfirst(x->isa(x, Core.ReturnNode), result.ir.stmts.inst)
        @assert !isnothing(i) && !isnothing(r)
        @test has_return_escape(result.state.arguments[2], r)
        @test has_global_escape(result.state.ssavalues[i])
        @test_broken !has_return_escape(result.state.ssavalues[i], r)
    end

    let
        result = analyze_escapes((String,)) do a # => ReturnEscape
            o1 = MutableSome(a) # => ReturnEscape
            o2 = MutableSome(o1) # no need to escape
            return getfield(o2, :value)
        end
        i1 = findfirst(==(MutableSome{String}), result.ir.stmts.type)
        i2 = findfirst(==(MutableSome{MutableSome{String}}), result.ir.stmts.type)
        r = findfirst(x->isa(x, Core.ReturnNode), result.ir.stmts.inst)
        @assert !isnothing(i1) && !isnothing(i2) && !isnothing(4)
        @test has_return_escape(result.state.arguments[2], r)
        @test has_return_escape(result.state.ssavalues[i1])
        @test_broken !has_return_escape(result.state.ssavalues[i2])
    end

    let
        result = analyze_escapes((Any,)) do a # => GlobalEscape
            t = tuple(a) # no need to escape
            global tt = t[1]
            return nothing
        end
        i = findfirst(t->t<:Tuple, result.ir.stmts.type) # allocation statement
        @assert !isnothing(i)
        @test has_global_escape(result.state.arguments[2])
        @test_broken !has_global_escape(result.state.ssavalues[i])
    end
end

# demonstrate a simple type level analysis can sometimes compensate the lack of yet unimplemented analysis
@testset "special-casing bitstype" begin
    let # lack of field analysis
        result = analyze_escapes((Int,)) do a
            o = MutableSome(a) # no need to escape
            f = getfield(o, :value)
            return f
        end
        i = findfirst(==(MutableSome{Int}), result.ir.stmts.type) # allocation statement
        @assert !isnothing(i)
        @test has_no_escape(result.state.ssavalues[i])
    end

    let # an escaped tuple stmt will not propagate to its Int argument (since Int is of bitstype)
        result = analyze_escapes((Int, Any, )) do a, b
            t = tuple(a, b)
            global tt = t
            return nothing
        end
        @test has_return_escape(result.state.arguments[2])
        @test has_global_escape(result.state.arguments[3])
    end
end

@testset "return flow-sensitivity" begin
    isa2(T) = x->isa(x,T)

    let
        result = analyze_escapes((Bool,)) do cond
            r = Ref("foo")
            if cond
                return cond
            end
            return r
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # allocation statement
        @assert !isnothing(i)
        rts = findall(isa2(Core.ReturnNode), result.ir.stmts.inst) # return statement
        @assert length(rts) == 2
        @test count(rt->has_return_escape(result.state.ssavalues[i], rt), rts) == 1
    end

    let
        result = analyze_escapes((Bool,)) do cond
            r = Ref("foo")
            cnt = 0
            while rand(Bool)
                cnt += 1
                rand(Bool) && return r
            end
            rand(Bool) && return r
            return cnt
        end
        i = findfirst(==(Base.RefValue{String}), result.ir.stmts.type) # allocation statement
        @assert !isnothing(i)
        rts = findall(isa2(Core.ReturnNode), result.ir.stmts.inst) # return statement
        @assert length(rts) == 3
        @test count(rt->has_return_escape(result.state.ssavalues[i], rt), rts) == 2
    end
end

# TODO
mutable struct WithFinalizer
    v
    function WithFinalizer(v)
        x = new(v)
        f(t) = @async println("Finalizing $t.")
        return finalizer(x, x)
    end
end
make_m(v = 10) = MyMutable(v)
function simple(cond)
    m = make_m()
    if cond
        # println(m.v)
        return nothing # <= insert `finalize` call here
    end
    return m
end

@testset "finalizer elision" begin
    @test can_elide_finalizer(EscapeAnalysis.NoEscape(), 1)
    @test !can_elide_finalizer(EscapeAnalysis.ReturnEscape(1), 1)
    @test can_elide_finalizer(EscapeAnalysis.ReturnEscape(1), 2)
    @test !can_elide_finalizer(EscapeAnalysis.ArgumentReturnEscape(), 1)
    @test !can_elide_finalizer(EscapeAnalysis.GlobalEscape(), 1)
    @test can_elide_finalizer(EscapeAnalysis.ThrownEscape(), 1)
end

@testset "code quality" begin
    # assert that our main routine are free from (unnecessary) runtime dispatches

    function function_filter(@nospecialize(ft))
        ft === typeof(Core.Compiler.widenconst) && return false # `widenconst` is very untyped, ignore
        ft === typeof(EscapeAnalysis.escape_builtin!) && return false # `escape_builtin!` is very untyped, ignore
        ft === typeof(isbitstype) && return false # `isbitstype` is very untyped, ignore
        return true
    end

    test_opt(only(methods(EscapeAnalysis.find_escapes)).sig; function_filter)
    for m in methods(EscapeAnalysis.escape_builtin!)
        Base._methods_by_ftype(m.sig, 1, Base.get_world_counter()) === false && continue
        test_opt(m.sig; function_filter)
    end
end

end # @testset "EscapeAnalysis" begin


text:
export  PhyHermite, degree,
        FamilyPhyHermite, FamilyScaledPhyHermite,
        DPhyPolyHermite,
        FamilyDPhyPolyHermite, FamilyDScaledPhyPolyHermite,
        FamilyD2PhyPolyHermite, FamilyD2ScaledPhyPolyHermite,
        derivative!, derivative,
        evaluate!, evaluate,
        vander!, vander

"""
PhyHermite <: Hermite

An immutable structure for physicist Hermite functions defined as ψm(x) = Hm(x)*exp(-x^2/2).

## Fields
-  `m` : order of the function
-  `Poly` : physicist Hermite polynomial of order m
- `scaled` : with rescaling to have unitary norm

## Constructors
- `PhyHermite(m, Poly, scaled)`
- `PhyHermite(m; scaled = false)``
"""
struct PhyHermite <: Hermite
    m::Int64
    Poly::PhyPolyHermite
    scaled::Bool
end

function Base.show(io::IO, P::PhyHermite)
    println(io, string(P.m)*"-th order physicist Hermite function, scaled = "*string(P.scaled))
end

PhyHermite(m::Int64; scaled::Bool = false) = PhyHermite(m, PhyPolyHermite(m; scaled = scaled), scaled)

"""
$(TYPEDSIGNATURES)

Return the degree of the physicist Hermite function `P`.
"""
degree(P::PhyHermite) = P.m

"""
$(TYPEDSIGNATURES)

Evaluates the physicist Hermite function `P`.
"""
(P::PhyHermite)(x) = P.Poly.P(x)*exp(-x^2/2)

const FamilyPhyHermite = ntuple(i->PhyHermite(i-1),CstMaxDegree+1)
const FamilyScaledPhyHermite = ntuple(i->PhyHermite(i-1; scaled = true),CstMaxDegree+1)

# Store P′n - Pn * X with Pn the n-th Physicist Hermite Polynomial
"""
$(TYPEDSIGNATURES)

Store P'ₙ - Pₙ * X with Pₙ the n-th Physicist Hermite Polynomial
"""
function DPhyPolyHermite(m::Int64; scaled::Bool)
    k = 1
    factor = 2^k*exp(loggamma(m+1) - loggamma(m+1-k))
    P = PhyPolyHermite(m; scaled = false)
    Pprime = derivative(P, 1)
    if scaled == false
        return Pprime - ImmutablePolynomial((0.0, 1.0))*P.P
    else
        C = 1/Cphy(m)
        return C*(Pprime - ImmutablePolynomial((0.0, 1.0))*P.P)
    end
end

const FamilyDPhyPolyHermite = ntuple(i->AbstractPhyHermite(DPhyPolyHermite(i-1; scaled = false), false), CstMaxDegree+1)
const FamilyDScaledPhyPolyHermite = ntuple(i->AbstractPhyHermite(DPhyPolyHermite(i-1; scaled = true), true), CstMaxDegree+1)

# Store P″n - 2 P′n * X  + Pn * (X^2 - 1) with Pn the n-th Physicist Hermite Polynomial
"""
$(TYPEDSIGNATURES)

Store P″ₙ  - 2 P′ₙ * X  + Pₙ * (X^2 - 1) with Pₙ the n-th Physicist Hermite Polynomial
"""
function D2PhyPolyHermite(m::Int64; scaled::Bool)
    k = 1
    factor = 2^k*exp(loggamma(m+1) - loggamma(m+1-k))
    P = PhyPolyHermite(m; scaled = false)
    Pprime  = derivative(P, 1)
    Ppprime = derivative(P, 2)
    if scaled == false
        return Ppprime - ImmutablePolynomial((0.0, 2.0))*Pprime +
               P.P*ImmutablePolynomial((-1.0, 0.0, 1.0))
    else
        C = 1/Cphy(m)
        return C*(Ppprime - ImmutablePolynomial((0.0, 2.0))*Pprime +
                  P.P*ImmutablePolynomial((-1.0, 0.0, 1.0)))
    end
end

const FamilyD2PhyPolyHermite = ntuple(i->AbstractPhyHermite(D2PhyPolyHermite(i-1; scaled = false), false), CstMaxDegree+1)
const FamilyD2ScaledPhyPolyHermite = ntuple(i->AbstractPhyHermite(D2PhyPolyHermite(i-1; scaled = true), true), CstMaxDegree+1)

"""
$(TYPEDSIGNATURES)

Computes in-place the k-th derivative of the physicist Hermite function `F`.
"""
function derivative!(dF, F::PhyHermite, k::Int64, x)
    m = F.m
    @assert k>-2 "anti-derivative is not implemented for k<-1"
    @assert size(dF,1) == size(x,1) "Size of dF and x don't match"
    N = size(x,1)
    if k==0
         # map!(xi -> F.Poly.P(xi)*exp(-xi^2/2), dF, x)
         @avx @. dF = F.Poly.P(x)*exp(-x^2/2)
         return  dF
    elseif k==1
        if F.scaled ==false
            Pprime = FamilyDPhyPolyHermite[m+1]
            @avx @. dF = Pprime.P(x)*exp(-x^2/2)
            return dF
        else
            Pprime = FamilyDScaledPhyPolyHermite[m+1]
            @avx @. dF = Pprime.P(x)*exp(-x^2/2)
            return dF
        end
    elseif k==2
        if F.scaled ==false
            Ppprime = FamilyD2PhyPolyHermite[m+1]
            @avx @. dF = Ppprime.P(x)*exp(-x^2/2)
            return dF
        else
            Ppprime = FamilyD2ScaledPhyPolyHermite[m+1]

            @avx @. dF = Ppprime.P(x)*exp(-x^2/2)
            return dF
        end
        # map!(xi->ForwardDiff.derivative(y->ForwardDiff.derivative(z->F(z), y), xi), dF, x)
    elseif k==-1 # Compute the anti-derivatives
        # extract monomial coefficients
        c = F.Poly.P.coeffs
        # Compute integral using monomial formula
        Cst = 0.0
        pi_erf_c = 0.0
        exp_c = zeros(m)

        for nn = 0:m
            if mod(nn,2)==0
                pi_erf_c += (c[nn+1]*fact2(nn-1))
                for k = 1:ceil(Int64, nn/2)
                    exp_c[nn-2*k+1+1] += (c[nn+1]*fact2(nn-1) / fact2(nn-2*k+1))
                end
            else
                Cst += c[nn+1] * fact2(nn-1)
                for k = 0:ceil(Int64, (nn-1)/2)
                    exp_c[nn-2*k-1+1] += (c[nn+1]*fact2(nn-1) / fact2(nn-2*k-1))
                end

            end
        end
        # Evaluate exponential terms
        ev_pi_erf = zeros(N)
        @. ev_pi_erf = sqrt(π/2)*erf(1/sqrt(2)*x)
        ev_exp = zeros(N)
        @. ev_exp = exp(-x^2/2)

        # Evaluate x^N matrix
        ntot = size(exp_c, 1)
        X = ones(N, ntot)
        @inbounds for i = 2:ntot
            X[:,i] = X[:,i-1] .*x
        end
        # Compute integral
        dF .= Cst .+ pi_erf_c .* ev_pi_erf - (X*exp_c) .* ev_exp
        return dF
    end
end

"""
$(TYPEDSIGNATURES)

Computes the k-th derivative of the physicist Hermite function `F`.
"""
derivative(F::PhyHermite, k::Int64, x::Array{Float64,1}) = derivative!(zero(x), F, k, x)

"""
$(TYPEDSIGNATURES)

Evaluates in-place the physicist Hermite function `P`.
"""
function evaluate!(dV, P::PhyHermite, x)
    m = P.m
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    # H_0(x) = 1, H_1(x) = 2*x
    col0 = view(dV,:,1)
    @avx @. col0 = exp(-0.5*x^2)

    if P.scaled
        rmul!(col0, 1/Cphy(0))
    end
    if m == 0
        return dV
    end

    if P.scaled
        rmul!(col0, Cphy(0))
    end

    col = view(dV,:,2)
    @avx @. col = 2.0*x*col0

    if P.scaled
        rmul!(col0, 1/Cphy(0))
        rmul!(col , 1/Cphy(1))
    end

    if m == 1
        return dV
    end

    if P.scaled
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = 2.0*(Cphy(i-1)* x * col - Cphy(i-2)*(i-1)*colm1)
            # dV[:,i+1] = 2.0*Cphy(i-1)*x .* dV[:,i] - 2.0*Cphy(i-2)*i*dV[:,i-1]
            rmul!(colp1, 1/Cphy(i))
        end
    else
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = 2.0*(x * col - (i-1)*colm1)
        end
    end
    return dV
end

"""
$(TYPEDSIGNATURES)

Evaluates the physicist Hermite function `P`.
"""
evaluate(P::PhyHermite, x::Array{Float64,1}) = evaluate!(zeros(size(x,1), P.m+1), P, x)

"""
$(TYPEDSIGNATURES)

Computes in-place the Vandermonde matrix of the k-th derivative of the physicist Hermite function `P`.
"""
function vander!(dV, P::PhyHermite, k::Int64, x)
    m = P.m

    if k==0
        evaluate!(dV, P, x)
        return dV
    elseif k==1
        # ψ′_n(x) = Q_n(x) exp(-x^2/2)
        # Q_0(x) = -x
        # Q_1(x) = -2x^2 - 2
        evaluate!(dV, P, x)
        if m==0
            col0 = view(dV,:,1)
            @avx @. col0 *= -x
            return dV
        end

        col0 = view(dV,:,1)
        col1 = view(dV,:,2)

        N = size(x,1)
        ψn = zeros(N)
        ψnp1 = zeros(N)

        copy!(ψn, col0)
        copy!(ψnp1, col1)

        @avx @. col0 *= -x

        if P.scaled == true
            @avx @. col1 = 2.0*Cphy(0)/Cphy(1)*(1.0 - x^2) * ψn
        else
            @avx @. col1 = 2.0*(1.0 - x^2) * ψn
        end

        @inbounds for i=2:m
            copy!(ψn, ψnp1)
            colp1 = view(dV,:,i+1)
            copy!(ψnp1, colp1)
            if P.scaled == true
                @avx @. colp1 = (1/Cphy(i))*(2*(i)*Cphy(i-1)*ψn - x * Cphy(i)*ψnp1)
            else
                @avx @. colp1 = (2*(i)*ψn - x * ψnp1)
            end
        end

        return dV
    elseif k==2
        # Use the relation ψn″(x)  + (2*n + 1 -x^2) ψn(x) = 0
        evaluate!(dV, P, x)
        @inbounds for i=0:m
            col = view(dV,:,i+1)
            @avx @. col *= (-2.0*i - 1.0 + x^2)
        end
        return dV
    elseif k==3
        # Use the relation ψn‴(x)  -2x ψn(x) + (2*n + 1 -x^2) ψn′(x) = 0
        # and ψn′(x) = 2*n ψ_{n-1}(x) - x ψ_{n}(x)
        # So ψn‴(x) = (2n +3 -x^2)*x*ψn(x) - (2*n +1 -x^2)*2*n*ψn-1(x)
        evaluate!(dV, P, x)
        col0 = view(dV,:,1)
        ψn   = zero(x)
        ψnm1 = zero(x)
        copy!(ψn, col0)
        @avx @. col0 *= (3*x -x^3)

        if m==0
            return dV
        end

        @inbounds for i=1:m
            copy!(ψnm1, ψn)
            col = view(dV,:,i+1)
            copy!(ψn, col)
            if P.scaled ==true
                col .= ((2*i+3)*x -x .^3) .* ψn -2*i*Cphy(i-1)/Cphy(i)*(2*i+1 .- x .^2) .* ψnm1
            else
                col .= ((2*i+3)*x -x .^3) .* ψn -2*i*(2*i+1 .- x .^2) .* ψnm1
            end
        end
        return dV
    end
end

"""
$(TYPEDSIGNATURES)

Computes the Vandermonde matrix of the k-th derivative of the physicist Hermite function `P`.
"""
vander(P::PhyHermite, k::Int64, x::Array{Float64,1}) = vander!(zeros(size(x,1), P.m+1), P, k, x)

ext: jl


text:
# The lexer included within the parser
# Every public variable starts with l

# --- Def ---
mutable struct LCtx
    state_id
    states
end

mutable struct LState
    id
    # Associates a token or epsilon to a state
    # Note that before determinizing ctx, it gathers states, not ints
    transitions
    # Whether we can accept the input in this state
    terminal
end

lctx_new() = LCtx(0, [])

function lstate_new!(ctx; transitions = [], terminal = false)
    s = lstate_create(ctx.state_id += 1, transitions, terminal)
    push!(ctx.states, s)

    return s
end

lstate_create(id = 0, transitions = [], terminal = false) = LState(id, transitions, terminal)

# Links a to b by the token tok (can be empty)
function lstate_link!(a, tok, b)
    push!(a.transitions, (tok, b))
end

Base.isequal(a::LState, b::LState) = a.id == b.id
Base.hash(a::LState) = hash(a.id)

function Base.print(io::IO, s::LState)
    println(io, "LState(")
    println(io, "    id = $(s.id)")
    println(io, "    terminal = $(s.terminal)")
    for (tok, state) in s.transitions
        println(io, "    '$tok' -> $(isa(state, LState) ? state.id : state)")
    end

    print(io, ")")
end

# --- Functions ---
lend!(ctx) = lstate_new!(ctx, terminal=true)

# One token
function ltok!(ctx, tok)
    start = lstate_new!(ctx)
    stop = lstate_new!(ctx)

    lstate_link!(start, tok, stop)

    return start
end

# Links all states with an epsilon transition
function llink!(states...)
    for i in 1:length(states) - 1
        lstate_link!(states[i], "", states[i + 1])
    end
end

# !!! The initial state must be the first created (with id 1)
function determinize!(ctx)
    # Epsilon closures
    closures = [Set() for _ in 1:length(ctx.states)]

    # eps-NFA to NFA
    for (i, s) in enumerate(ctx.states)
        # Bfs
        visited = Set([s])
        q = [s]
        while length(q) > 0
            state = popfirst!(q)
            push!(closures[i], state.id)

            for (tok, next_state) in state.transitions
                # Epsilon
                if next_state ∉ visited && tok == ""
                    push!(q, next_state)
                    push!(visited, next_state)
                end
            end
        end
    end

    # NFA to DFA
    new_states = [closures[1]]
    visited_states = Set()
    states = Set([])
    init_state = nothing
    while length(new_states) != 0
        # Now an id is a set of int ids
        state_ids = pop!(new_states)
        push!(visited_states, state_ids)

        transitions = Dict()
        terminal = false

        for state_id in state_ids
            state = ctx.states[state_id]

            # This state contains a terminal state
            if state.terminal
                terminal = true
            end

            for (tok, next_state) in state.transitions
                # Ignore epsilons
                if tok == ""
                    continue
                end

                next_state_closure = closures[next_state.id]

                if !haskey(transitions, tok)
                    transitions[tok] = next_state_closure
                else
                    union!(transitions[tok], next_state_closure)
                end
            end
        end

        # Parse new states
        for (tok, next_state_ids) in transitions
            # If not parsed, add it to the stack
            if next_state_ids ∉ visited_states
                push!(new_states, next_state_ids)
            end
        end

        state = lstate_create(state_ids, transitions, terminal)

        if length(visited_states) != 1
            push!(states, state)
        else
            init_state = state
        end
    end

    # Convert states to linear and put init_state at the first position
    states = [init_state, states...]

    # Change states ids
    stateids2stateid = Dict()
    for (i, s) in enumerate(states)
        stateids2stateid[s.id] = i
    end

    for s in states
        s.id = stateids2stateid[s.id]
        for (tok, nstate) in s.transitions
            s.transitions[tok] = stateids2stateid[nstate]
        end
    end

    ctx.states = states
end

# tok2index maps a token (character) to an index, returns nothing if invalid token
function maketable(ctx, tok2index, ntoks)
    actions = Matrix{Any}(nothing, length(ctx.states), ntoks)

    for state in ctx.states
        for (tok, nstate) in state.transitions
            @assert isa(tok, Char) "Tokens must be chars (got $tok)"

            t = tok2index(tok)

            @assert t != nothing "Invalid character '$tok'"

            actions[state.id, t] = nstate
        end
    end

    terminal_states = [s.terminal for s in ctx.states]

    return ltable_new(actions, terminal_states, tok2index)
end

# --- Parsing ---
# Intermediate automaton representation
function parsereg(reg)
    states = []

    i = 1
    while i <= length(reg)
        char = reg[i]
        if char in "*?"
            @assert length(states) != 0 "Invalid $char at column $i for regex $reg"

            last_state = pop!(states)
            new_state = ("$char", last_state)
            push!(states, new_state)
        elseif char == '['
            start = i
            while i <= length(reg) && reg[i] != ']'
                i += 1
            end

            @assert i <= length(reg) "Unmatched '[' at column $start for regex $reg"

            push!(states, ("group", reg[start + 1 : i - 1]))
        else
            push!(states, ("char", reg[i]))
        end

        i += 1
    end

    return states
end

reggroups = Dict(
                 "num" => Set(String('0':'9')),
                 "alpha" => Set("_" * String('a':'z') * String('A':'Z')),
                 "\\n" => Set("\n"),
                 "\\t" => Set("\t"),
                 "\\s" => Set("\t "),
                 "\\b" => Set("\t \n"),
                 "\\:" => Set(":"),
                 "\\\"" => Set("\""),
                 "\\*" => Set("*"),
                 "\\[" => Set("["),
                 "\\\\" => Set("\\"),
                 "." => Set([Char(code) for code in 1:128]),
        )
pop!(reggroups["."], '\n')

reggroups["alnum"] = union(reggroups["num"], reggroups["alpha"])

reggroups["^'"] = Set(reggroups["."])
pop!(reggroups["^'"], '\'')

reggroups["^\\\""] = Set(reggroups["."])
pop!(reggroups["^\\\""], '"')

# Returns start automaton state and last one
# Thompson's construction basicaly
function state2nfa!(ctx, state)
    global reggroups

    (class, data) = state

    if class == "char"
        start = lstate_new!(ctx)
        stop = lstate_new!(ctx)

        lstate_link!(start, data, stop)

        return start, stop
    elseif class == "group"
        @assert haskey(reggroups, data) "No group named '$data'"

        start = lstate_new!(ctx)
        stop = lstate_new!(ctx)

        for char in reggroups[data]
            lstate_link!(start, char, stop)
        end

        return start, stop
    elseif class == "*"
        start = lstate_new!(ctx)
        mid_start, mid_stop = state2nfa!(ctx, data)
        suff = lstate_new!(ctx)
        stop = lstate_new!(ctx)

        # start ---------------------------------> stop
        #   +-> mid_start -> mid_stop -> suff +^
        #          ^----------------------+
        llink!(start, stop)
        llink!(start, mid_start)
        llink!(mid_stop, suff, stop)
        llink!(suff, mid_start)

        return start, stop
    elseif class == "?"
        start = lstate_new!(ctx)
        mid_start, mid_stop = state2nfa!(ctx, data)
        stop = lstate_new!(ctx)

        # start -------------------> stop
        #   +-> mid_start -> mid_stop ^
        llink!(start, stop)
        llink!(start, mid_start)
        llink!(mid_stop, stop)

        return start, stop
    else
        error("Invalid state '$class'")
    end
end

# Intermediate to epsilon NFA
function states2nfa(states)
    ctx = lctx_new()
    init_state = lstate_new!(ctx)
    last_state = init_state

    for state in states
        start, stop = state2nfa!(ctx, state)

        transition = lstate_new!(ctx)
        llink!(last_state, start)
        llink!(stop, transition)
        last_state = transition
    end

    terminal_state = lstate_new!(ctx, terminal=true)
    llink!(last_state, terminal_state)

    return ctx
end

# Regex to table
function compilereg(reg)
    # Compile
    states = parsereg(reg)

    # Create automaton
    ctx = states2nfa(states)
    determinize!(ctx)

    # Create table (default char to index mapping)
    table = maketable(ctx, ldefault_chartok2index, ldefault_chartok2index_length)

    return table
end

function generate_lexer(io, regs)
    parserdir = dirname(@__FILE__)

    # Compile regexes
    for (i, (name, reg, rule)) in enumerate(regs)
        regs[i] = (name, compilereg(reg), rule)
    end

    # Regs
    println(io, "regs = [")
    for (name, reg, rule) in regs
        println(io, "    ($(repr(name)), ltable_new($(reg.actions), $(reg.terminal_states), ldefault_chartok2index), $rule),")
    end
    println(io, "]")

    # Template
    println(io, read("$parserdir/lexer_template.jl", String))
end


text:
#=
Example:
- Julia version: 1.5.3
- Java version: 1.8.0_271
- Authors: Ivan Jericevich, Dharmesh Sing, Tim Gebbie
- Structure:
    1. Preliminaries
    2. Login and start session
    3. Submit orders
    4. Market data updates
    5. End session and logout
=#
#---------------------------------------------------------------------------------------------------

#----- Preliminaries -----#
# Import the Java-Julia interface package
using JavaCall
cd(@__DIR__); clearconsole() # pwd()
# Add the path to java classes as well as the path to the ".jar" files containing the required java dependencies
JavaCall.addClassPath("/home/ivanjericevich/CoinTossX/ClientSimulator/build/install/ClientSimulator/lib/*.jar") # JavaCall.getClassPath()
# Initialize JVM
JavaCall.init()
# Import the class containing the reqired methods
utilities = @jimport example.Utilities # JavaCall.listmethods(utilities)
#---------------------------------------------------------------------------------------------------

#----- Login and start session -----#
# Define the client ID corresponding the client to be logged in as well as the security ID corresponding to the security in which they will trade
clientId = 1; securityId = 1
# Load the simulation settings as well as all client data (ports, passwords and IDs) and create/initialize client
client = jcall(utilities, "loadClientData", JavaObject{Symbol("client.Client")}, (jint, jint), clientId, securityId)
# Start trading session by Logging in client to the gateways and connecting to ports
jcall(client, "sendStartMessage", Nothing, ())
#---------------------------------------------------------------------------------------------------

#----- Submit orders -----#
# Arguments for "submitOrder": order ID, volume, price, side, order type, time in force, display quantity, min execution size, stop price
jcall(client, "submitOrder", Nothing, (jlong, jlong, JString, JString, JString, jlong, jlong, jlong), "1", 1000, 99, "Buy", "Limit", "Day", 1000, 0, 0) # Buy limit order
jcall(client, "submitOrder", Nothing, (jlong, jlong, JString, JString, JString, jlong, jlong, jlong), "2", 1000, 101, "Sell", "Limit", "Day", 1000, 0, 0) # Sell limit order
jcall(client, "submitOrder", Nothing, (jlong, jlong, JString, JString, JString, jlong, jlong, jlong), "3", 1000, 0, "Buy", "Market", "Day", 1000, 0, 0) # Buy market order
jcall(client, "submitOrder", Nothing, (jlong, jlong, JString, JString, JString, jlong, jlong, jlong), "4", 1000, 0, "Buy", "StopLimit", "Day", 1000, 0, 0) # Stop buy limit order
jcall(client, "submitOrder", Nothing, (jlong, jlong, JString, JString, JString, jlong, jlong, jlong), "5", 1000, 0, "Buy", "Stop", "Day", 1000, 0, 0) # Stop buy market order
# Arguments for "cancelOrder": order id, side, price
jcall(client, "cancelOrder", Nothing, (jlong, jlong), "1", "Buy", 99) # Cancel limit order
#---------------------------------------------------------------------------------------------------

#----- Market data updates -----#
jcall(client, "calcVWAP", jlong, (JString,), "Buy") # VWAP of buy side of LOB
jcall(client, "getBid", jlong, ()) # Best bid price
jcall(client, "getBidQuantity", jlong, ()) # Best bid volume
jcall(client, "getOffer", jlong, ()) # Best ask price
jcall(client, "getOfferQuantity", jlong, ()) # Best ask volume
jcall(client, "isAuction", jboolean, ()) # Current trading session
jcall(client, "lobSnapshot", JavaObject{Symbol("java.util.ArrayList")}, ()) # Snapshot of the entire LOB
jcall(client, "getStaticPriceReference", jlong, ()) # Static price reference at the end of a trading session
#---------------------------------------------------------------------------------------------------

#----- End session and logout -----#
# End trading session by logging out client and closing connections
jcall(client, "sendEndMessage", Nothing, ()); jcall(client, "close", Nothing, ())
#---------------------------------------------------------------------------------------------------

max_stars_repo_name: specialforcea/CoinTossX

max_issues_repo_name: IvanJericevich/CoinTossX

max_forks_repo_name: IvanJericevich/CoinTossX


text:
"""
    const ghcn_root_dir

The location of the GHCN data.
"""
const ghcn_root_dir = "https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/"

"""
    const station_metadata_fname

Name of the file containing station metadata.
"""
const station_metadata_fname = "ghcnd-stations.txt"

"""
    const countries_fname

Name of the file containing country-code -- code pairing data set. See
`download_countries` and `parse_countries_metadata` for functionality.
"""
const countries_fname = "ghcnd-countries.txt"

"""
    const inventory_fname

Name of the inventory file in the GHCN daily data.
"""
const inventory_fname = "ghcnd-inventory.txt"

"""
    const ghcnd_data

Name of the file containing all of the ghcn daily data.
"""
const ghcnd_data = "ghcnd_all.tar.gz"
const ghcnd_data_dir = "ghcnd_all"

# As per the suggestion in the DataDeps doc, we perform DataDep registration in __init__.
function __init__()
    register(DataDep(
        "ghcn-stations",
        "GHCN station metadata.",
        ghcn_root_dir * station_metadata_fname,
    ))
    register(DataDep(
        "ghcn-countries",
        "GHCN country codes.",
        ghcn_root_dir * countries_fname,
    ))
    register(DataDep(
        "ghcn-inventories",
        "GHCN country codes.",
        ghcn_root_dir * inventory_fname,
    ))
    register(DataDep(
        "ghcn-data",
        "GHCN country codes.",
        ghcn_root_dir * ghcnd_data,
        post_fetch_method=unpack,
    ))
end


text:
using ClusterTrees, BEAST, CompScienceMeshes

a, h = 1.0, 0.061
#a, h = 1.0, 0.151
Γ1 = meshsphere(a, h)
Γ2 = CompScienceMeshes.translate(Γ1, point(2.10,0,0))

X1 = raviartthomas(Γ1)
X2 = raviartthomas(Γ2)

@show numfunctions(X1)
@show numfunctions(X2)

κ = 1.0; γ = κ*im;
t = Maxwell3D.singlelayer(gamma=γ)

A = ClusterTrees.h1compress(t,X1,X2)

sz = storedentries(A)
println("Compression: ",  sz/numfunctions(X1)/numfunctions(X2))

@time Z1 = assemble(t,X1,X2)
Z2 = ClusterTrees.reconstruct(A)

norm(Z1-Z2)


text:
""" Toolbelt
Location for various useful tools that might not have a better place to live
within the code base yet. Definitions listed here are likely to move at some
point.
"""
module Toolbelt
#include("mcmcsampler.jl")

end # Module Toolbelt


text:
using OpticSim
using OpticSimVis
using Test

@testset "general tests" begin
    t = OpticSim.identitytransform()
    tt = transform(t)
    @info typeof(tt)
end