text:
using SymbolicUtils: FnType, Sym, metadata
using Setfield

const IndexMap = Dict{Char,Char}(
            '-' => '₋',
            '0' => '₀',
            '1' => '₁',
            '2' => '₂',
            '3' => '₃',
            '4' => '₄',
            '5' => '₅',
            '6' => '₆',
            '7' => '₇',
            '8' => '₈',
            '9' => '₉')

struct VariableDefaultValue end
struct VariableSource end

function recurse_and_apply(f, x)
    if symtype(x) <: AbstractArray
        getindex_posthook(x) do r,x,i...
            recurse_and_apply(f, r)
        end
    else
        f(x)
    end
end

function setdefaultval(x, val)
    if symtype(x) <: AbstractArray
        if val isa AbstractArray
            getindex_posthook(x) do r,x,i...
                setdefaultval(r, val[i...])
            end
        else
            getindex_posthook(x) do r,x,i...
                setdefaultval(r, val)
            end
        end
    else
        setmetadata(x,
                    VariableDefaultValue,
                    val)
    end
end

struct GetindexParent end

function scalarize_getindex(x, parent=Ref{Any}(x))
    if symtype(x) <: AbstractArray
        parent[] = getindex_posthook(x) do r,x,i...
            scalarize_getindex(r, parent)
        end
    else
        xx = scalarize(x)
        xx = metadata(xx, metadata(x))
        if symtype(xx) <: FnType
            setmetadata(CallWithMetadata(xx, metadata(xx)), GetindexParent, parent[])
        else
            setmetadata(xx, GetindexParent, parent[])
        end
    end
end

function map_subscripts(indices)
    str = string(indices)
    join(IndexMap[c] for c in str)
end

function unwrap_runtime_var(v)
    isruntime = Meta.isexpr(v, :$) && length(v.args) == 1
    isruntime && (v = v.args[1])
    return isruntime, v
end

# Build variables more easily
function _parse_vars(macroname, type, x, transform=identity)
    ex = Expr(:block)
    var_names = Symbol[]
    # if parsing things in the form of
    # begin
    #     x
    #     y
    #     z
    # end
    x = x isa Tuple && first(x) isa Expr && first(x).head == :tuple ? first(x).args : x # tuple handling
    x = flatten_expr!(x)
    cursor = 0
    isoption(ex) = Meta.isexpr(ex, [:vect, :vcat, :hcat])
    while cursor < length(x)
        cursor += 1
        v = x[cursor]

        # We need lookahead to the next `v` to parse
        # `@variables x [connect=Flow,unit=u]`
        nv = cursor < length(x) ? x[cursor+1] : nothing
        val = unit = connect = options = nothing

        # x = 1, [connect = flow; unit = u"m^3/s"]
        if Meta.isexpr(v, :(=))
            v, val = v.args
            if Meta.isexpr(val, :tuple) && length(val.args) == 2 && isoption(val.args[2])
                options = val.args[2].args
                val = val.args[1]
            end
        end

        type′ = type

        if Meta.isexpr(v, :(::))
            v, type′ = v.args
            type′ = type′ === :Complex ? Complex{type} : type′
        end

        # x [connect = flow; unit = u"m^3/s"]
        if isoption(nv)
            options = nv.args
            cursor += 1
        end

        isruntime, v = unwrap_runtime_var(v)
        iscall = Meta.isexpr(v, :call)
        isarray = Meta.isexpr(v, :ref)
        issym  = v isa Symbol
        @assert iscall || isarray || issym "@$macroname expects a tuple of expressions or an expression of a tuple (`@$macroname x y z(t) v[1:3] w[1:2,1:4]` or `@$macroname x y z(t) v[1:3] w[1:2,1:4] k=1.0`)"

        if iscall
            isruntime, fname = unwrap_runtime_var(v.args[1])
            call_args = map(last∘unwrap_runtime_var, @view v.args[2:end])
            var_name, expr = construct_vars(macroname, fname, type′, call_args, val, options, transform, isruntime)
        else
            var_name, expr = construct_vars(macroname, v, type′, nothing, val, options, transform, isruntime)
        end

        push!(var_names, var_name)
        push!(ex.args, expr)
    end
    rhs = build_expr(:vect, var_names)
    push!(ex.args, rhs)
    return ex
end

function construct_vars(macroname, v, type, call_args, val, prop, transform, isruntime)
    issym  = v isa Symbol
    isarray = isa(v, Expr) && v.head == :ref
    if isarray
        var_name = v.args[1]
        if Meta.isexpr(var_name, :(::))
            var_name, type′ = var_name.args
            type = type′ === :Complex ? Complex{type} : type′
        end
        isruntime, var_name = unwrap_runtime_var(var_name)
        indices = v.args[2:end]
        expr = _construct_array_vars(macroname, isruntime ? var_name : Meta.quot(var_name), type, call_args, val, prop, indices...)
    else
        var_name = v
        if Meta.isexpr(v, :(::))
            var_name, type′ = v.args
            type = type′ === :Complex ? Complex{type} : type′
        end
        expr = construct_var(macroname, isruntime ? var_name : Meta.quot(var_name), type, call_args, val, prop)
    end
    lhs = isruntime ? gensym(var_name) : var_name
    rhs = :($transform($expr))
    lhs, :($lhs = $rhs)
end

function option_to_metadata_type(::Val{opt}) where {opt}
    throw(Base.Meta.ParseError("unknown property type $opt"))
end

function setprops_expr(expr, props, macroname, varname)
    expr = :($setmetadata($expr, $VariableSource, ($(Meta.quot(macroname)), $varname,)))
    isnothing(props) && return expr
    for opt in props
        if !Meta.isexpr(opt, :(=))
            throw(Base.Meta.ParseError(
                "Variable properties must be in " *
                "the form of `a = b`. Got $opt."))
        end

        lhs, rhs = opt.args

        @assert lhs isa Symbol "the lhs of an option must be a symbol"
        expr = :($setmetadata($expr,
                              $(option_to_metadata_type(Val{lhs}())),
                       $rhs))
    end
    expr
end

struct CallWithMetadata{T,M} <: Symbolic{T}
    f::Symbolic{T}
    metadata::M
end

for f in [:istree, :operation, :arguments]
    @eval SymbolicUtils.$f(x::CallWithMetadata) = $f(x.f)
end

SymbolicUtils.Code.toexpr(x::CallWithMetadata, st) = SymbolicUtils.Code.toexpr(x.f, st)

CallWithMetadata(f) = CallWithMetadata(f, nothing)

function Base.show(io::IO, c::CallWithMetadata)
    show(io, c.f)
    print(io, "⋆")
end

function (f::CallWithMetadata)(args...)
    wrap(metadata(f.f(args...), metadata(f)))
end

function construct_var(macroname, var_name, type, call_args, val, prop)
    expr = if call_args === nothing
        :($Sym{$type}($var_name))
    elseif !isempty(call_args) && call_args[end] == :..
        :($CallWithMetadata($Sym{$FnType{Tuple, $type}}($var_name)))
    else
        :($Sym{$FnType{NTuple{$(length(call_args)), Any}, $type}}($var_name)($(map(x->:($value($x)), call_args)...)))
    end

    if val !== nothing
        expr = :($setdefaultval($expr, $val))
    end

    :($wrap($(setprops_expr(expr, prop, macroname, var_name))))
end

struct CallWith
    args
end

(c::CallWith)(f) = unwrap(f(c.args...))

function _construct_array_vars(macroname, var_name, type, call_args, val, prop, indices...)
    # TODO: just use Sym here
    ndim = length(indices)

    need_scalarize = false
    expr = if call_args === nothing
        ex = :($Sym{Array{$type, $ndim}}($var_name))
        :($setmetadata($ex, $ArrayShapeCtx, ($(indices...),)))
    elseif !isempty(call_args) && call_args[end] == :..
        need_scalarize = true
        ex = :($Sym{Array{$FnType{Tuple, $type}, $ndim}}($var_name))
        ex = :($setmetadata($ex, $ArrayShapeCtx, ($(indices...),)))
        :($map($CallWithMetadata, $ex))
    else
        # [(R -> R)(R) ....]
        need_scalarize = true
        ex = :($Sym{Array{$FnType{Tuple, $type}, $ndim}}($var_name))
        ex = :($setmetadata($ex, $ArrayShapeCtx, ($(indices...),)))
        :($map($CallWith(($(call_args...),)), $ex))
    end

    if val !== nothing
        expr = :($setdefaultval($expr, $val))
    end
    expr = setprops_expr(expr, prop, macroname, var_name)
    if need_scalarize
        expr = :($scalarize_getindex($expr))
    end

    expr = :($wrap($expr))

    return expr
end

"""

Define one or more unknown variables.

```julia
@variables t α σ(..) β[1:2]
@variables w(..) x(t) y z(t, α, x)

expr = β[1]* x + y^α + σ(3) * (z - t) - β[2] * w(t - 1)
```

`(..)` signifies that the value should be left uncalled.

Symbolics supports creating variables that denote an array of some size.

```julia
julia> @variables x[1:3]
1-element Vector{Symbolics.Arr{Num, 1}}:
 x[1:3]

julia> @variables y[1:3, 1:6] # support for  tensors
1-element Vector{Symbolics.Arr{Num, 2}}:
 y[1:3,1:6]

julia> @variables t z[1:3](t) # also works for dependent variables
2-element Vector{Any}:
 t
  (map(#5, z))[1:3]
```

A symbol or expression that represents an array can be turned into an array of
symbols or expressions using the `scalarize` function.

```julia
julia> Symbolics.scalarize(z)
3-element Vector{Num}:
 z[1](t)
 z[2](t)
 z[3](t)
```

Note that `@variables` returns a vector of all the defined variables.

`@variables` can also take runtime symbol values by the `\$` interpolation
operator, and in this case, `@variables` doesn't automatically assign the value,
instead, it only returns a vector of symbolic variables. All the rest of the
syntax also applies here.

```julia
julia> a, b, c = :runtime_symbol_value, :value_b, :value_c
:runtime_symbol_value

julia> vars = @variables t \$a \$b(t) \$c[1:3](t)
4-element Vector{Any}:
      t
 runtime_symbol_value
   value_b(t)
       (map(#9, value_c))[1:3]

julia> (t, a, b, c)
(t, :runtime_symbol_value, :value_b, :value_c)
```
"""
macro variables(xs...)
    esc(_parse_vars(:variables, Real, xs))
end

TreeViews.hastreeview(x::Sym) = true

function TreeViews.treelabel(io::IO,x::Sym,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Text(x.name))
end

const _fail = Dict()

_getname(x, _) = nameof(x)
_getname(x::Symbol, _) = x
function _getname(x::Symbolic, val)
    ss = getsource(x, nothing)
    if ss === nothing
        ss = getsource(getparent(x), val)
    end
    ss === _fail && throw(ArgumentError("Variable $x doesn't have a source defined."))
    ss[2]
end

getsource(x, val=_fail) = getmetadata(unwrap(x), VariableSource, val)

getname(x, val=_fail) = _getname(unwrap(x), val)

function getparent(x, val=_fail)
    maybe_parent = getmetadata(x, Symbolics.GetindexParent, nothing)
    if maybe_parent !== nothing
        return maybe_parent
    else
        if istree(x) && operation(x) === getindex
            return arguments(x)[1]
        end
    end
    val === _fail && throw(ArgumentError("Cannot find the parent of $x."))
    return val
end

function getdefaultval(x, val=_fail)
    x = unwrap(x)
    val = getmetadata(x, VariableDefaultValue, val)
    if val !== _fail
        return val
    else
        error("$x has no default value")
    end
end

"""
    variables(name::Symbol, indices...)

Create a multi-dimensional array of individual variables named with subscript
notation. Use `@variables` instead to create symbolic array variables (as
opposed to array of variables). See `variable` to create one variable with
subscripts.

```julia-repl
julia> Symbolics.variables(:x, 1:3, 3:6)
3×4 Matrix{Num}:
 x₁ˏ₃  x₁ˏ₄  x₁ˏ₅  x₁ˏ₆
 x₂ˏ₃  x₂ˏ₄  x₂ˏ₅  x₂ˏ₆
 x₃ˏ₃  x₃ˏ₄  x₃ˏ₅  x₃ˏ₆
```
"""
function variables(name, indices...; T=Real)
    [variable(name, ij...; T=T) for ij in Iterators.product(indices...)]
end

"""
    variable(name::Symbol, idx::Integer...; T=Real)

Create a variable with the given name along with subscripted indices with the
`symtype=T`. When `T=FnType`, it creates a symbolic function.

```julia-repl
julia> Symbolics.variable(:x, 4, 2, 0)
x₄ˏ₂ˏ₀

julia> Symbolics.variable(:x, 4, 2, 0, T=Symbolics.FnType)
x₄ˏ₂ˏ₀⋆
```

Also see `variables`.
"""
function variable(name, idx...; T=Real)
    name_ij = Symbol(name, join(map_subscripts.(idx), "ˏ"))
    if T <: FnType
        first(@variables $name_ij(..))
    else
        first(@variables $name_ij::T)
    end
end

##### Renaming #####
# getname
# rename
# getindex parent
# calls
# symbolic function x[1:3](..)
#
# x_t
# sys.x

function rename_getindex_source(x, parent=x)
    getindex_posthook(x) do r,x,i...
        hasmetadata(r, GetindexParent) ? setmetadata(r, GetindexParent, parent) : r
    end
end

function rename_metadata(from, to, name)
    if hasmetadata(from, VariableSource)
        s = getmetadata(from, VariableSource)
        to = setmetadata(to, VariableSource, (s[1], name))
    end
    if hasmetadata(from, GetindexParent)
        s = getmetadata(from, GetindexParent)
        to = setmetadata(to, GetindexParent, rename(s, name))
    end
    return to
end

function rename(x::Sym, name)
    xx = @set! x.name = name
    xx = rename_metadata(x, xx, name)
    symtype(xx) <: AbstractArray ? rename_getindex_source(xx) : xx
end

rename(x::Union{Num, Arr}, name) = wrap(rename(unwrap(x), name))
function rename(x::ArrayOp, name)
    t = x.term
    args = arguments(t)
    # Hack:
    @assert operation(t) === (map) && (args[1] isa CallWith || args[1] == CallWithMetadata)
    rn = rename(args[2], name)

    xx = metadata(operation(t)(args[1], rn), metadata(x))
    rename_getindex_source(rename_metadata(x, xx, name))
end

function rename(x::CallWithMetadata, name)
    rename_metadata(x, CallWithMetadata(rename(x.f, name), x.metadata), name)
end

function rename(x::Symbolic, name)
    if istree(x) && operation(x) === getindex
        rename(arguments(x)[1], name)[arguments(x)[2:end]...]
    elseif istree(x) && symtype(operation(x)) <: FnType || operation(x) isa CallWithMetadata
        @assert x isa Term
        xx = @set x.f = rename(operation(x), name)
        @set! xx.hash = Ref{UInt}(0)
        return rename_metadata(x, xx, name)
    else
        error("can't rename $x to $name")
    end
end

# Deprecation below

struct Variable{T} end

function (::Type{Variable{T}})(s, i...) where {T}
    Base.depwarn("Variable{T}(name, idx...) is deprecated, use variable(name, idx...; T=T)", :Variable, force=true)
    variable(s, i...; T=T)
end

(::Type{Variable})(s, i...) = Variable{Real}(s, i...)

function (::Type{Sym{T}})(s, x, i...) where {T}
    Base.depwarn("Sym{T}(name, x, idx...) is deprecated, use variable(name, x, idx...; T=T)", :Variable, force=true)
    variable(s, x, i...; T=T)
end
(::Type{Sym})(s, x, i...) = Sym{Real}(s, x, i...)


text:
using Espresso: simplify, matchex, rewrite, set_default_placeholders

using Match: @match

set_default_placeholders(Set([:x, :y, :a, :b, :n]))

integrate(ex::Value, x::Symbolic = :x) = simplify ∘ _integrate ∘ simplify

function _integrate(n, x)
    @match n begin
        0 => 0
        1 => x
        n::Number => Expr(:call, :*, n, x)
        _ => tableIntegral(n, x)
    end
end

function _integrate(ex::Expr, x)
    if ex.head ≠ :call
        ArgumentError("ex must be a call expression") |> throw
    end

    op = ex.args[1]
    @match op begin
        :+ || :- => Expr(:call, op, [integrate(a, x) for a ∈ ex.args[2:end]]...)
        :* || :/ =>
            @match consts(ex) begin
                (0, _)  => 0
                (1, ex) => tableIntegral(ex, x)
                (n, ex) => Expr(:call, :*, n, tableIntegral(ex, x))
            end
        _ => tableIntegral(ex, x)
    end
end

function tableIntegral(ex::Value, x)
    int = :($(ex) * d($(x)))

    for rule ∈ integralTable
        pat = Expr(:call, :*, rule.pat, :(d(x)))
        match = matchex(pat, int; allow_ex = true, exact = false)

        if match ≠ nothing && rule.check(match)
            return rewrite(ex, rule.pat, rule.rpat)
        end
    end

    nothing
end

ext: jl


text:
module PersistentCohomology

using StructArrays
using StructArrays: finduniquesorted
using SparseArrays, LinearAlgebra
using IntervalSets

export Cochain, vietorisrips, persistentcocycles

include("simplexinterface.jl")
include("cochain.jl")
include("vietorisrips.jl")
include("cocycle.jl")
include("itertools.jl")

end # module


text:
using SpecialFunctions

# computes $\frac{J_1(2\pi r)}{2\pi r}$ where $J_1$ - Bessel function of first kind
jinc(r) = real(2*besselj(1,2π*r) / (2π*r))

using Roots

jinc_zeros = fzeros(jinc, 0, 5)
println("jinc_zeros : $jinc_zeros")

using Cairo, Fontconfig
using Gadfly

my_plot = plot(jinc, 0, 5
        , Guide.XLabel("X"), Guide.YLabel("jinc")
        , Guide.xticks(ticks=[0, 0.61, 1.12, 1.62, 2.12, 2.62, 3.12, 3.62, 4.12, 4.62] ))

draw(PDF("images/jinc.pdf", 5inch, 3inch), my_plot)
draw(PNG("images/jinc.png", 5inch, 3inch), my_plot)

max_stars_repo_name: alsam/julia-samples

max_issues_repo_name: alsam/julia-samples

max_forks_repo_name: alsam/julia-samples


text:
##############
#  Polygons  #
##############

struct ConvexPolygon{T} <: Surface
    vertices::Vector{Point{T}}
end

ConvexPolygon(v::Vector{Tuple{T, T}}) where T = ConvexPolygon{T}([SVector(i1, i2) for (i1, i2) in v])

nb_vertices(p::ConvexPolygon) = length(p.vertices)
nb_sides(p::ConvexPolygon) = length(p.vertices)

==(p1::ConvexPolygon, p2::ConvexPolygon) = nb_sides(p1) == nb_sides(p2) && all(p1.vertices .== p2.vertices)

vertices(p::ConvexPolygon) = p.vertices

isempty(p::ConvexPolygon) = length(p.vertices) <= 2

translate(c::ConvexPolygon{T}, v) where T = ConvexPolygon{T}(map(x -> translate(x, v), c.vertices))
rotate(c::ConvexPolygon{T}, ϕ; kw...) where T = ConvexPolygon{T}(map(x -> rotate(x, ϕ; kw...), c.vertices))
scale(c::ConvexPolygon{T}, λ; center=Point(0.0, 0.0)) where T = ConvexPolygon{T}(map(x -> scale(x, λ; center), c.vertices))

function area(c::ConvexPolygon{T}) where T
	if isempty(c)
        return zero(T)*zero(T)
	else
        # Shoelace formula
        a = zero(T)*zero(T)
        last_vertex = c.vertices[end]
        for vertex in c.vertices
            a += last_vertex[1] * vertex[2] - last_vertex[2] * vertex[1]
            last_vertex = vertex
        end
        return abs(a)/2
	end
end

# REPRESENTATION AS INTERSECTION OF HALF-PLANES

function convert(::Type{Intersection{HalfPlane{T}}}, p::ConvexPolygon{U}) where {T, U}
    if isempty(p); return Intersection{HalfPlane{T}}(HalfPlane{T}[]); end
    half_planes = Intersection{HalfPlane{T}}(Vector{HalfPlane{T}}(undef, nb_sides(p)))
    last_i = nb_vertices(p)
    for i in 1:nb_vertices(p)
      half_planes.content[i] = HalfPlane{T}(convert(Point{T}, p.vertices[last_i]), convert(Point{T}, p.vertices[i]), :right)
        last_i = i
    end
    return half_planes
end
convert(::Type{Reunion{Intersection{HalfPlane{T}}}}, c::ConvexPolygon{U}) where {T, U} = convert(Reunion{Intersection{HalfPlane{T}}}, convert(Intersection{HalfPlane{T}}, c))

convert(::Type{ConvexPolygon{T}}, p::ConvexPolygon{U}) where {T, U} = ConvexPolygon{T}(convert.(Point{T}, p.vertices))
promote_rule(::Type{ConvexPolygon{T}}, ::Type{ConvexPolygon{U}}) where {T, U} = ConvexPolygon{promote_type(T, U)}

in(p, c::ConvexPolygon{T}) where T = in(p, convert(Intersection{HalfPlane{T}}, c))
complement(c::ConvexPolygon{T}) where T = complement(convert(Intersection{HalfPlane{T}}, c))

# INTERSECTIONS BETWEEN CONVEX POLYGONS

# Basically Sutherland-Hodgeman algorithm
function intersect(c::ConvexPolygon{T}, h::HalfPlane{U}) where {T, U}
    if isempty(c); return c; end
    S = promote_type(T, U)
    poly_in = ConvexPolygon{S}(Point{S}[])
    previous_corner = c.vertices[end]
    for corner in c.vertices
        if corner in h
            if previous_corner in h
                # Stay in
                push!(poly_in.vertices, corner)
            else
                # Entering
                push!(poly_in.vertices, corner_point(HalfPlane{T}(previous_corner, corner, :right), h))
                push!(poly_in.vertices, corner)
            end
        elseif previous_corner in h
                # Exiting
                push!(poly_in.vertices, corner_point(HalfPlane{T}(previous_corner, corner, :right), h))
        end
        previous_corner = corner
    end
    return poly_in
end
intersect(h::HalfPlane, c::ConvexPolygon) = intersect(c, h)

function intersect(c::ConvexPolygon, ih::Intersection{HalfPlane{T}}) where T
	for h in ih.content
		c = c ∩ h
		if isempty(c)
			break
		end
	end
    return c
end
intersect(ih::Intersection{HalfPlane{T}}, c::ConvexPolygon) where T = intersect(c, ih)

intersect(c1::ConvexPolygon{U}, c2::ConvexPolygon{T})  where {U, T} = intersect(c1, convert(Intersection{HalfPlane{promote_type(U, T)}}, c2))

# REUNION OF CONVEX POLYGONS

convert(::Type{Reunion{ConvexPolygon{T}}}, c::ConvexPolygon{U}) where {T, U} = Reunion{ConvexPolygon{promote_type(T, U)}}([c])
function convert(::Type{Reunion{Intersection{HalfPlane{T}}}}, u::Reunion{ConvexPolygon{U}}) where {T, U}
    S = promote_type(T, U)
    Reunion{Intersection{HalfPlane{S}}}(map(c -> convert(Intersection{HalfPlane{S}}, c), u.content))
end

promote_rule(::Type{ConvexPolygon{T}}, ::Type{Reunion{ConvexPolygon{U}}}) where {T, U} = Reunion{ConvexPolygon{promote_type(T, U)}}

complement(c::Reunion{ConvexPolygon{T}}) where T = complement(convert(Reunion{Intersection{HalfPlane{T}}}, c))

# _cut(c::ConvexPolygon{T}, h::HalfPlane{T}) where T = (c ∩ h, c ∩ complement(h))
# function _cut(c::ConvexPolygon{T}, ih::Intersection{HalfPlane{T}}) where T
# 	for h in ih.content
#         c, rest = _cut(c, h)
# 		if isempty(c)
# 			break
# 		end
# 	end
#     return c
#     (c ∩ h, c ∩ complement(h))
# end

function intersect(c::ConvexPolygon{T}, uh::Reunion{HalfPlane{U}}) where {T, U}
    S = promote_type(T, U)
    u_poly = Reunion{ConvexPolygon{S}}(ConvexPolygon{S}[])
    rest = c
	for h in uh.content
        ci = rest ∩ h
        rest = rest ∩ complement(h)
        if !isempty(ci)
            push!(u_poly.content, ci)
        end
	end
    return u_poly
end
intersect(uh::Reunion{HalfPlane{T}}, c::ConvexPolygon{U}) where {T, U} = intersect(c, uh)

function intersect(c::ConvexPolygon{T}, uih::Reunion{Intersection{HalfPlane{U}}}) where {T, U}
    S = promote_type(T, U)
    u_poly = Reunion{ConvexPolygon{S}}(ConvexPolygon{S}[])
    rest = c
	for ih in uih.content
        ci = rest ∩ ih
        rest = rest ∩ complement(ih)
        if !isempty(ci)
            u_poly = u_poly ∪ ci
        end
	end
    return u_poly
end
intersect(uih::Reunion{Intersection{HalfPlane{T}}}, c::ConvexPolygon{U}) where {T, U} = intersect(c, uih)

function intersect(uc::Reunion{ConvexPolygon{T}}, s::Surface) where T
    u_poly = Reunion{ConvexPolygon{T}}(ConvexPolygon{T}[])
    for ci in uc.content
        cii = ci ∩ s
		if !isempty(cii)
            u_poly = u_poly ∪ cii
		end
    end
    return u_poly
end
intersect(s::Surface, uc::Reunion{ConvexPolygon{T}}) where T = intersect(uc, s)
intersect(uc1::Reunion{ConvexPolygon{T}}, uc2::Reunion{ConvexPolygon{U}}) where {T, U} = intersect(uc1, convert(Reunion{Intersection{HalfPlane{T}}}, uc2))

union(c1::ConvexPolygon{T}, c2::ConvexPolygon{U}) where {T, U} = Reunion{ConvexPolygon{promote_type(T, U)}}([c1, c2])
union(u1::Reunion{ConvexPolygon{T}}, u2::Reunion{ConvexPolygon{U}}) where {T, U} = Reunion{ConvexPolygon{promote_type(T, U)}}(vcat(u1.content, u2.content))

# AREA OF REUNION OF CONVEX POLYGON
disjoint(c::ConvexPolygon) = c
function disjoint(u::Reunion{ConvexPolygon{T}}) where T
    isempty(u) && return u
    l = Reunion{ConvexPolygon{T}}(ConvexPolygon{T}[u.content[1]])
    for c in u.content[2:end]
        rest = c \ l
        if !isempty(rest)
            l = l ∪ rest
        end
    end
    return l
end
area(u::Reunion{ConvexPolygon{T}}) where T = isempty(u) ? zero(T)*zero(T) : sum(area(c) for c in u.content)


text:
println("Loading Libraries onto main node ...")
include("post_process_imports.jl")
println("Finished loading Libraries onto main node ...")

using Distributed
addprocs(20)

@everywhere include("post_process_defs.jl")
include("../Util/Util.jl")

if length(ARGS) > 0
  RESULTS_DIR = ARGS[1]
  AWS_SIM_BUCKET = ARGS[2]
else
  RESULTS_DIR = "results/test-sims"
  AWS_SIM_BUCKET = "s3://cytof-repfam/test-sims"
end

# Collecting paths to output
paths_to_output = ["$(root)/output.bson"
                   for (root, dirs, files) in walkdir(RESULTS_DIR)
                   if "output.bson" in files]

# Post processing
println("Doing post processing...")
@everywhere function makeplots(path_to_output)
  println("Processing: $(path_to_output)")
  pathdir = getpath(path_to_output)
  path_to_simdat = "$(pathdir)/simdat.bson"
  post_process(path_to_output, path_to_simdat=path_to_simdat)
end

println("Make graphs in parallel ...")
status = pmap(makeplots, paths_to_output, on_error=identity)

println("Send results to S3 ...")
Util.s3sync(RESULTS_DIR, AWS_SIM_BUCKET)

println("Success / Failure status: ")
simsucceeded(x) = (x == nothing)
simfailed(x) = !simsucceeded(x)
numsuccess = count(simsucceeded, status)
println("Number of simulations successfully processed: $(numsuccess)")
failures_indices = findall(simfailed, status)
println("Simulations unsuccessfully processed:")
foreach(i -> println(paths_to_output[i]), failures_indices)

println("DONE!")

#= TEST
path_to_output = "results/test-sims/KMCMC2/z1/scale0/output.bson"
path_to_simdat = "results/test-sims/KMCMC2/z1/scale0/simdat.bson"
include("post_process_defs.jl")
post_process(path_to_output, path_to_simdat=path_to_simdat)
=#


text:
using MultivectorAlgebras
using Base.Test: @test, @test_throws

# Multivector2: same typed arguments
@test begin
    l = Multivector2(1.2, 3.4, 5.6, 7.8)
    r = Multivector2{Float64}(Multivector1{Float64}(1.2, 3.4), Multivector1{Float64}(5.6, 7.8))
    l == r
end

# Multivector2: mixed typed arguments
@test begin
    l = Multivector2(1, 2.3, 4, 5.6)
    r = Multivector2{Float64}(Multivector1{Float64}(1.0, 2.3), Multivector1{Float64}(4.0, 5.6))
    l == r
end

# Multivector2: single argument
@test begin
    l = Multivector2(2.3)
    r = Multivector2{Float64}(Multivector1{Float64}(2.3, 0.0), Multivector1{Float64}(0.0, 0.0))
    l == r
end

# Multivector2: two one-dimensionals
@test begin
    l = Multivector2(Multivector1(1, 2), Multivector1(3, 4))
    r = Multivector2(1, 2, 3, 4)
    l == r
end

# Multivector2: one one-dimensionals
@test begin
    l = Multivector2(Multivector1(1, 2))
    r = Multivector2(1, 2, 0, 0)
    l == r
end

# Multivector2: three arguments
@test begin
    l = Multivector2(1, 2, 3)
    r = Multivector2(1, 2, 3, 0)
    l == r
end

# real
@test begin
    l = real(Multivector2(1, 2, 3, 4))
    r = 1
    l == r
end

# unreal: Multivector2
@test begin
    l = unreal(Multivector2(3, 4, 5, 6))
    r = [4, 5, 6]
    l == r
end

# zero: zero-like
@test begin
    l = zero(Multivector2(1, 2, 3, 4))
    r = Multivector2(0, 0, 0, 0)
    l == r
end

# zero: type argument
@test begin
    l = zero(Multivector2{Float64})
    r = Multivector2(0.0, 0.0, 0.0, 0.0)
    l == r
end

# one: one-like
@test begin
    l = one(Multivector2(1, 2, 3, 4))
    r = Multivector2(1, 0, 0, 0)
    l == r
end

# one: type argument
@test begin
    l = one(Multivector2{Float64})
    r = Multivector2(1.0, 0.0, 0.0, 0.0)
    l == r
end

# conj
@test begin
    l = conj(Multivector2(5, 6, 7, 8))
    r = Multivector2(5, -6, -7, -8)
    l == r
end

# Anti-distributivity of conj
@test begin
    a = Multivector2(1, 2, 3, 4)
    b = Multivector2(5, 6, 7, 8)
    l = conj(a * b)
    r = conj(b) * conj(a)
    l == r
end

# cloak
@test begin
    l = cloak(Multivector2(5, 6, 7, 8))
    r = Multivector2(-5, 6, 7, -8)
    l == r
end

# Anti-distributivity of cloak
@test begin
    a = Multivector2(1, 2, 5, 6)
    b = Multivector2(3, 4, 7, 8)
    l = -cloak(a * b)
    r = cloak(a) * cloak(b)
    l == r
end

# cloak: equivalent to -dagger
@test begin
    a = Multivector2(7, 8, 5, 6)
    l = cloak(a)
    r = -dagger(a)
    l == r
end

# dagger
@test begin
    l = dagger(Multivector2(1, 2, 3, 4))
    r = Multivector2(1, -2, -3, 4)
    l == r
end

# Anti-distributivity of dagger
@test begin
    a = Multivector2(1, 2, 5, 6)
    b = Multivector2(3, 4, 7, 8)
    l = dagger(a * b)
    r = dagger(a) * dagger(b)
    l == r
end

# star
@test begin
    l = star(Multivector2(1, 2, 3, 4))
    r = Multivector2(4, -3, 2, 1)
    l == r
end

# Involutivity of conj
@test begin
    l = conj(conj(Multivector2(1, 2, 3, 4)))
    r = Multivector2(1, 2, 3, 4)
    l == r
end

# Involutivity of cloak
@test begin
    l = cloak(cloak(Multivector2(1, 2, 3, 4)))
    r = Multivector2(1, 2, 3, 4)
    l == r
end

# Involutivity of dagger
@test begin
    l = dagger(dagger(Multivector2(1, 2, 3, 4)))
    r = Multivector2(1, 2, 3, 4)
    l == r
end

# star twice is dagger
@test begin
    l = star(star(Multivector2(1, 2, 3, 4)))
    r = dagger(Multivector2(1, 2, 3, 4))
    l == r
end

# +
@test begin
    l = Multivector2(1, 2, 5, 6) + Multivector2(3, 4, 7, 8)
    r = Multivector2(4, 6, 12, 14)
    l == r
end

# +: addition of real
@test begin
    l = 2 + Multivector2(8, 9, 1, 2)
    r = Multivector2(10, 9, 1, 2)
    l == r
end

# +: addition of real
@test begin
    l = Multivector2(8, 9, 2, 1) + 1
    r = Multivector2(9, 9, 2, 1)
    l == r
end

# zero is addition identity
@test begin
    l = Multivector2(5, 6, 1, 2) + zero(Multivector2(5, 6, 1, 2))
    r = Multivector2(5, 6, 1, 2)
    l == r
end

# Commutativity of addition
@test begin
    l = Multivector2(1, 2, 5, 6) + Multivector2(3, 4, 7, 8)
    r = Multivector2(3, 4, 7, 8) + Multivector2(1, 2, 5, 6)
    l == r
end

# Associativity of addition
@test begin
    l = (Multivector2(1, 2, 7, 8) + Multivector2(3, 4, 7, 8)) + Multivector2(5, 6, 7, 8)
    r = Multivector2(1, 2, 7, 8) + (Multivector2(3, 4, 7, 8) + Multivector2(5, 6, 7, 8))
    l == r
end

# -
@test begin
    l = Multivector2(1, 2, 7, 8) - Multivector2(3, 4, 5, 6)
    r = Multivector2(-2, -2, 2, 2)
    l == r
end

# -: subtraction from real
@test begin
    l = 2 - Multivector2(8, 9, 7, 8)
    r = Multivector2(-6, -9, -7, -8)
    l == r
end

# -: subtraction of real
@test begin
    l = Multivector2(8, 9, 7, 8) - 1
    r = Multivector2(7, 9, 7, 8)
    l == r
end

# -: minus
@test begin
    l = -Multivector2(7, 8, 7, 8)
    r = Multivector2(-7, -8, -7, -8)
    l == r
end

# Anti-commutativity of subtraction
@test begin
    l = Multivector2(1, 2, 7, 8) - Multivector2(3, 4, 7, 8)
    r = -(Multivector2(3, 4, 7, 8) - Multivector2(1, 2, 7, 8))
    l == r
end

# Associativity of subtraction
@test begin
    l = (Multivector2(1, 2, 3, 7) - Multivector2(3, 4, 3, 7)) - Multivector2(5, 6, 3, 7)
    r = Multivector2(1, 2, 3, 7) - (Multivector2(3, 4, 3, 7) + Multivector2(5, 6, 3, 7))
    l == r
end

# *
@test begin
    l = Multivector2(2, 3, 6, 7) * Multivector2(4, 5, 8, 9)
    r = Multivector2(8, 22, 40, 40)
    l == r
end

# *: multiplication by real
@test begin
    l = 5 * Multivector2(2, 3, 4, 5)
    r = Multivector2(10, 15, 20, 25)
    l == r
end

# *: multiplication by real
@test begin
    l = Multivector2(3, 4, 5, 6) * 3
    r = Multivector2(9, 12, 15, 18)
    l == r
end

# one is multiplication identity
@test begin
    l = Multivector2(2, 3, 4, 5) * one(Multivector2(2, 3, 4, 5))
    r = Multivector2(2, 3, 4, 5)
    l == r
end

# Non-commutativity of multiplication
@test begin
    l = commutator(Multivector2(1, 2, 3, 4), Multivector2(5, 6, 7, 8))
    r = Multivector2(0, 0, 0, 0)
    l != r
end

# Associativity of multiplication
@test begin
    l = associator(Multivector2(1, 2, 3, 4), Multivector2(5, 6, 7, 8), Multivector2(9, 10, 11, 12))
    r = Multivector2(0, 0, 0, 0)
    l == r
end

# Nilpotency
@test begin
    l = Multivector2(0, 1, 2, 3) * Multivector2(0, 1, 2, 3)
    r = Multivector2(0, 0, 0, 0)
    l == r
end

# Positivity
@test begin
    l = abs2(Multivector2(6, 7, 8, 9)) > 0
    r = true
    l == r
end

# Composition
@test begin
    a = Multivector2(1, 2, 5, 6)
    b = Multivector2(3, 4, 7, 8)
    l = abs2(a) * abs2(b)
    r = abs2(a * b)
    l == r
end

# \
@test begin
    a = Multivector2(1 // 2, 3 // 4, 5 // 6, 7 // 8)
    b = Multivector2(9 // 10, 11 // 12, 13 // 14, 15 // 16)
    c = a * b
    l = a \ c
    r = b
    l == r
end

# \: zero divisor error
@test_throws ErrorException begin
    Multivector2(0, 1, 3, 4) \ Multivector2(2, 3, 9, 8)
end

# \: zero divisor error
@test_throws ErrorException begin
    Multivector2(0, 1, 2, 3) \ 2
end

# /
@test begin
    a = Multivector2(1 // 2, 3 // 4, 5 // 6, 7 // 8)
    b = Multivector2(9 // 10, 11 // 12, 13 // 14, 15 // 16)
    c = a * b
    l = c / b
    r = a
    l == r
end

# /: zero divisor error
@test_throws ErrorException begin
    Multivector2(1, 2) / Multivector2(0, 3)
end

# /: zero divisor error
@test_throws ErrorException begin
    2 / Multivector2(0, 3)
end

# /: division by real
@test begin
    l = Multivector2(1.2, 3.4) / 5
    r = Multivector2(1.2 / 5, 3.4 / 5)
    l == r
end

# /: division of real
@test begin
    l = 2 / Multivector2(3 // 4, 5 // 6)
    r = 2 * inv(Multivector2(3 // 4, 5 // 6))
    l == r
end

# /: zero denominator
@test_throws ErrorException begin
    Multivector2(1, 2, 3, 4) / 0
end

# \: division by real
@test begin
    l = 5 \ Multivector2(1.2, 3.4, 5.6, 7.8)
    r = Multivector2(5 \ 1.2, 5 \ 3.4, 5 \ 5.6, 5 \ 7.8)
    l == r
end

# \: division of real
@test begin
    l = Multivector2(3 // 4, 5 // 6, 7 // 8, 9 // 10) \ 2
    r = inv(Multivector2(3 // 4, 5 // 6, 7 // 8, 9 // 10)) * 2
    l == r
end

# \: zero denominator
@test_throws ErrorException begin
    0 \ Multivector2(1, 2, 3, 4)
end

# inv
@test begin
    a = Multivector2(1 // 2, 3 // 4, 7 // 8, 9 // 10)
    b = Multivector2(5 // 6, 7 // 8, 7 // 8, 9 // 10)
    l = a / b
    r = inv(b / a)
    l == r
end

# inv: zero divisor error
@test_throws ErrorException begin
    inv(Multivector2(0, 1, 2, 3))
end

# isreal: false case
@test begin
    l = isreal(Multivector2(1, 2, 3, 4))
    r = false
    l == r
end

# isreal: true case
@test begin
    l = isreal(Multivector2(3, 0, 0, 0))
    r = true
    l == r
end

# asarray
@test begin
    l = asarray(Multivector2(4, 5, 6, 7))
    r = [4, 5, 6, 7]
    l == r
end

# iszero: false case
@test begin
    l = iszero(Multivector2(1, 2, 3, 4))
    r = false
    l == r
end

# iszero: true case
@test begin
    l = iszero(Multivector2(0, 0, 0, 0))
    r = true
    l == r
end

# abs: negative case
@test begin
    l = abs(Multivector2(-2, 3, 5, 4))
    r = 2
    l == r
end

# abs: positive case
@test begin
    l = abs(Multivector2(2, 3, 4, 5))
    r = 2
    l == r
end

# abs2
@test begin
    l = abs2(Multivector2(4, 5, 8, 9))
    r = 16
    l == r
end

# iszerodivisor: false case
@test begin
    l = iszerodivisor(Multivector2(2, 3, 4, 5))
    r = false
    l == r
end

# iszerodivisor: true case
@test begin
    l = iszerodivisor(Multivector2(0, 3, 5, 6))
    r = true
    l == r
end

# selfstar
@test begin
    a = 1.2
    b = 3.4
    l = selfstar(Multivector2(a, 5.0, 6.0, b))
    r = Multivector2((a + b) / 2, 0.0, 0.0, (a + b) / 2)
    l == r
end

# Involutivity of selfstar
@test begin
    z = Multivector2(1 // 2, 3 // 4, 5 // 6, 7 // 8)
    l = selfstar(selfstar(z))
    r = selfstar(z)
    l == r
end

# antiselfstar
@test begin
    a = 1 // 2
    b = 3 // 4
    l = antiselfstar(Multivector2(a, 3, 4, b))
    r = Multivector2((a - b) / 2, 0, 0, (b - a) / 2)
    l == r
end

# Involutivity of antiselfstar
@test begin
    z = Multivector2(1 // 2, 3 // 4, 5 // 6, 7 // 8)
    l = antiselfstar(antiselfstar(z))
    r = antiselfstar(z)
    l == r
end

# show
@test begin
    io = IOBuffer()
    show(io, Multivector2(1, 2, 3, 4))
    l = String(take!(io))
    r = "[1: 1, A: 2, B: 3, AB: 4]"
    l == r
end

# crossratioL and möbiusR
@test begin
    a = 1
    b = 2
    c = 3
    d = 4
    x1 = Multivector2(1 // 2, 3 // 4, 5 // 6, 7 // 8)
    x2 = Multivector2(9 // 2, 10 // 4, 11 // 6, 12 // 8)
    x3 = Multivector2(13 // 2, 14 // 4, 15 // 6, 16 // 8)
    x4 = Multivector2(17 // 2, 18 // 4, 19 // 6, 20 // 8)
    y1 = möbiusR(x1, a, b, c, d)
    y2 = möbiusR(x2, a, b, c, d)
    y3 = möbiusR(x3, a, b, c, d)
    y4 = möbiusR(x4, a, b, c, d)
    l = real(crossratioL(x1, x2, x3, x4))
    r = real(crossratioL(y1, y2, y3, y4))
    l == r
end

# crossratioR and möbiusL
@test begin
    a = 1
    b = 2
    c = 3
    d = 4
    x1 = Multivector2(1 // 2, 3 // 4, 5 // 6, 7 // 8)
    x2 = Multivector2(9 // 2, 10 // 4, 11 // 6, 12 // 8)
    x3 = Multivector2(13 // 2, 14 // 4, 15 // 6, 16 // 8)
    x4 = Multivector2(17 // 2, 18 // 4, 19 // 6, 20 // 8)
    y1 = möbiusL(x1, a, b, c, d)
    y2 = möbiusL(x2, a, b, c, d)
    y3 = möbiusL(x3, a, b, c, d)
    y4 = möbiusL(x4, a, b, c, d)
    l = real(crossratioR(x1, x2, x3, x4))
    r = real(crossratioR(y1, y2, y3, y4))
    l == r
end

# crossratioL and möbiusR
@test begin
    a = Multivector2(98, -42, -13, 118)
    b = Multivector2(-54, 102, 86, 6)
    c = Multivector2(-126, -18, 33, -74)
    d = Multivector2(58, -86, 0, 20)
    x1 = Multivector2(Int128(1) // 2, 3 // 4, 5 // 6, 7 // 8)
    x2 = Multivector2(Int128(9) // 2, 10 // 4, 11 // 6, 12 // 8)
    x3 = Multivector2(Int128(13) // 2, 14 // 4, 15 // 6, 16 // 8)
    x4 = Multivector2(Int128(17) // 2, 18 // 4, 19 // 6, 20 // 8)
    y1 = möbiusR(x1, a, b, c, d)
    y2 = möbiusR(x2, a, b, c, d)
    y3 = möbiusR(x3, a, b, c, d)
    y4 = möbiusR(x4, a, b, c, d)
    l = real(crossratioL(x1, x2, x3, x4))
    r = real(crossratioL(y1, y2, y3, y4))
    l == r
end

# crossratioR and möbiusL
@test begin
    a = Multivector2(98, -42, -13, 118)
    b = Multivector2(-54, 102, 86, 6)
    c = Multivector2(-126, -18, 33, -74)
    d = Multivector2(58, -86, 0, 20)
    x1 = Multivector2(Int128(1) // 2, 3 // 4, 5 // 6, 7 // 8)
    x2 = Multivector2(Int128(9) // 2, 10 // 4, 11 // 6, 12 // 8)
    x3 = Multivector2(Int128(13) // 2, 14 // 4, 15 // 6, 16 // 8)
    x4 = Multivector2(Int128(17) // 2, 18 // 4, 19 // 6, 20 // 8)
    y1 = möbiusL(x1, a, b, c, d)
    y2 = möbiusL(x2, a, b, c, d)
    y3 = möbiusL(x3, a, b, c, d)
    y4 = möbiusL(x4, a, b, c, d)
    l = real(crossratioR(x1, x2, x3, x4))
    r = real(crossratioR(y1, y2, y3, y4))
    l == r
end

# random
@test begin
    l = random(Multivector2{Int})
    r = random(Multivector2{Int})
    l != r
end


text:
#

# Task
# move linear at constant velocity along x axis, 10m/s, starting from 0

# Track the state separate like a filter by using underlying Factor Graph operations instead.
# the next example will keep the data in the same factor, and achieve objective more efficiently.

using IncrementalInference

# state
# ContinuousScalar

# Prediction
# LinearConditional

# measurement
# Prior

function forwardConvolve(X0::Array{Float64,2}, model)

  fg = initfg()

  addVariable!(fg, :x0, ContinuousScalar)
  initManual!(fg, :x0, X0)
  addVariable!(fg, :x1, ContinuousScalar)

  addFactor!(fg, [:x0;:x1], model)

  ## TODO -- dont use name here, add API to just use z2 here
  return approxConv(fg, :x0x1f1, :x1)
end

## Start

# first numerical values -- samples from the marginal of X0
z1 = Normal(0,0.1)
X0 = rand(z1, 1,100)

## predict -- project / conv
# 0 -> 1 seconds
# make approx function
z2 = Normal(11,1.0) # odo
statemodel = LinearConditional( z2 )
X1_ = forwardConvolve(X0, statemodel)

## measure -- product of beliefs, using `ApproxManifoldProducts.jl`

predX1 = manikde!(X1_, ContinuousScalar)
z3 = Normal(9.5,0.75)
measX1 = manikde!(reshape(rand(z3,100),1,:), ContinuousScalar)

# do actual product
posterioriX1 = predX1 * measX1
X1 = getPoints(posterioriX1)

## predict, 1->2 seconds
z4 = Normal(8,2.0) # odo
statemodel = LinearConditional( z4 )
X2_ = forwardConvolve(X1, statemodel)

# measure (product)
z5 = Normal(20.8,2.0)
measX2 = manikde!(reshape(rand(z5,100),1,:), ContinuousScalar)
posterioriX2 = predX2 * measX2
X2 = getPoints(posterioriX2)

# predict
# ...

# measure

# predict

# measure

#


text:
using LinearAlgebra

export Domain

# Geometric tolerance
geps=1e-18

"""
Struct to define a finite element domain

* `interp`: `InterpolationFunction` instance
* `restrictions`: definitions of domain limits in terms of auxiliary coordinates (`interp.coords`)
"""
mutable struct Domain
    interp::InterpolationFunction
    restrictions::Vector{
        Tuple{Any, Any}
    }
    _restriction_funcs::Vector{
        Tuple{Function, Function}
    }
    _restriction_derivatives::Vector{
        Tuple{Function, Function}
    }
    _regressor::Regressor
end

"""
Function to obtain restriction function
"""
function _get_restrfunc(re::Union{Expr, Number}, coords::Vector{Symbol}, ncoord::Int64)
    q=quote
        bfun=($(coords...),) -> $(re)
    end

    try
        eval(q)

        return bfun
    catch
        throw(error("_get_restrfunc:ERROR:unable to obtain restriction function"))
    end
end

"""
Function to obtain restriction gradient function
"""
function _get_gradfunc(re::Expr, coords::Vector{Symbol}, ncoord::Int64)
    der_exprs=:([$(differentiate(re, coords)...)])

    q=quote
        bfun=($(coords...),) -> $(der_exprs)
    end

    try
        eval(q)

        return bfun
    catch
        throw(error("_get_gradfunc:ERROR:unable to obtain restriction gradient. Please check for correctness with respect to doc. directives"))
    end
end

"""
Alternative implementation for literals
"""
function _get_gradfunc(re::Number, coords::Vector{Symbol}, ncoord::Int64)
    q=quote
        bfun=($(coords...),) -> ( 0.0 )
    end

    try
        eval(q)

        return bfun
    catch
        throw(error("_get_gradfunc:ERROR:error parsing restriction expression of literal type"))
    end
end

"""
Constructor for Domain struct

* `interp`: `InterpolationFunction` instance
* `restrictions`: definitions of domain limits in terms of auxiliary coordinates (`interp.coords`),
    which are negative within the domain and positive within it
"""
function Domain(
    interp::InterpolationFunction,
    restrictions,
    points::Vector{Vector{Float64}}
)
    ncoords=length(restrictions)

    _restriction_derivatives=[
        (_get_gradfunc(:($(vr[1])-$(c)), interp.coords, i), _get_gradfunc(:($(c)-$(vr[2])), interp.coords, i))
            for (i, (c, vr)) in enumerate(zip(interp.coords, restrictions))
    ]

    _restriction_functions=[
        (_get_restrfunc(vr[1], interp.coords, i), _get_restrfunc(vr[2], interp.coords, i))
            for (i, vr) in enumerate(restrictions)
    ]

    _restrictions=[
        Tuple{Any, Any}(vr)
                for vr in restrictions
    ]

    Domain(
        interp,
        _restrictions,
        _restriction_functions,
        _restriction_derivatives,
        Regressor(
            interp,
            points
        )
    )
end

"""
Same as `regress(regr::Regressor, Vector{Float64})`, but using a domain's regressor
"""
function regress(dmn::Domain, vals::Vector{Float64})
    return regress(dmn._regressor, vals)
end

"""
Get projection of a vector onto another (`u` onto `v`)
"""
function _getproj(u::Vector{Float64}, v::Vector{Float64})
    v*(dot(u, v)/(norm(v)^2+geps))
end

"""
Get determinant module by orthonormalizing the matrix at hand
"""
function orthonormal_det!(M::Matrix{Float64})
    d=1.0

    nd=size(M, 2)

    locn=1.0
    for i=1:nd
        for j=1:(i-1)
            M[:, i].-=dot(M[:, j], M[:, i]).*M[:, j]
        end

        locn=norm(M[:, i])+geps

        M[:, i]./=locn

        d*=locn
    end

    d
end

"""
Get component normal to the Jacobian space
"""
function normal_volspace!(M::Matrix{Float64}, nrestr::Int64, restr_derivs, isfirst::Bool)
    d=1.0

    nd=size(M, 2)

    locn=1.0
    for i=1:nd
        if i!=nrestr
            M[:, i].-=M[:, nrestr].*(restr_derivs[i]/restr_derivs[nrestr])

            for j=1:(i-1)
                if j!=nrestr
                    M[:, i].-=dot(M[:, j], M[:, i]).*M[:, j]
                end
            end

            locn=norm(M[:, i])+geps

            M[:, i]./=locn

            d*=locn
        end
    end

    for i=1:nd
        if i!=nrestr
            M[:, nrestr].-=dot(M[:, i], M[:, nrestr]).*M[:, i]
        end
    end

    M[:, nrestr].*=d/(norm(M[:, nrestr])+geps)

    if isfirst
        @. M[:, nrestr]=-M[:, nrestr]
    end

    M[:, nrestr]
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.isomorphism.vf2" begin

@testset "subgraph" begin
    null = immutableplaingraph()
    @test !isgraphmatch(null, null)
    @test !issubgraphmatch(null, null)
    @test isempty(graphmatches(null, null))
    @test isempty(subgraphmatches(null, null))

    p5 = pathgraph(5)
    k5 = completegraph(5)
    k6 = completegraph(6)
    @test !issubgraphmatch(p5, k5)
    @test issubgraphmatch(k6, k5)
    @test !issubgraphmatch(k5, k6)
    @test isgraphmatch(k5, k5)
    @test !isgraphmatch(k6, k5)

    p7 = pathgraph(7)
    disconn1 = immutableplaingraph(6, [(1,2), (2,3), (4,5), (5,6)])
    @test issubgraphmatch(p7, disconn1)
    disconn2 = immutableplaingraph(6, [(1,2), (3,4), (5,6)])
    @test !issubgraphmatch(p7, disconn2)

    cd6 = circularladder(6)
    mb6 = moebiusladder(6)
    @test !isgraphmatch(cd6, mb6)
    # TODO: relabeled mb6
    # @test isgraphmatch(rmb6, mb6)
end

@testset "edgesubgraph" begin
    null = immutableplaingraph()
    @test !isedgesubgraphmatch(null, null)
    @test isempty(edgesubgraphmatches(null, null))

    c4 = cyclegraph(4)
    k4 = completegraph(4)
    star = immutableplaingraph(6, [(1,2), (1,3), (1,4), (1,5), (1,6)])
    @test isedgesubgraphmatch(k4, c4)
    @test !isedgesubgraphmatch(k4, star)

    p4 = pathgraph(4)
    disconn1 = immutableplaingraph(4, [(1,2), (3,4)])
    @test isedgesubgraphmatch(p4, disconn1)
    disconn2 = immutableplaingraph(5, [(1,2), (3,4), (4,5)])
    @test !isedgesubgraphmatch(p4, disconn2)

    # Delta-Y exchange
    tri = cyclegraph(3)
    star = immutableplaingraph(4, [(1,2), (1,3), (1,4)])
    @test isedgesubgraphmatch(tri, tri)
    @test isedgesubgraphmatch(star, star)
    @test !isedgesubgraphmatch(tri, star)
    @test !isedgesubgraphmatch(star, tri)

    k4 = completegraph(4)
    〼 = immutableplaingraph(4, [(1,2), (2,3), (1,3), (3,4), (4,2)])
    butterfly = immutableplaingraph(
        5, [(1,2), (2,3), (1,3), (3,4), (4,5), (5,3)])
    @test !isedgesubgraphmatch(k4, butterfly)
    matches = edgesubgraphmatches(k4, 〼)
    @test length(collect(matches)) == 24
end

@testset "mandatory" begin
    path = pathgraph(7)
    subp = pathgraph(3)
    eiso = edgesubgraphmatches(path, subp)
    @test length(collect(eiso)) == 10

    restricted = edgesubgraphmatches(path, subp, mandatory=Dict(3 => 1))
    @test length(collect(restricted)) == 2
end

end # graph.isomorphism.vf2


text:
using Revise
using HPD

path = joinpath(dirname(pathof(HPD)), "..", "examples")
domain = load_domain(joinpath(path, "domain.hpd"))
problem = load_problem(joinpath(path, "problem.hpd"))
state = initstate(domain, problem)
goal = goalstate(domain, problem)
1


text:
"""
Replaces an expression value in the expression container if the key exists
"""
function set_expression!(
    container::OptimizationContainer,
    ::Type{S},
    cost_expression::JuMP.AbstractJuMPScalar,
    component::T,
    time_period::Int,
) where {S <: CostExpressions, T <: PSY.Component}
    if has_container_key(container, S, T)
        device_cost_expression = get_expression(container, S(), T)
        component_name = PSY.get_name(component)
        device_cost_expression[component_name, time_period] = cost_expression
    end
    return
end


text:
"""
	$(SIGNATURES)

Constructors for ParamVector
"""
function ParamVector(id :: ObjectId)
    return ParamVector(objId = id)
end

# With default parameter transformation
# Providing vector of Param or an OrderedDict.
function ParamVector(id :: ObjectId, pv)
    return ParamVector(id,
        make_ordered_dict(pv),
        LinearTransformation{ValueType}())
end

function validate_pvec(pvec :: ParamVector; silent = true)
    isValid = true;
    for p in pvec
        isValid = isValid &&  validate(p; silent);
    end
    return isValid
end

function make_ordered_dict(v)
    d = OrderedDict{Symbol, AbstractParam}();
    for p in v
        d[p.name] = p;
    end
    return d
end

make_ordered_dict(v :: OrderedDict{Symbol, AbstractParam}) = v;

function Base.show(io :: IO,  pvec :: ParamVector)
    n = length(pvec);
    idStr = make_string(pvec.objId);
    print(io,  "ParamVector $idStr of length $n");
    # if n > 0
    #     for p in pvec
    #         println(io, "  ", p);
    #     end
    # end
    return nothing
end

ModelObjectsLH.get_object_id(pv :: ParamVector) = pv.objId;
# ModelObjectsLH.has_object_id(pv :: ParamVector) = true;
Base.length(pvec :: ParamVector) = Base.length(pvec.pv);
Base.isempty(pvec :: ParamVector) = Base.isempty(pvec.pv);

# Check that param vector matches model object
function object_id_matches(pvec :: ParamVector, objId :: ObjectId)
    return isequal(pvec.objId, objId)
end

## ----------  Retrieve

# Allows to access values simply as `p[j]`.
# But access by name is far more efficient.
function getindex(pvec :: ParamVector, j :: Integer)
    @argcheck j <= Base.length(pvec.pv);
    k = collect(keys(pvec.pv))[j];
    return pvec.pv[k]
end

function Base.iterate(pvec :: ParamVector)
    if isempty(pvec)
        return nothing
    else
        return pvec[1], 1
    end
end

function Base.iterate(pvec :: ParamVector, s)
    if s >= length(pvec)
        return nothing
    else
        return pvec[s+1], s+1
    end
end

Base.eltype(pvec :: ParamVector) = Param;

"""
    $(SIGNATURES)

Returns a named parameter and its index in the `ParamVector`.
First occurrence. Returns `nothing` if not found.
"""
function retrieve(pvec :: ParamVector, pName :: Symbol)
    return get(pvec.pv, pName, nothing);
end

# param_index(pvec :: ParamVector, pName :: Symbol) =
#     findfirst(x -> x == pName,  collect(keys(pvec.pv)));
    # findfirst(x -> isequal(x.name, pName),  pvec.pv);

"""
	$(SIGNATURES)

Is this parameter calibrated?
"""
function is_calibrated(pvec :: ParamVector, pName :: Symbol)
    @assert param_exists(pvec, pName)  "Not found: $pName";
    return is_calibrated(retrieve(pvec, pName))
end

"""
	$(SIGNATURES)

Does a `Param` named `pName` exist in `ParamVector pvec`?
"""
param_exists(pvec :: ParamVector, pName :: Symbol) =
    haskey(pvec.pv, pName);

"""
	$(SIGNATURES)

Return the value of a parameter from a `ParamVector`.
Returns `nothing` if not found.
"""
function param_value(pvec :: ParamVector, pName :: Symbol)
    p = retrieve(pvec, pName);
    if isnothing(p)
        return nothing
    else
        return value(p)
    end
end

"""
	$(SIGNATURES)

Return the default value of a parameter from a `ParamVector`.
Returns `nothing` if not found.
"""
function param_default_value(pvec :: ParamVector, pName :: Symbol)
    p = retrieve(pvec, pName);
    if isnothing(p)
        return nothing
    else
        return default_value(p)
    end
end

function calibrated_value(pvec :: ParamVector, pName :: Symbol)
    return calibrated_value(retrieve(pvec, pName));
end

# """
# 	$(SIGNATURES)

# Return indices of all parameters with a given calibration status.
# """
# function indices_calibrated(pvec :: ParamVector, isCalibrated :: Bool)
#     n = length(pvec);
#     if n > 0
#         isCalV = [pvec.pv[k].isCalibrated  for k in keys(pvec.pv)];
#         # isCalV = [pvec.pv[j].isCalibrated  for j in 1 : n];
#         idxV = findall(isCalV .== isCalibrated);
#     else
#         idxV = Vector{Int}();
#     end
#     return idxV
# end

"""
	$(SIGNATURES)

List of calibrated or not calibrated parameters. Returns vector of `Param`.
"""
function calibrated_params(pvec :: ParamVector; isCalibrated :: Bool = true)
    return all_params(pvec; isCalibrated)
end

"""
	$(SIGNATURES)

Vector of all parameters, including nested objects. Optionally filter by calibration status.
"""
function all_params(pvec :: ParamVector; isCalibrated = nothing)
    pList = Vector{AbstractParam}();
    for (k, p) in pvec.pv
        if isnothing(isCalibrated)  ||   (is_calibrated(p) == isCalibrated)
            push!(pList, p);
        end
    end
    return pList
end

"""
    $(SIGNATURES)

Number of calibrated or fixed parameters and their total element count.
"""
function n_calibrated_params(pvec :: ParamVector; isCalibrated :: Bool = true)
    pList = calibrated_params(pvec; isCalibrated);
    nParams = length(pList);
    nElem = 0;
    for p in pList
        if isCalibrated
            pVal = calibrated_value(p);
            @assert !isnothing(pVal)  "calibrated value nothing: $pvec / $p";
            nCal = length(pVal);
        else
            nCal = 0;
        end
        if isCalibrated
            nElem += nCal;
        else
            nElem += (length((value(p))) - nCal);
        end
    end
    return nParams, nElem
end

## ------------  Modify

"""
    $(SIGNATURES)

Append a `Param` to a `ParamVector`
"""
function append!(pvec :: ParamVector,  p :: AbstractParam)
    @assert !param_exists(pvec, p.name)  "$(p.name) already exists";
    pvec.pv[p.name] = p;
    return nothing
end

"""
	$(SIGNATURES)

Remove the parameter names `pName` from `pvec`.
"""
function remove!(pvec :: ParamVector, pName :: Symbol)
    @assert haskey(pvec.pv, pName)  "$pName does not exist";
    delete!(pvec.pv, pName);
    return nothing
end

"""
	$(SIGNATURES)

Replace a parameter with a new parameter `p`. Without changing the order.
"""
function replace!(pvec :: ParamVector, p :: AbstractParam)
    pvec.pv[p.name] = p;
    return nothing
end

function set_calibration_status_all_params!(pvec :: ParamVector, isCalibrated :: Bool)
    for p in calibrated_params(pvec; isCalibrated = !isCalibrated)
        set_calibration_status!(p, isCalibrated);
    end
end

"""
	$(SIGNATURES)

Set default values for all parameters to values.
Useful for fixing parameters at previously calibrated values.
"""
function set_default_values_all_params!(pvec :: ParamVector)
    for p in pvec
        set_default_value!(p, value(p));
    end
end

function set_default_value!(pvec :: ParamVector, pName :: Symbol, v)
    p = retrieve(pvec, pName);
    set_default_value!(p, v);
end

"""
	$(SIGNATURES)

Set whether or not a parameter is calibrated.
"""
function change_calibration_status!(pvec :: ParamVector, pName :: Symbol,
    doCal :: Bool)

    @assert param_exists(pvec, pName)  "$pName does not exist";
    p = retrieve(pvec, pName);
    if doCal
        calibrate!(p);
    else
        fix!(p);
    end
end

calibrate!(pvec :: ParamVector, pName :: Symbol) =
    change_calibration_status!(pvec, pName, true);
fix!(pvec :: ParamVector, pName :: Symbol) =
    change_calibration_status!(pvec, pName, false);

"""
	$(SIGNATURES)

Change the value of parameter `pName`.
"""
function change_value!(pvec :: ParamVector, pName :: Symbol, newValue;
    skipInvalidSize :: Bool = false
    )
    @assert param_exists(pvec, pName)  "$pName does not exist";
    p = retrieve(pvec, pName);
    oldValue = set_value!(p, newValue; skipInvalidSize);

    # if size(default_value(p)) == size(newValue)
    #     oldValue = set_value!(p, newValue; skipInvalidSize);
    # else
    #     @warn("""
    #         Wrong size for $p in $pvec
    #         Given: $(size(newValue))
    #         Expected: $(size(default_value(p)))
    #         """);
    #     if skipInvalidSize
    #         oldValue = value(p);
    #     else
    #         error("Stopped");
    #     end
    # end
    return oldValue
end

## ------------------  Reporting

"""
    $(SIGNATURES)

Reports calibrated (or fixed) parameters for one ParamVector as a `PrettyTable`.
"""
function report_params(pvec :: ParamVector, isCalibrated :: Bool;
    io :: IO = stdout,  closeToBounds :: Bool = false)

    objId = make_string(pvec.objId);
    dataM = param_table(pvec, isCalibrated;  closeToBounds = closeToBounds);

    if isnothing(dataM)
        println(io, "\t$objId:  Nothing to report");
    else
        # Reports description, name, value (not symbol)
        pretty_table(io,
            [get_descriptions(dataM) get_names(dataM) get_values(dataM)];
            header = [objId, " ", " "]);
    end
    return nothing
end

"""
	$(SIGNATURES)

Table with calibrated (or fixed) parameters.
Optionally reports params that are close to bounds.
Columns are name, description, value, symbol.
This is for reporting during a computation.
"""
function param_table(pvec :: ParamVector, isCalibrated :: Bool;
    closeToBounds :: Bool = false)

    if closeToBounds
        pList = find_close_to_bounds(pvec);
    else
        pList = calibrated_params(pvec; isCalibrated);
    end
    n = length(pList);

    if isempty(pList)
        dataM = nothing;
    else
        # dataM = Matrix{String}(undef, n, 4);
        dataM = ParamTable(n);
        for (j, p) in enumerate(pList)
            set_row!(dataM, j, string(name(p)), lsymbol(p),
                p.description, formatted_value(value(p)));
        end
    end
    return dataM
end

# Find parameters that are close to bounds
function find_close_to_bounds(pvec :: ParamVector; rtol = 0.01)
    pList = calibrated_params(pvec);
    pCloseV = Vector{AbstractParam}();
    for p in pList
        if close_to_bounds(p; rtol = rtol)
            push!(pCloseV, p);
        end
    end
    return pCloseV
end

## -------------  Compare

function compare_params(pvec1 :: ParamVector, pvec2 :: ParamVector;
    ignoreCalibrationStatus :: Bool = true)

    pMiss1 = missing_params(pvec1, pvec2);
    pMiss2 = missing_params(pvec2, pvec1);
    pDiff = find_param_diffs(pvec1, pvec2; ignoreCalibrationStatus = true);
    return pMiss1, pMiss2, pDiff
end

# List params that are in p2V but not in p1V
function missing_params(pvec1 :: ParamVector, pvec2 :: ParamVector)
    missList = Vector{Symbol}();
    for p in pvec2
        pName = name(p);
        if !param_exists(pvec1, pName)
            push!(missList, pName);
        end
    end
    return missList
end

# List all params that are in both vectors, but differ
function find_param_diffs(pvec1 :: ParamVector, pvec2 :: ParamVector;
    ignoreCalibrationStatus :: Bool = true)

    diffList = Dict{Symbol, Any}();
    for p1 in pvec1
        pName = name(p1);
        if param_exists(pvec2, pName)
            diffs = param_diffs(p1, retrieve(pvec2, pName); ignoreCalibrationStatus);
            if !isempty(diffs)
                diffList[pName] = diffs;
            end
        end
    end
    return diffList
end

param_diffs(p1, p2) = [:type];

function param_diffs(p1 :: T1, p2 :: T1;
    ignoreCalibrationStatus :: Bool = true) where T1 <: AbstractParam

    diffs = Vector{Symbol}();
    if !isequal(value(p1), value(p2))
        push!(diffs, :value);
    end
    if !ignoreCalibrationStatus
        if is_calibrated(p1) != is_calibrated(p2)
            push!(diffs, :calibration);
        end
    end
    return diffs
end

## --------------  Dicts and Vectors

"""
    $(SIGNATURES)

Collect calibrated values or default values into Dict.
Used to go back and forth between guess and model parameters.
All values are stored, even if only some are calibrated.

# Arguments
- `valueType`: `:defaultValue`, `:value`, or `:calibratedValue`. Determines which of those is stored.
"""
function make_dict(pvec :: ParamVector; isCalibrated :: Bool, valueType :: Symbol)

    pd = Dict{Symbol, Any}()
    pList = calibrated_params(pvec; isCalibrated);
    for p in pList
        if valueType == :value
            v = value(p);
        elseif valueType == :defaultValue
            v = default_value(p);
        elseif valueType == :calibratedValue
            v = calibrated_value(p);
        else
            error("Unknown valueType: $valueType");
        end
        # Copy is important. Otherwise changing the values later will change the dict.
        pd[p.name] = copy(v);
    end
    return pd
end

# Typically, useValues when calibrated; defaults otherwise
# make_dict(pvec :: ParamVector; isCalibrated :: Bool) =
#     make_dict(pvec; isCalibrated, useValues = isCalibrated)

# """
#     $(SIGNATURES)

#     Make a vector of values into a Dict

#     The inverse of `make_guess`.
#     Used to go back from vector to model parameters.

#     Undoes the parameter transformation from `make_guess`.

#     OUT
#         pd :: Dict
#             maps param names (symbols) to values
#         iEnd :: Integer
#             last element of `v` used up
# """
# function vector_to_dict(pvec :: ParamVector, vVec :: ValueVector,
#     isCalibrated :: Bool; startIdx = 1)

#     n = length(pvec);
#     @assert n > 0  "$pvec vector is empty"
#     pList = calibrated_params(pvec, isCalibrated);
#     v = get_values(vVec);
#     pNameV = pnames(vVec);

#     pd = Dict{Symbol, Any}();
#     # Last index of `v` used so far
#     iEnd = startIdx - 1;
#     for p in pList
#         nElem = length(default_value(p));
#         if isa(default_value(p), AbstractFloat)
#             vIdxV = iEnd + 1;
#             pValue = v[vIdxV];
#         elseif isa(default_value(p), AbstractArray)
#             vIdxV = iEnd .+ (1 : nElem);
#             pValue = reshape(v[vIdxV], size(default_value(p)));
#         else
#             pType = typeof(default_value(p));
#             error("Unexpected type: $pType")
#         end
#         iEnd += nElem;
#         pd[name(p)] = untransform_param(pvec.pTransform, p, pValue);
#         @assert all(isequal.(name(p),  pNameV[vIdxV]))  """
#             Name mismatch:  $(get_object_id(pvec)) $(name(p))
#             Values: $pValue
#             Indices: $vIdxV
#             Names expected: $(pNameV[vIdxV])
#             $startIdx
#             """
#     end
#     return pd, iEnd
# end

"""
    $(SIGNATURES)

Set values in param vector from dictionary.
"""
function set_values_from_dict!(pvec :: ParamVector, d :: D1;
    skipInvalidSize :: Bool = false) where D1 <: AbstractDict
    for (pName, newValue) in d
        change_value!(pvec, pName :: Symbol, newValue;  skipInvalidSize);
    end
    return nothing
end

"""
	$(SIGNATURES)

Set values in `pvecOld` from another `ParamVector` `pvecNew`.
Only for values that are in both `ParamVector`s and that are `isCalibrated` in both.
Only if the size matches.
"""
function set_own_values_from_pvec!(pvecOld :: ParamVector,  pvecNew :: ParamVector,
    isCalibrated :: Bool;
    skipInvalidSize :: Bool = false
    )

    for (pName, p) in pvecNew.pv
        if param_exists(pvecOld, pName)  &&  is_calibrated(p)
            pOld = retrieve(pvecOld, pName);
            if is_calibrated(pOld)
                newValue = calibrated_value(p);
                set_calibrated_value!(pOld, newValue; skipInvalidSize);
            end
        end
    end

    # dOld = make_dict(pvecOld; isCalibrated, useValues = true);
    # dNew = make_dict(pvecNew; isCalibrated, useValues = true);
    # pNameV = intersect(keys(dOld), keys(dNew));
    # for pName in pNameV
    #     change_value!(pvecOld, pName, dNew[pName]; skipInvalidSize);
    # end
    return nothing
end

## -----------  Testing

# Make a ParamVector for testing
# Alternates between calibrated and fixed parameters
function make_test_pvector(n :: Integer; objId :: Symbol = :obj1,
    offset :: Float64 = 0.0)

    pv = ParamVector(ObjectId(objId));
    for i1 = 1 : n
        if i1 == 1
            p = make_test_cal_array(:p1, 2; offset);
        else
            p = init_parameter(i1; offset = offset);
        end
        if isodd(i1)
            calibrate!(p);
        end
        ModelParams.append!(pv, p);
    end
    return pv
end

function init_parameter(i1 :: Integer;  offset :: Float64 = 0.0)
    pSym = Symbol("p$i1");
    pName = "param$i1";
    pDescr = "sym$i1";
    valueM = (offset + i1) .+ collect(1 : i1) * [1.0 2.0];
    return Param(pSym, pName, pDescr, valueM)
end

# ------------------


text:
# Those will be necessary for all test related files

# default feature set used in parser
FEATS=["s1c","s1v","s1p","s1A","s1a","s1B","s1b",
       "s1rL", # "s1rc","s1rv","s1rp",
       "s0c","s0v","s0p","s0A","s0B","s0a","s0b","s0d",
       "s0rL", # "s0rc","s0rv","s0rp",
       "n0lL", # "n0lc","n0lv","n0lp",
       "n0c","n0v","n0p","n0A","n0a",
       "n1c","n1v","n1p",
       ]

# first step to close buffer related features
FEATS1=["s1c","s1v","s1p","s1A","s1a","s1B","s1b",
       "s1rL", # "s1rc","s1rv","s1rp",
       "s0c","s0v","s0p","s0A","s0B","s0a","s0b","s0d",
       "s0rL" # "s0rc","s0rv","s0rp",
       #"n0lL", # "n0lc","n0lv","n0lp",
       #"n0c","n0v","n0p","n0A","n0a",
       #"n1c","n1v","n1p",
       ]

# second step to close stack related features
FEATS2=[#"s1c","s1v","s1p","s1A","s1a","s1B","s1b",
       #"s1rL", # "s1rc","s1rv","s1rp",
       #"s0c","s0v","s0p","s0A","s0B","s0a","s0b","s0d",
       #"s0rL", # "s0rc","s0rv","s0rp",
       "n0lL", # "n0lc","n0lv","n0lp",
       "n0c","n0v","n0p","n0A","n0a",
       "n1c","n1v","n1p",
       ]

using JLD, Knet
include("../src/header.jl")

const language_model  = "/scratch/users/okirnap/ud-treebanks-v2.2/chmodel_converted/english_chmodel.jld"
const small_data_file = "/scratch/users/okirnap/ud-treebanks-v2.2/UD_English-LinES/en_lines-ud-train.conllu"
const small_dev_file  = "/scratch/users/okirnap/ud-treebanks-v2.2/UD_English-LinES/en_lines-ud-dev.conllu"

const odict = Dict{Symbol, Any}();
odict[:embed] = (128, 32, 16); odict[:optimization] = Adam; odict[:hidden] = [2048]; odict[:feats]=FEATS1; odict[:arctype] = ArcHybridR1;
odict[:stembed]=odict[:bufembed]=128+950; odict[:sthidden]=odict[:bufhidden]=odict[:acthidden]=256; odict[:actembed]=32; odict[:posembed]=128;


text:
const ED_DEFAULT_SCALE  = 0.2
const ED_DEFAULT_SIGMA  = 2
const ED_DEFAULT_ITER   = 1
const ED_DEFAULT_BORDER = false
const ED_DEFAULT_NORM   = true

@doc """
    ElasticDistortion <: Augmentor.ImageOperation

Description
--------------

Distorts the given image using a randomly (uniform) generated
vector field of the given grid size. This field will be stretched
over the given image when applied, which in turn will morph the
original image into a new image using a linear interpolation of
both the image and the vector field.

In contrast to [`RandomDistortion`], the resulting vector
field is also smoothed using a Gaussian filter with of parameter
`sigma`. This will result in a less chaotic vector field and thus
resemble a more natural distortion.

Usage
--------------

    ElasticDistortion(gridheight, gridwidth, scale, sigma, [iter=$ED_DEFAULT_ITER], [border=$ED_DEFAULT_BORDER], [norm=$ED_DEFAULT_NORM])

    ElasticDistortion(gridheight, gridwidth, scale; [sigma=$ED_DEFAULT_SIGMA], [iter=$ED_DEFAULT_ITER], [border=$ED_DEFAULT_BORDER], [norm=$ED_DEFAULT_NORM])

    ElasticDistortion(gridheight, [gridwidth]; [scale=$ED_DEFAULT_SCALE], [sigma=$ED_DEFAULT_SIGMA], [iter=$ED_DEFAULT_ITER], [border=$ED_DEFAULT_BORDER], [norm=$ED_DEFAULT_NORM])

Arguments
--------------

- **`gridheight`** : The grid height of the displacement vector
    field. This effectively specifies the number of vertices
    along the Y dimension used as landmarks, where all the
    positions between the grid points are interpolated.

- **`gridwidth`** : The grid width of the displacement vector
    field. This effectively specifies the number of vertices
    along the Y dimension used as landmarks, where all the
    positions between the grid points are interpolated.

- **`scale`** : Optional. The scaling factor applied to all
    displacement vectors in the field. This effectively defines
    the "strength" of the deformation. There is no theoretical
    upper limit to this factor, but a value somewhere between
    `0.01` and `1.0` seem to be the most reasonable choices.
    Default to `$ED_DEFAULT_SCALE`.

- **`sigma`** : Optional. Sigma parameter of the Gaussian filter.
    This parameter effectively controls the strength of the
    smoothing. Defaults to `$ED_DEFAULT_SIGMA`.

- **`iter`** : Optional. The number of times the smoothing
    operation is applied to the displacement vector field. This
    is especially useful if `border = false` because the border
    will be reset to zero after each pass. Thus the displacement
    is a little less aggressive towards the borders of the image
    than it is towards its center. Defaults to
    `$ED_DEFAULT_ITER`.

- **`border`** : Optional. Specifies if the borders should be
    distorted as well. If `false`, the borders of the image will
    be preserved. This effectively pins the outermost vertices on
    their original position and the operation thus only distorts
    the inner content of the image. Defaults to
    `$ED_DEFAULT_BORDER`.

- **`norm`** : Optional. If `true`, the displacement vectors of
    the field will be normalized by the norm of the field. This
    will have the effect that the `scale` factor should be more
    or less independent of the grid size. Defaults to
    `$ED_DEFAULT_NORM`.

See also
--------------

[`augment`](@ref)

Examples
--------------

```julia
using Augmentor
img = testpattern()

# distort with pinned borders
augment(img, ElasticDistortion(15, 15; scale = 0.1))

# distort everything more smoothly.
augment(img, ElasticDistortion(10, 10; sigma = 4, iter=3, border=true))
```
""" ->
struct ElasticDistortion <: ImageOperation
    gridheight::Int
    gridwidth::Int
    scale::Float64
    sigma::Float64
    iterations::Int
    border::Bool
    normalize::Bool

    function ElasticDistortion(
            gridheight::Integer, gridwidth::Integer,
            scale, sigma, iterations = ED_DEFAULT_ITER,
            border::Bool = ED_DEFAULT_BORDER, normalize::Bool = ED_DEFAULT_NORM)
        (gridwidth > 2 && gridheight > 2) || throw(ArgumentError("gridwidth and gridheight need to be greater than 2"))
        sigma > 0 || throw(ArgumentError("sigma needs to be greater than 0"))
        iterations > 0 || throw(ArgumentError("iterations needs to be greater than 0"))
        new(Int(gridheight), Int(gridwidth), Float64(scale), Float64(sigma), Int(iterations), border, normalize)
    end
end

function ElasticDistortion(
        gridheight::Integer, gridwidth::Integer, scale;
        sigma = ED_DEFAULT_SIGMA, iter = ED_DEFAULT_ITER,
        border::Bool = ED_DEFAULT_BORDER, norm::Bool = ED_DEFAULT_NORM)
    ElasticDistortion(gridheight, gridwidth, scale, sigma, iter, border, norm)::ElasticDistortion
end

function ElasticDistortion(
        gridheight::Integer, gridwidth::Integer = gridheight;
        scale = ED_DEFAULT_SCALE,
        sigma = ED_DEFAULT_SIGMA, iter = ED_DEFAULT_ITER,
        border::Bool = ED_DEFAULT_BORDER, norm::Bool = ED_DEFAULT_NORM)
    ElasticDistortion(gridheight, gridwidth, scale, sigma, iter, border, norm)::ElasticDistortion
end

@inline supports_eager(::Type{<:ElasticDistortion}) = false
@inline supports_lazy(::Type{<:ElasticDistortion})  = true

function randparam(op::ElasticDistortion, img)
    gaussian_field(op.gridheight, op.gridwidth, op.scale, op.sigma, op.iterations, op.border, op.normalize)
end

function applylazy(op::ElasticDistortion, img::AbstractArray, field)
    DistortedView(img, field)
end

function showconstruction(io::IO, op::ElasticDistortion)
    str_size   = join(map(string,(op.gridheight,op.gridwidth)), ", ")
    str_scale  = op.scale      == ED_DEFAULT_SCALE  ? "" : "scale=$(op.scale)"
    str_sigma  = op.sigma      == ED_DEFAULT_SIGMA  ? "" : "sigma=$(op.sigma)"
    str_iter   = op.iterations == ED_DEFAULT_ITER   ? "" : "iter=$(op.iterations)"
    str_border = op.border     == ED_DEFAULT_BORDER ? "" : "border=$(op.border)"
    str_norm   = op.normalize  == ED_DEFAULT_NORM   ? "" : "norm=$(op.normalize)"
    str_all = filter(str->str!="", [str_size, str_scale, str_sigma, str_iter, str_border, str_norm])
    print(io, typeof(op).name.name, '(', join(str_all, ", "), ')')
end

function Base.show(io::IO, op::ElasticDistortion)
    if get(io, :compact, false)
        str_size   = join(map(string,(op.gridheight,op.gridwidth)), "×")
        str_sigma  = op.sigma      == ED_DEFAULT_SIGMA ? "" : "(σ=$(op.sigma)) "
        str_iter   = op.iterations == 1     ? "" : "$(op.iterations)-times "
        str_norm   = op.normalize  == true  ? "and normalized " : ""
        str_border = op.border     == false ? " with pinned border" : ""
        print(io, "Distort using a $(str_iter)smoothed $(str_norm)$(str_size) grid$(str_border)")
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end


text:
export extend, NfToNfMor, automorphisms

struct NfMorSet{T}
  field::T
end

function show(io::IO, S::NfMorSet{T}) where {T}
  print(io, "Set of automorphisms of ", S.field)
end

mutable struct NfToNfMor <: Map{AnticNumberField, AnticNumberField, HeckeMap, NfToNfMor}
  header::MapHeader{AnticNumberField, AnticNumberField}
  prim_img::nf_elem
  prim_preimg::nf_elem

  function NfToNfMor()
    z = new()
    z.header = MapHeader()
    return r
  end

  function NfToNfMor(K::AnticNumberField, L::AnticNumberField, y::nf_elem, isomorphism::Bool = false)
    z = new()
    z.prim_img = y

    function _image(x::nf_elem)
      g = parent(K.pol)(x)
      return evaluate(g, y)
    end

    if !isomorphism
      z.header = MapHeader(K, L, _image)
      return z
    end

    M = zero_matrix(FlintQQ, degree(L), degree(L))
    b = basis(K)
    for i = 1:degree(L)
      c = _image(b[i])
      for j = 1:degree(L)
        M[j, i] = coeff(c, j - 1)
      end
    end
    t = zero_matrix(FlintQQ, degree(L), 1)
    if degree(L) == 1
      t[1, 1] = coeff(gen(L), 0)
    else
      t[2, 1] = fmpq(1) # coefficient vector of gen(L)
    end

    s = solve(M, t)
    z.prim_preimg = K(parent(K.pol)([ s[i, 1] for i = 1:degree(K) ]))

    function _preimage(x::nf_elem)
      g = parent(L.pol)(x)
      return evaluate(g, z.prim_preimg)
    end

    z.header = MapHeader(K, L, _image, _preimage)
    return z
  end

  function NfToNfMor(K::AnticNumberField, L::AnticNumberField, y::nf_elem, y_inv::nf_elem)
    z = new()
    z.prim_img = y
    z.prim_preimg = y_inv

    function _image(x::nf_elem)
      g = parent(K.pol)(x)
      return evaluate(g, y)
    end

    function _preimage(x::nf_elem)
      g = parent(L.pol)(x)
      return evaluate(g, y_inv)
    end

    z.header = MapHeader(K, L, _image, _preimage)
    return z
  end
end

parent(f::NfToNfMor) = NfMorSet(domain(f))

################################################################################
#
#  NfToNfRelMor
#
################################################################################

mutable struct NfToNfRel <: Map{AnticNumberField, NfRel{nf_elem}, HeckeMap, NfToNfRel}
  header::MapHeader{AnticNumberField, NfRel{nf_elem}}

  function NfToNfRel(L::AnticNumberField, K::NfRel{nf_elem}, a::nf_elem, b::nf_elem, c::NfRelElem{nf_elem})
    # let K/k, k absolute number field
    # k -> L, gen(k) -> a
    # K -> L, gen(K) -> b
    # L -> K, gen(L) -> c

    k = K.base_ring
    Ly, y = PolynomialRing(L, cached = false)
    R = parent(k.pol)
    S = parent(L.pol)

    function image(x::nf_elem)
      # x is an element of L
      f = S(x)
      res = evaluate(f, c)
      return res
    end

    function preimage(x::NfRelElem{nf_elem})
      # x is an element of K
      f = data(x)
      # First evaluate the coefficients of f at a to get a polynomial over L
      # Then evaluate at b
      r = Vector{nf_elem}(undef, degree(f) + 1)
      for  i = 0:degree(f)
        r[i+1] = evaluate(R(coeff(f, i)), a)
      end
      return evaluate(Ly(r), b)
    end

    z = new()
    z.header = MapHeader(L, K, image, preimage)
    return z
  end
end

function show(io::IO, h::NfToNfRel)
  println(io, "Morphism between ", domain(h), "\nand ", codomain(h))
end

################################################################################
#
#  Generic groups to set of homomorphisms
#
################################################################################

mutable struct GrpGenToNfMorSet{T} <: Map{GrpGen, NfMorSet{T}, HeckeMap, GrpGenToNfMorSet{T}}
  header::MapHeader{GrpGen, NfMorSet{T}}

  function GrpGenToNfMorSet(G::GrpGen, S::NfMorSet{T}) where {T}
    z = new{T}()
    z.header = MapHeader(G, S)
    return z
  end
end

function GrpGenToNfMorSet(G::GrpGen, K::AnticNumberField)
  return GrpGenToNfMorSet(G, NfMorSet(K))
end

function image(f::GrpGenToNfMorSet, g::GrpGenElem)
  K = codomain(f).field
  return automorphisms(K)[g[]]
end

function (f::GrpGenToNfMorSet)(g::GrpGenElem)
  return image(f, g)
end

function preimage(f::GrpGenToNfMorSet, a::NfToNfMor)
  K = codomain(f).field
  aut = automorphisms(K, copy = false)
  for i in 1:length(aut)
    if a == aut[i]
      return domain(f)[i]
    end
  end
  error("something wrong")
end

function inv(f::NfToNfMor)
  if degree(domain(f)) != degree(codomain(f))
    error("The map is not invertible")
  end
  if isdefined(f, :prim_preimg)
    return hom(codomain(f), domain(f), f.prim_preimg, check = false)
  end
  L = codomain(f)
  K = domain(f)
  M = zero_matrix(FlintQQ, degree(L), degree(L))
  a = f.prim_img
  el = one(L)
  M[1, 1] = 1
  for i = 2:degree(L)
    mul!(el, el, a)
    for j = 1:degree(L)
      M[j, i] = coeff(el, j - 1)
    end
  end
  t = zero_matrix(FlintQQ, degree(L), 1)
  if degree(L) == 1
    t[1, 1] = coeff(gen(L), 0)
  else
    t[2, 1] = fmpq(1) # coefficient vector of gen(L)
  end

  s = solve(M, t)
  img = K(parent(K.pol)([ s[i, 1] for i = 1:degree(K) ]))
  return hom(L, K, img, check = false)
end

function _compute_preimg(m::NfToNfMor)
  # build the matrix for the basis change
  K = domain(m)
  L = codomain(m)
  M = zero_matrix(FlintQQ, degree(L), degree(L))
  b = basis(K)
  for i = 1:degree(L)
    c = m(b[i])
    for j = 1:degree(L)
      M[j, i] = coeff(c, j - 1)
    end
  end
  t = zero_matrix(FlintQQ, degree(L), 1)
  t[2, 1] = fmpq(1) # coefficient vector of gen(L)
  s = solve(M, t)
  m.prim_preimg = K(parent(K.pol)([ s[i, 1] for i = 1:degree(K) ]))
  return m.prim_preimg
end

function Base.:(==)(f::NfToNfMor, g::NfToNfMor)
  if (domain(f) != domain(g)) || (codomain(f) != codomain(g))
    return false
  end

  return f.prim_img == g.prim_img
end

#_D = Dict()

function evaluate(f::fmpq_poly, a::nf_elem)
  R = parent(a)
  if iszero(f)
    return zero(R)
  end
  l = length(f) - 1
  s = R(coeff(f, l))
  for i in l-1:-1:0
    #s = s*a + R(coeff(f, i))
    mul!(s, s, a)
    # TODO (easy): Once fmpq_poly_add_fmpq is improved in flint, remove the R(..)
    add!(s, s, R(coeff(f, i)))
  end
  return s
end

function *(f::NfToNfMor, g::NfToNfMor)
#  global _D
#  _s = Base.stacktrace()[2:3]
#  if !(_s in keys(_D))
#    _D[_s] = true
#    println("Called ...")
#    Base.show_backtrace(stdout, Base.stacktrace()[2:3])
#  end
  codomain(f) == domain(g) || throw("Maps not compatible")

  a = gen(domain(f))
  y = g(f(a))

  return hom(domain(f), codomain(g), y, check = false)
end

function ^(f::NfToNfMor, b::Int)
  K = domain(f)
  @assert K == codomain(f)
  d = degree(K)
  b = mod(b, d)
  if b == 0
    return NfToNfMor(K, K, gen(K))
  elseif b == 1
    return f
  else
    bit = ~((~UInt(0)) >> 1)
    while (UInt(bit) & b) == 0
      bit >>= 1
    end
    z = f
    bit >>= 1
    while bit != 0
      z = z * z
      if (UInt(bit) & b) != 0
        z = z * f
      end
      bit >>= 1
    end
    return z
  end
end

Base.copy(f::NfToNfMor) = f

Base.hash(f::NfToNfMor, h::UInt) = Base.hash(f.prim_img, h)

function show(io::IO, h::NfToNfMor)
  if domain(h) == codomain(h)
    println(io, "Automorphism of ", domain(h))
  else
    println(io, "Injection of ", domain(h), " into ", codomain(h))
  end
  println(io, "defined by ", gen(domain(h)), " -> ", h.prim_img)
end

function _automorphisms(K::AnticNumberField)
  if degree(K) == 1
    auts = NfToNfMor[hom(K, K, one(K))]
  else
    f = K.pol
    Kt, t = PolynomialRing(K, "t", cached = false)
    f1 = change_ring(f, Kt)
    divpol = Kt(nf_elem[-gen(K), K(1)])
    f1 = divexact(f1, divpol)
    lr = roots(f1, max_roots = div(degree(K), 2))
    Aut1 = Vector{NfToNfMor}(undef, length(lr)+1)
    for i = 1:length(lr)
      Aut1[i] = hom(K, K, lr[i], check = false)
    end
    Aut1[end] = id_hom(K)
    auts = closure(Aut1, degree(K))
  end
  return auts
end

function automorphisms(K::AnticNumberField; copy::Bool = true)
  try
    Aut = _get_automorphisms_nf(K)::Vector{NfToNfMor}
    if copy
      return Base.copy(Aut)
    else
      return Aut
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
  end
  auts = _automorphisms(K)
  _set_automorphisms_nf(K, auts)
  if copy
    return Base.copy(auts)
  else
    return auts
  end
end

function automorphism_group(K::AnticNumberField)
  aut = automorphisms(K)
  n = degree(K)
  #First, find a good prime
  p = 11
  d = numerator(discriminant(K.pol))
  while mod(d, p) == 0
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  pols = gfp_poly[Rx(g.prim_img) for g in aut]
  #if g in aut
  #  push!(pols, Rx(g.prim_img))
  #end
  Dcreation = Vector{Tuple{gfp_poly, Int}}(undef, length(pols))
  for i = 1:length(pols)
    Dcreation[i] = (pols[i], i)
  end
  D = Dict{gfp_poly, Int}(Dcreation)
  @assert length(D) == n
  mult_table = Array{Int, 2}(undef, n, n)
  #permutations = Array{Array{Int, 1},1}(undef, n)
  for s = 1:n
    #perm = Array{Int, 1}(undef, length(aut))
    for i = 1:length(aut)
      mult_table[s, i] = D[Hecke.compose_mod(pols[s], pols[i], fmod)]
    end
  end
  G = GrpGen(mult_table)
  return G, GrpGenToNfMorSet(G, K)
  #return _perm_to_gap_grp(permutations)
end

function hom(K::AnticNumberField, L::AnticNumberField, a::nf_elem; check::Bool = true, compute_inverse::Bool = false)
 if check
   if !iszero(evaluate(K.pol, a))
     error("The data does not define a homomorphism")
   end
 end
 return NfToNfMor(K, L, a, compute_inverse)
end

function hom(K::AnticNumberField, L::AnticNumberField, a::nf_elem, a_inv::nf_elem; check::Bool = true)
 if check
   if !iszero(evaluate(K.pol, a))
     error("The data does not define a homomorphism")
   end
   if !iszero(evaluate(L.pol, a_inv))
     error("The data does not define a homomorphism")
   end
 end
 return NfToNfMor(K, L, a, a_inv)
end

id_hom(K::AnticNumberField) = hom(K, K, gen(K), check = false)

morphism_type(::Type{AnticNumberField}) = NfToNfMor

isinjective(m::NfToNfMor) = true

issurjective(m::NfToNfMor) = (degree(domain(m)) == degree(codomain(m)))

isbijective(m::NfToNfMor) = issurjective(m)

###############################################################################
#
#  NfToNfMor closure
#
###############################################################################

function closure(S::Vector{NfToNfMor}, final_order::Int = -1)

  K = domain(S[1])
  d = numerator(discriminant(K.pol))
  p = 11
  while mod(d, p) == 0
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)

  t = length(S)
  order = 1
  elements = [id_hom(K)]
  pols = gfp_poly[x]
  gpol = Rx(S[1].prim_img)
  if gpol != x
    push!(pols, gpol)
    push!(elements, S[1])
    order += 1

    gpol = compose_mod(gpol, pols[2], fmod)

    while gpol != x
      order = order +1
      push!(elements, S[1]*elements[end])
      push!(pols, gpol)
      gpol = compose_mod(gpol, pols[2], fmod)
    end
  end
  if order == final_order
    return elements
  end

  for i in 2:t
    if !(S[i] in elements)
      pi = Rx(S[i].prim_img)
      previous_order = order
      order = order + 1
      push!(elements, S[i])
      push!(pols, Rx(S[i].prim_img))
      for j in 2:previous_order
        order = order + 1
        push!(pols, compose_mod(pols[j], pi, fmod))
        push!(elements, elements[j]*S[i])
      end
      if order == final_order
        return elements
      end
      rep_pos = previous_order + 1
      while rep_pos <= order
        for k in 1:i
          s = S[k]
          po = Rx(s.prim_img)
          att = compose_mod(pols[rep_pos], po, fmod)
          if !(att in pols)
            elt = elements[rep_pos]*s
            order = order + 1
            push!(elements, elt)
            push!(pols, att)
            for j in 2:previous_order
              order = order + 1
              push!(pols, compose_mod(pols[j], att, fmod))
              push!(elements, elements[j] *elt)
            end
            if order == final_order
              return elements
            end
          end
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  return elements
end

function small_generating_set(Aut::Array{NfToNfMor, 1})
  K = domain(Aut[1])
  a = gen(K)
  Identity = Aut[1]
  for i in 1:length(Aut)
    Au = Aut[i]
    if Au.prim_img == a
      Identity = Aut[i]
      break
    end
  end
  return small_generating_set(Aut, *, Identity)
end

function generic_group(G::Vector{NfToNfMor}, ::typeof(*))
  K = domain(G[1])
  n = length(G)
  #First, find a good prime
  p = 11
  d = numerator(discriminant(K.pol))
  while mod(d, p) == 0
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  pols = gfp_poly[Rx(g.prim_img) for g in G]
  Dcreation = Vector{Tuple{gfp_poly, Int}}(undef, length(pols))
  for i = 1:length(pols)
    Dcreation[i] = (pols[i], i)
  end
  D = Dict{gfp_poly, Int}(Dcreation)
  @assert length(D) == degree(K)
  permutations = Array{Array{Int, 1},1}(undef, n)

  m_table = Array{Int, 2}(undef, n, n)

  for s = 1:n
    for i = 1:n
      m_table[s, i] =  D[Hecke.compose_mod(pols[s], pols[i], fmod)]
    end
  end

  Gen = GrpGen(m_table)
  GentoG = Dict{GrpGenElem, eltype(G)}(Gen[i] => G[i] for i in 1:length(G))
  GtoGen = Dict{eltype(G), GrpGenElem}(G[i] => Gen[i] for i in 1:length(G))
  return Gen, GtoGen, GentoG
end

################################################################################
#
#  Induced image
#
################################################################################

# For Carlo:
(f::NfToNfMor)(x::NfOrdIdl) = induce_image(f, x)

function induce_image(f::NfToNfMor, x::NfOrdIdl)
  domain(f) !== codomain(f) && throw(error("Map must be an automorphism"))
  OK = order(x)
  K = nf(OK)
  I = ideal(OK)
  if isdefined(x, :gen_two)
    I.gen_two = OK(f(K(x.gen_two)))
  end
  if isdefined(x, :princ_gen)
    I.princ_gen = OK(f(K(x.princ_gen)))
  end
  for i in [:gen_one, :is_prime, :gens_normal, :gens_weakly_normal, :is_principal,
          :iszero, :minimum, :norm, :splitting_type]
    if isdefined(x, i)
      setfield!(I, i, getfield(x, i))
    end
  end
  if !has_2_elem(I)
    #I need to translate the basis matrix
    bb = Vector{NfOrdElem}(undef, degree(K))
    B = basis(x, copy = false)
    for i = 1:length(bb)
      bb[i] = OK(f(K(B[i])))
    end
    I.basis = bb
    M = zero_matrix(FlintZZ, degree(K), degree(K))
    for i = 1:degree(K)
      el = coordinates(I.basis[i])
      for j = 1:degree(K)
        M[i, j] = el[j]
      end
    end
    I.basis_matrix = M
  end
  return I
end

################################################################################
#
#  Maps to algebras
#
################################################################################

# Embedding of a number field into an algebra over Q.
mutable struct NfAbsToAbsAlgAssMor{S} <: Map{AnticNumberField, S, HeckeMap, NfAbsToAbsAlgAssMor}
  header::MapHeader{AnticNumberField, S}
  mat::fmpq_mat
  t::fmpq_mat

  function NfAbsToAbsAlgAssMor{S}(K::AnticNumberField, A::S, M::fmpq_mat) where { S <: AbsAlgAss{fmpq} }
    z = new{S}()
    z.mat = M
    z.t = zero_matrix(FlintQQ, 1, degree(K))

    function _image(x::nf_elem)
      for i = 1:degree(K)
        z.t[1, i] = coeff(x, i - 1)
      end
      s = z.t*z.mat
      return A([ s[1, i] for i = 1:dim(A) ])
    end

    z.header = MapHeader{AnticNumberField, S}(K, A, _image)
    return z
  end
end

function NfAbsToAbsAlgAssMor(K::AnticNumberField, A::S, M::fmpq_mat) where { S <: AbsAlgAss{fmpq} }
  return NfAbsToAbsAlgAssMor{S}(K, A, M)
end

function haspreimage(m::NfAbsToAbsAlgAssMor, a::AbsAlgAssElem)
  A = parent(a)
  t = matrix(FlintQQ, 1, dim(A), coeffs(a))
  b, p = can_solve(m.mat, t, side = :left)
  if b
    return true, domain(m)([ p[1, i] for i = 1:nrows(m.mat) ])
  else
    return false, zero(domain(m))
  end
end


text:
##############################################################################
#
# REFERENCES: Forbes et al. "Statistical Distributions"
#
##############################################################################

immutable BetaPrime <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    function BetaPrime(a::Real, b::Real)
        (a > zero(a) && b > zero(b)) || error("Alpha and beta must be positive")
        new(float64(a), float64(b))
    end
end

BetaPrime() = BetaPrime(1.0, 1.0)

@continuous_distr_support BetaPrime 0.0 Inf

function mean(d::BetaPrime)
    d.beta > 1.0 ? d.alpha / (d.beta - 1.0) : NaN
end

mode(d::BetaPrime) = d.alpha > 1.0 ? (d.alpha - 1.0) / (d.beta + 1.0) : 0.0

function pdf(d::BetaPrime, x::Real)
    α, β = d.alpha, d.beta
    (x^(α - 1.0) * (1.0 + x)^(-(α + β))) / beta(α, β)
end

cdf(d::BetaPrime, q::Real) = cdf(Beta(d.alpha, d.beta), q / (one(q) + q))
ccdf(d::BetaPrime, q::Real) = ccdf(Beta(d.alpha, d.beta), q / (one(q) + q))
logcdf(d::BetaPrime, q::Real) = logcdf(Beta(d.alpha, d.beta), q / (one(q) + q))
logccdf(d::BetaPrime, q::Real) = logccdf(Beta(d.alpha, d.beta), q / (one(q) + q))

quantile(d::BetaPrime, p::Real) = (x = quantile(Beta(d.alpha,d.beta),p); x / (1.0-x))
cquantile(d::BetaPrime, p::Real) = (x = cquantile(Beta(d.alpha,d.beta),p); x / (1.0-x))
invlogcdf(d::BetaPrime, p::Real) = (x = invlogcdf(Beta(d.alpha,d.beta),p); x / (1.0-x))
invlogccdf(d::BetaPrime, p::Real) = (x = invlogccdf(Beta(d.alpha,d.beta),p); x / (1.0-x))

function rand(d::BetaPrime)
    x = rand(Gamma(d.alpha, 1))
    y = rand(Gamma(d.beta, 1))
    x/y
end

function skewness(d::BetaPrime)
    α, β = d.alpha, d.beta
    β > 3.0 ? (2.0 * (2.0 * α + β - 1))/(β - 3.0) * sqrt((β - 2.0)/(α * (α + β - 1.0))) : NaN
end

function var(d::BetaPrime)
    α, β = d.alpha, d.beta
    β > 2.0 ? (α * (α + β - 1.0)) / ((β - 2.0) * (β - 1.0)^2) : NaN
end


text:
using Yao, YaoExtensions
using Yao.ConstGate
using LinearAlgebra, Test, Random
using Yao.AD: Rotor, generator

@testset "Diff Block" begin
    reg = rand_state(4)
    block = put(4, 2=>rot(X, 0.3))
    df = Diff(block)
    @test nqubits(df) == 4

    df2 = Diff(rot(CNOT, 0.3))
    @test nqubits(df2) == 2

    reg = rand_state(4)
    df2 = Diff(rot(CNOT, 0.3))
    @test nqubits(df2) == 2

    @test df2' isa Diff
    @test mat(df2) == mat(df2')'

    circuit = chain(put(4, 1=>Rx(0.1)), control(4, 2, 1=>Ry(0.3)))
    c2 = circuit |> markdiff
    @test c2[1].content isa Diff
    @test !(c2[2] isa Diff)
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
"Package demontrating how to make your code into a Julia pacakge"
module ExamplePkg

using Distributions

export pdf_normal

"""
   pdf_normal(x, mu, sigma)

Return pdf of a normal distribution with mean mu and standard deviation sigma at point x.
"""
function pdf_normal(x::Number, mu::Number, sigma::Number)
  dist = Normal(mu,sigma)
  pdf(dist,x)
end

# The code below is optional.  It serves to demonstrate how to make a program use the module code above.

"Parse command-line arguements for demo program"
function parse_args_for_demo(argv::Vector{String})
  if length(argv) != 3
     @error("Usage:  demo x mu sigma")
     exit()
  end
  args = try
     (; zip((:x,:mu,:sigma), parse.(Float64, argv) )... )
  catch ex
     @error("Error parsing arguments.\nx mu and sigma must all be floating point numbers.")
     exit()
  end
  @assert(!isnan(args[:x]) && !isinf(args[:x]))
  @assert(!isnan(args[:mu]) && !isinf(args[:mu]))
  @assert(!isnan(args[:sigma]) && !isinf(args[:sigma]))
  @assert(args[:sigma]>0)
  return args
end

"""
Demo program to evaluate normal PDF from command line.
Prepare executable using [PackageCompiler.jl](https://julialang.github.io/PackageCompiler.jl/stable/apps.html)
"""
Base.@ccallable function julia_main()::Cint
  x, mu, sigma = parse_args_for_demo(ARGS)
  result = pdf_normal(x,mu,sigma)
  println("prob(x | mu, sigma) = ", result)
  return 0 # if things finished successfully
end

println("# Precompiling common signatures...")
# Instruct Julia to pre-compile a functions with concrete types
precompile(pdf_normal,(Float64,Float64,Float64,))

end # module


text:
# Hack around Flux to allow multiple inputs to `Chain`
Flux.applychain(::Tuple{}, x1, xs...) = (x1, xs...)
Flux.applychain(::Tuple{}, x::Tuple) = x
Flux.applychain(fs::Tuple, x1, xs...) = Flux.applychain(Base.tail(fs), first(fs)(x1, xs...))
Flux.applychain(fs::Tuple, x::Tuple) = Flux.applychain(Base.tail(fs), first(fs)(x...))

(c::Chain)(x1, xs...) = Flux.applychain(Tuple(c.layers), x1, xs...)
(c::Chain)(x::Tuple) = c(x...)

# Some wrappers around common functions for prettier printing
"""
    ReshapeLayer(dims)

Reshapes the passed arrays to have a size of `(dims..., :)`
"""
struct ReshapeLayer{D}
    dims::D
end

(r::ReshapeLayer)(x) = reshape(x, r.dims..., :)
(r::ReshapeLayer)(xs...) = r.(xs)

"""
    FlattenLayer()

Flattens the passed arrays to 2D matrices.
"""
struct FlattenLayer end

(f::FlattenLayer)(x::AbstractArray{T,N}) where {T,N} = reshape(x, :, size(x, N))
(f::FlattenLayer)(xs...) = f.(xs)

"""
    SelectDim(dim, i)

See the documentation for `selectdim` for more information.
"""
struct SelectDim{I}
    dim::Int
    i::I
end

(s::SelectDim)(x) = selectdim(x, s.dim, s.i)
(s::SelectDim)(xs...) = s.(xs)

"""
    NoOpLayer()

As the name suggests. Does nothing. Better to use this than
`identity` since it takes in multiple inputs
"""
struct NoOpLayer end

(noop::NoOpLayer)(x...) = x
(noop::NoOpLayer)(x) = x


text:
struct Surface{T} <: LavaAbstraction
    handle::Vk.SurfaceKHR
    target::T
end

vk_handle_type(::Type{Surface{T}}) where {T} = (Vk.SurfaceKHR,T)

function Surface(win::XCBWindow)
    handle = create(Surface{XCBWindow}, Vk.XcbSurfaceCreateInfo(win.conn.h, win.id))
    Surface(handle, win)
end

struct Swapchain{T}
    handle::Vk.SwapchainKHR
    surface::Surface{T}
end

vk_handle_type(::Type{<:Swapchain}) = Vk.SwapchainKHR


text:
# pseudo dataset fot testing

function getdataset(dims, num)
  X = [rand(dims...) for _ in 1:num]
  Y = rand(num)
  collect(zip(X, Y))
end


text:
function Relationship(e::EzXML.Node) :: Relationship
    @assert EzXML.nodename(e) == "Relationship" "Unexpected XMLElement: $(EzXML.nodename(e)). Expected: \"Relationship\"."

    return Relationship(
        e["Id"],
        e["Type"],
        e["Target"]
    )
end

function parse_relationship_target(prefix::String, target::String) :: String
    @assert !isempty(prefix) && !isempty(target)

    if target[1] == '/'
        @assert sizeof(target) > 1 "Incomplete target path $target."
        return target[2:end]
    else
        return prefix * '/' * target
    end
end

function get_relationship_target_by_id(prefix::String, wb::Workbook, Id::String) :: String
    for r in wb.relationships
        if Id == r.Id
            return parse_relationship_target(prefix, r.Target)
        end
    end
    error("Relationship Id=$(Id) not found")
end

function get_relationship_target_by_type(prefix::String, wb::Workbook, _type_::String) :: String
    for r in wb.relationships
        if _type_ == r.Type
            return parse_relationship_target(prefix, r.Target)
        end
    end
    error("Relationship Type=$(_type_) not found")
end

function has_relationship_by_type(wb::Workbook, _type_::String) :: Bool
    for r in wb.relationships
        if _type_ == r.Type
            return true
        end
    end
    false
end

function get_package_relationship_root(xf::XLSXFile) :: EzXML.Node
    xroot = xmlroot(xf, "_rels/.rels")
    @assert EzXML.nodename(xroot) == "Relationships" "Malformed XLSX file $(xf.filepath). _rels/.rels root node name should be `Relationships`. Found $(EzXML.nodename(xroot))."
    @assert EzXML.namespaces(xroot) == Pair{String,String}[""=>"http://schemas.openxmlformats.org/package/2006/relationships"]
    return xroot
end

function get_workbook_relationship_root(xf::XLSXFile) :: EzXML.Node
    xroot = xmlroot(xf, "xl/_rels/workbook.xml.rels")
    @assert EzXML.nodename(xroot) == "Relationships" "Malformed XLSX file $(xf.filepath). xl/_rels/workbook.xml.rels root node name should be `Relationships`. Found $(EzXML.nodename(xroot))."
    @assert EzXML.namespaces(xroot) == Pair{String,String}[""=>"http://schemas.openxmlformats.org/package/2006/relationships"]
    return xroot
end

"""
Adds new relationship. Returns new generated rId.
"""
function add_relationship!(wb::Workbook, target::String, _type::String) :: String
    xf = get_xlsxfile(wb)
    @assert is_writable(xf) "XLSXFile instance is not writable."
    local rId :: String

    let
        got_unique_id = false
        id = 1

        while !got_unique_id
            got_unique_id = true
            rId = string("rId", id)
            for r in wb.relationships
                if r.Id == rId
                    got_unique_id = false
                    id += 1
                    break
                end
            end
        end
    end

    # adds to relationship vector
    new_relationship = Relationship(rId, _type, target)
    push!(wb.relationships, new_relationship)

    # adds to XML tree
    xroot = get_workbook_relationship_root(xf)
    el = EzXML.addelement!(xroot, "Relationship")
    el["Id"] = rId
    el["Target"] = target
    el["Type"] = _type

    return rId
end


text:
using Distributions
using LightBayes

function simulate(r::Float64, n::Int)
    # Simulate Δb for H0 & H1
    #
    # generate two clusters of samples, respectively centered
    # at (r, 0) and (-r, 0), where each cluster contains n
    # samples.
    #
    # This function computes the values of Δb, respectively
    # for the following two hypotheses:
    #
    # - H0: all samples are generated by a single component
    # - H1: each cluster of samples are generated by a different
    #       component
    #
    # Prior model: θ ~ N(0, σ0^2) with σ0 = 10.0
    # Likelihood:  x ~ N(θ, σ^2) with σ = 1.0
    #

    # data
    n = 1000
    X1 = [r, 0.0] .+ rand(2, 1000)
    X2 = [-r, 0.0] .+ rand(2, 1000)

    X = hcat(X1, X2)
    @assert size(X) == (2, 2n)

    # prior
    β0 = zeros(2)
    σ0 = 10.0
    κ0 = inv(σ0^2)
    pri = MvNormalCanon(β0, κ0)

    # likelihood model
    σ = 1.0
    g = IsoGaussModel(σ)

    # H0: single component
    pp0 = posterior(pri, suffstats(g, X))
    Δb0 = logpar(pp0) - logpar(pri)

    # H1: two components:
    pp1 = posterior(pri, suffstats(g, X1))
    pp2 = posterior(pri, suffstats(g, X2))
    Δb1 = (logpar(pp1) - logpar(pri)) + (logpar(pp2) - logpar(pri))

    return (Δb0, Δb1)
end

# main script

n = 1000

for r in [0.0, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0]
    Δb0, Δb1 = simulate(r, n)
    @printf("With r = %5.1f:  H0.Δb = %8.2f,  H1.Δb = %8.2f\n", r, Δb0, Δb1)
end