text:
using ReachabilityAnalysis, Symbolics, Plots

LazySets.set_ztol(Float64, 1e-15)

function multistable_oscillator(; X0 = Interval(0.0, 0.05),
                                  V₊ = +13.5, V₋ = -13.5,
                                  R = 20.E3, C = 5.5556E-8,
                                  R1 = 20.E3, R2 = 20.E3)

    @variables x
    τ = 1/(R*C)
    α = R2/(R1+R2)

    A = -τ
    b = (τ/α) * V₊
    I₊ = HalfSpace(x <= α*V₊)
    m1 = @system(x' = Ax + b, x ∈ I₊)

    b = (τ/α) * V₋
    I₋ = HalfSpace(x >= α*V₋)
    m2 = @system(x' = Ax + b, x ∈ I₋)

    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    g1 = Hyperplane(x == α*V₊)
    r1 = ConstrainedIdentityMap(1, g1)

    add_transition!(automaton, 2, 1, 2)
    g2 = Hyperplane(x == α*V₋)
    r2 = ConstrainedIdentityMap(1, g2)

    modes = [m1, m2]
    resetmaps = [r1, r2]
    H = HybridSystem(automaton, modes, resetmaps)

    # initial condition in mode 1
    X0e = [(1, X0)]
    return IVP(H, X0e)
end

prob = multistable_oscillator()

sol = solve(prob, T=100e-4, alg=INT(δ=1.E-6), fixpoint_check=false);

plot(sol, vars=(0, 1), xlab="t", ylab="v-")

tspan.(sol)

location.(sol)

plot(sol[1][end-10:end], vars=(0, 1), xlab="t", ylab="v-")
plot!(x -> 6.75, xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red)

Xc = cluster(sol[1], [318, 319, 320], BoxClustering(1))

plot(sol[1][end-10:end], vars=(0, 1))
plot!(sol[2][1:10], vars=(0, 1))
plot!(x -> 6.75, xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red)
plot!(Xc[1], vars=(0, 1), c=:grey)

sol = solve(prob, T=100e-4, alg=INT(δ=1.E-6), fixpoint_check=true);

plot(sol, vars=(0, 1), xlab="t", ylab="v-")

tspan(sol)


text:
using .RCall
import .RCall: rcopy, rcopytype, RClass
function rcopy(::Type{PointPattern}, s::Ptr{VecSxp})
    x = rcopy(Vector, s[:x])
    y = rcopy(Vector, s[:y])
    window = rcopy(Dict{Symbol, Any}, s[:window])

    PointPattern(collect(zip(x,y)), (x=window[:xrange], y=window[:yrange]))
end

rcopytype(::Type{RClass{:ppp}}, ::Ptr{VecSxp}) = PointPattern

import .RCall: sexp, protect, unprotect, setclass!, RClass, sexpclass
function sexp(::Type{RClass{:ppp}}, pp::PointPattern)
    # win = protect(sexp(Dict(:type => "regtangle",
    #     :xrange => [pp.window.x...],
    #     :yrange => [pp.window.y...])))
    # setclass!(win, sexp("owin"))
    # r = protect(sexp(Dict(
    #     :n => length(pp.data),
    #     :x => getx.(pp.data),
    #     :y => gety.(pp.data),
    #     :window => win)))
    # setclass!(r, sexp("ppp"))
    # unprotect(2)
    r = R"spatstat.geom::ppp"(getx.(pp.data), gety.(pp.data),
        window=R"spatstat.geom::owin"([pp.window.x...], [pp.window.y...]))
    r
end

sexpclass(f::PointPattern) = RClass{:ppp}

function rcopy(::Type{NamedTuple{(:x, :y), Tuple{Float64, Float64}}}, s::Ptr{VecSxp})
    if !rcopy(R"spatstat.geom::is.rectangle"(s))
        @warn "Non-rectangle windows not supported. Using the bounding box as the window."
    end

    window = rcopy(Dict{Symbol, Any}, s)

    (x=window[:xrange], y=window[:yrange])
end

rcopytype(::Type{RClass{:owin}}, ::Ptr{VecSxp}) = NamedTuple{(:x, :y), Tuple{Float64, Float64}}

ext: jl


text:
using DataFrames

"""
	detect_ignore(df, missing_threshold=0.5, detect_id=true)

Return list of column names, where the column has too many missing, or the name is "Id" and the like.
Set `missing_thresold` to change the rate of missing allowed.
"""
function detect_ignore(df; missing_threshold=0.5, detect_id=true)
	filter(df |> names) do name
		if lowercase(name) == "id" && detect_id
			true
		else
			col = df[!, name]
			n = length(col)
			nm = count(ismissing, n)
			(n / nm) ≥ 0.5
		end
	end
end

"""
	detect(df, T::Type; exclude=String[], use_eltype=true, inv=false)

Return list of column names, where the `eltype(df[!, name]) <: T`.
Use `inv=true` for the inverse.
Use `use_eltype=false` for the whole column type instead of eltype.
Use `exclude` to exclude specific names.
"""
function detect(df, T::Type; exclude=String[], use_eltype=true, inv=false)
	return filter(names(df)) do name
		col = df[!, name]
		if name in exclude
			return false
		else
			S = use_eltype ? eltype(col) : typeof(col)
			cmp = S <: T
			return inv ? !cmp : cmp
		end
	end
end

function detect_categorical(df; exclude=String[])
	n1 = detect(df, Union{Number, Missing}; exclude=exclude, inv=true)
	n2 = detect(df, CategoricalArray; exclude=exclude, use_eltype=false)
	n1 ∪ n2
end

function detect_numeric(df; exclude=String[])
	detect(df, Union{Number, Missing}; exclude=exclude)
end

struct DummyCoding
	columns
	new_columns
	transformation
end

"""
	DummyCoding(df, categoricals)

Return a DummyCoding transformation, which transform categorical columns to dummy coding fake columns and remove the original one.
"""
function DummyCoding(df::AbstractDataFrame, categoricals)
	transforms = mapreduce(vcat, categoricals) do name
		col = df[!, name]
		vals = unique(col)
		@. name => ByRow(isequal(vals)) .=> Symbol(name * "_", vals)
	end
	new_columns = last.(last.(transforms))
	DummyCoding(categoricals, new_columns, transforms)
end

function transform_(df, t::DummyCoding)
	df = transform(df, t.transformation)
	select(df, Not(t.columns))
end

function HierarchicalCoding(H::Dict)
	function maybe_lowercase(x)
		ismissing(x) ? missing : lowercase(x)
	end
	map(collect(H)) do (k, v)
		lv = maybe_lowercase.(v)
		function f(x)
			y = @chain begin
				maybe_lowercase.(x)
				indexin(_, lv)
			end
			@assert !(nothing in y) "Missing level code for $k: $(unique(x))"
			convert(Vector{Int}, y)
		end
		k => f => k
	end
end


text:
#
println("Ba) Tests of the dipole, em multipole and momentum-transfer amplitudes.")
#
grid = Radial.Grid("grid: exponential")
wa = Atomic.Computation("xx",  Nuclear.Model(26.); grid=grid, properties=JAC.AtomicLevelProperty[],
                        configs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")] )

wxa  = perform(wa; output=true)
wma  = wxa["multiplet:"]

flow = 6;    fup = 8;   ilow = 1;   iup = 3
println("\n\nDipole amplitudes:\n")
for  finalLevel in wma.levels
    for  initialLevel in wma.levels
        if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
            ##x println("Levels with indices f = $(finalLevel.index) and i = $(initialLevel.index)")
            MultipoleMoment.dipoleAmplitude(finalLevel, initialLevel, grid; display=true)
        end
    end
end

println("\n\nElectromagnetic multipole-transition (absorption) amplitudes from PhotoEmission.amplitude(..):\n")
for  finalLevel in wma.levels
    for  initialLevel in wma.levels
        if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
            PhotoEmission.amplitude("absorption", E1, Basics.Coulomb,  1.0, finalLevel, initialLevel, grid; display=true)
            PhotoEmission.amplitude("absorption", M1, Basics.Magnetic, 1.0, finalLevel, initialLevel, grid; display=true)
        end
    end
end

println("\n\nElectromagnetic multipole-transition (absorption) amplitudes from MultipoleMoment.transitionAmplitude(..):\n")
for  finalLevel in wma.levels
    for  initialLevel in wma.levels
        if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
            MultipoleMoment.transitionAmplitude(E1, Basics.Velocity, 1.0, finalLevel, initialLevel, grid; display=true)
            MultipoleMoment.transitionAmplitude(M1, Basics.Magnetic, 1.0, finalLevel, initialLevel, grid; display=true)
        end
    end
end

println("\n\nElectromagnetic multipole-moment amplitudes from MultipoleMoment.amplitude(..):\n")
for  finalLevel in wma.levels
    for  initialLevel in wma.levels
        if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
            MultipoleMoment.amplitude(E1, Basics.Velocity, 1.0, finalLevel, initialLevel, grid; display=true)
            MultipoleMoment.amplitude(M1, Basics.Magnetic, 1.0, finalLevel, initialLevel, grid; display=true)
        end
    end
end

println("\n\nMomentum transfer amplitudes:\n")
for  finalLevel in wma.levels
    for  initialLevel in wma.levels
        if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
            FormFactor.amplitude(1.0, finalLevel, initialLevel, grid; display=true)
        end
    end
end


text:
# Analysis for the Dobble Drilling variant. (Drilling = Triplet)
# In this variant, there are 9 open dobble cards and a player has to find three
# cards which share a symbol. They then take the three cards and three new cards
# are added to the game.

# https://en.wikipedia.org/wiki/Dobble

# This is using a projective space over Z/7Z, cards correspond to points and
# symbols correspond to lines.

using JuMP
import GLPK

# Question: How many cards can we choose such that no Drilling appears?

maximumOnALine = 2

model = Model(GLPK.Optimizer)

# Regular points
@variable(model, p[1:7, 1:7], Bin)

# Infinite points, on slope (1, d)
@variable(model, pi[1:7], Bin)

# Infinite infinite on slope (0, 1)
@variable(model, pii, Bin)

# Lines are Constraints
# There is the infinite line which does not touch any regular point
@constraint(model, sum(pi) + pii <= maximumOnALine)

# Then there are lines with a delta of (0, 1) through p[x, 1] and pii
for x in 1:7
    @constraint(model, sum(p[x, 1:7]) + pii <= maximumOnALine)
end

# # All other lines have a slope of (1, d), start at p[1, y] and go through pi[d]

for y in 1:7, d in 0:6
    line = pi[(d + 6) % 7 + 1]
    for t in 0:6
        line = line + p[1 + t, (y + t * d - 1) % 7 + 1]
    end
    @constraint(model, line <= maximumOnALine)
end

# Now the objective is to place as many cards as possible

@objective(model, Max, sum(p) + sum(pi) + pii)

@show "Starting to optimize"

optimize!(model)

@show termination_status(model)
@show objective_value(model)

# Running this shows, that there is a possibility to choose 8 cards such that
# no Drilling occurs, but that this is also optimal. This means 9 cards always
# have a Drilling.

max_stars_repo_name: kreibaum/dobble

max_issues_repo_name: kreibaum/dobble

max_forks_repo_name: kreibaum/dobble


text:
#=
test_xmlparser:
- Julia version: 1.1.0
- Author: bramb
- Date: 2019-02-20
=#
using Test

@testset "xmlparser" begin
    include("../src/xmlparser.jl")

    xml = "<xml>Mondays are <del>well good</del><add>def bad</add>!</xml>"
    root = xp_parse(xml)
    # dump(root)

    all_nodes = create_an_array_of_the_xml_nodes(root)
    blocks = convert_to_xml_blocks(all_nodes)
    println(blocks)

    # We mappen de XML blocks naar een combo van tag en nul of meer text nodes.

    tuples_of_tag_and_text_nodes::Array{NamedTuple} = []
    for block in blocks
        named_tuple = transform_block_into_text_nodes_named_tuple(block)
        push!(tuples_of_tag_and_text_nodes, named_tuple)
    end

    println(tuples_of_tag_and_text_nodes)

#     textNodes = create_nodes_for_each_xml_block(blocks)
#     println(textNodes)

    # de partities doen we later wel.
    #partitions = partition_block_into_groups(blocks)
    #println(partitions)
end


text:
function cite_package(name, bib; cite_commands)
    isempty(bib) && return ""
    if length(bib) == 1
        key = only(keys(bib))
        cmd = get(cite_commands, name, DEFAULT_CITE)

        cmd * "[" * name * "]" * "{" * key * "}"
    else
        ks = keys(bib)
        # determine common cite command
        cmd = get(cite_commands, first(ks), DEFAULT_CITE)
        for key in ks
            cmdᵢ = get(cite_commands, name, DEFAULT_CITE)
            # all the keys must have the same command
            if cmdᵢ ≠ cmd
                # fall back to DEFAULT_CITE
                cmd = DEFAULT_CITE
            end
        end

        cmd * "[" * name * "]" * "{" * join(ks, ',') * "}"
    end
end

function start_sentence(n, cite_commands=Dict{String,String}())
    bib = get_julia_bib()
    julia = cite_package("Julia v$VERSION", bib; cite_commands)
    if n < 3
        "This work was done in $julia and made use of the "
    else
        "This work was done in $julia and made use of the following packages: "
    end
end

function sentence_ending(n)
    if n == 1
        " package."
    elseif n == 2
        " packages."
    else
        "."
    end
end

function make_sentence(pkg_citations; cite_commands=Dict{String,String}(), jl=true, texttt=false)
    pkgs = keys(pkg_citations)
    n = length(pkgs)

    start = start_sentence(n, cite_commands)
    ending = sentence_ending(n)

    citations = String[]
    for pkg in pkgs
        pkg_name = jl ? pkg * ".jl" : pkg
        pkg_name = texttt ? add_texttt(pkg_name) : pkg_name
        c = pkg_name * cite_package(pkg, pkg_citations[pkg]; cite_commands)
        push!(citations, c)
    end

    middle = n == 2 ? citations[1] * " and " * citations[2] : join(citations, ", ", " and ")

    return start * middle * ending
end

"""
    get_tool_citation(io::IO=stdout; jl = true, texttt = false, copy = true, cite_commands=Dict{String,String}(), filename="julia_citations.bib"; badge=false)

Print a sentence describing the packages used in the current environment.
If you only want to consider the direct dependencies, you can set `only_direct=true`.
The sentence is automatically copied to the clipboard(you can avoid this by using `copy=false`).
The package names have the .jl ending by default. You can ommit it with `jl=false`.
Package names can be wrapped in `texttt` by setting `texttt=true` and you can also customize
the cite command used for each package by using `cite_commands=Dict("PackageName"=>"custom_cite")`.
The filename of the .bib file can be passed via the `filename` keyword.
Use `badge = true` to get the citations from packages without
a `Citation.bib` file, but with a DOI badge.
"""
function get_tool_citation(io::IO=stdout;
    jl = true,
    texttt = false,
    copy = true,
    only_direct = false,
    cite_commands=Dict{String,String}(),
    filename="julia_citations.bib",
    badge=false)

    pkg_citations = collect_citations(only_direct, badge=badge)

    cite_sentence = make_sentence(pkg_citations; cite_commands, jl, texttt)
    println(io, cite_sentence)

    try
        if copy
            clipboard(cite_sentence)
            @info "The above sentence was copied to your clipboard."
        end
    catch e
        @error e
    end

    julia_bib = get_julia_bib()
    all_citations = merge!(Dict("Julia"=>julia_bib), pkg_citations)
    export_citations(filename, all_citations)
end


text:
module operators_sparse

export SparseOperator, diagonaloperator

import Base: ==, *, /, +, -, Broadcast
import ..operators
import ..operators_dense: DataOperatorStyle, DenseOperatorStyle
import SparseArrays: sparse

using ..bases, ..states, ..operators, ..operators_dense, ..sparsematrix
using SparseArrays, LinearAlgebra

"""
    SparseOperator(b1[, b2, data])

Sparse array implementation of Operator.

The matrix is stored as the julia built-in type `SparseMatrixCSC`
in the `data` field.
"""
mutable struct SparseOperator{BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} <: DataOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    data::T
    function SparseOperator{BL,BR,T}(b1::Basis, b2::Basis, data::T) where {BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}}
        if length(b1) != size(data, 1) || length(b2) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(b1, b2, data)
    end
end

SparseOperator(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} = SparseOperator{BL,BR,T}(b1, b2, data)
SparseOperator{BL,BR}(b1::BL, b2::BR, data::T) where {BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} = SparseOperator{BL,BR,T}(b1, b2, data)
SparseOperator(b1::Basis, b2::Basis, data) = SparseOperator(b1, b2, convert(SparseMatrixCSC{ComplexF64,Int}, data))
SparseOperator(b::Basis, data::SparseMatrixCSC{ComplexF64, Int}) = SparseOperator(b, b, data)
SparseOperator(b::Basis, data::Matrix{ComplexF64}) = SparseOperator(b, sparse(data))
SparseOperator(op::DenseOperator) = SparseOperator(op.basis_l, op.basis_r, sparse(op.data))

SparseOperator(b1::Basis, b2::Basis) = SparseOperator(b1, b2, spzeros(ComplexF64, length(b1), length(b2)))
SparseOperator{BL,BR}(b1::BL, b2::BR) where {BL<:Basis,BR<:Basis} = SparseOperator{BL,BR}(b1, b2, spzeros(ComplexF64, length(b1), length(b2)))
SparseOperator(b::Basis) = SparseOperator(b, b)

Base.copy(x::SparseOperator) = SparseOperator(x.basis_l, x.basis_r, copy(x.data))
operators.dense(a::SparseOperator) = DenseOperator(a.basis_l, a.basis_r, Matrix(a.data))

"""
    sparse(op::AbstractOperator)

Convert an arbitrary operator into a [`SparseOperator`](@ref).
"""
sparse(a::AbstractOperator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) instead."))
sparse(a::SparseOperator) = copy(a)
sparse(a::DenseOperator) = SparseOperator(a.basis_l, a.basis_r, sparse(a.data))

==(x::SparseOperator, y::SparseOperator) = false
==(x::T, y::T) where T<:SparseOperator = (x.data == y.data)

# Arithmetic operations
+(a::SparseOperator, b::SparseOperator) = throw(bases.IncompatibleBases())
+(a::T, b::T) where T<:SparseOperator = T(a.basis_l, a.basis_r, a.data+b.data)
+(a::SparseOperator, b::DenseOperator) = throw(bases.IncompatibleBases())
+(a::SparseOperator{B1,B2}, b::DenseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data+b.data)
+(a::DenseOperator, b::SparseOperator) = throw(bases.IncompatibleBases())
+(a::DenseOperator{B1,B2}, b::SparseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data+b.data)

-(a::T) where T<:SparseOperator = T(a.basis_l, a.basis_r, -a.data)
-(a::SparseOperator, b::SparseOperator) = throw(bases.IncompatibleBases())
-(a::T, b::T) where T<:SparseOperator = SparseOperator(a.basis_l, a.basis_r, a.data-b.data)
-(a::SparseOperator, b::DenseOperator) = throw(bases.IncompatibleBases())
-(a::SparseOperator{B1,B2}, b::DenseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data-b.data)
-(a::DenseOperator, b::SparseOperator) = throw(bases.IncompatibleBases())
-(a::DenseOperator{B1,B2}, b::SparseOperator{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperator{B1,B2}(a.basis_l, a.basis_r, a.data-b.data)

*(a::SparseOperator{B1,B2}, b::SparseOperator{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = SparseOperator{B1,B3}(a.basis_l, b.basis_r, a.data*b.data)
*(a::SparseOperator, b::SparseOperator) = throw(bases.IncompatibleBases())
*(a::T, b::Number) where T<:SparseOperator = T(a.basis_l, a.basis_r, complex(b)*a.data)
*(a::Number, b::T) where T<:SparseOperator = T(b.basis_l, b.basis_r, complex(a)*b.data)

/(a::T, b::Number) where T<:SparseOperator = T(a.basis_l, a.basis_r, a.data/complex(b))

operators.dagger(x::SparseOperator) = SparseOperator(x.basis_r, x.basis_l, x.data')
operators.ishermitian(A::SparseOperator) = false
operators.ishermitian(A::SparseOperator{B,B}) where B<:Basis = ishermitian(A.data)

operators.tensor(a::SparseOperator, b::SparseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))
operators.tensor(a::DenseOperator, b::SparseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))
operators.tensor(a::SparseOperator, b::DenseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data))

operators.tr(op::SparseOperator{B,B}) where B<:Basis = tr(op.data)
operators.tr(op::SparseOperator) = throw(bases.IncompatibleBases())

operators.conj(op::T) where T<:SparseOperator = T(op.basis_l, op.basis_r, conj(op.data))
operators.conj!(op::SparseOperator) = conj!(op.data)

operators.transpose(op::SparseOperator{BL,BR,T}) where {BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} = SparseOperator{BR,BL,T}(op.basis_r, op.basis_l, T(transpose(op.data)))

function operators.ptrace(op::SparseOperator, indices::Vector{Int})
    operators.check_ptrace_arguments(op, indices)
    shape = [op.basis_l.shape; op.basis_r.shape]
    data = sparsematrix.ptrace(op.data, shape, indices)
    b_l = ptrace(op.basis_l, indices)
    b_r = ptrace(op.basis_r, indices)
    SparseOperator(b_l, b_r, data)
end

function operators.expect(op::SparseOperator{B,B}, state::Ket{B}) where B<:Basis
    state.data' * op.data * state.data
end

function operators.expect(op::SparseOperator{B1,B2}, state::DenseOperator{B2,B2}) where {B1<:Basis,B2<:Basis}
    result = ComplexF64(0.)
    @inbounds for colindex = 1:op.data.n
        for i=op.data.colptr[colindex]:op.data.colptr[colindex+1]-1
            result += op.data.nzval[i]*state.data[colindex, op.data.rowval[i]]
        end
    end
    result
end

function operators.permutesystems(rho::SparseOperator{B1,B2}, perm::Vector{Int}) where {B1<:CompositeBasis,B2<:CompositeBasis}
    @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm)
    @assert isperm(perm)
    shape = [rho.basis_l.shape; rho.basis_r.shape]
    data = sparsematrix.permutedims(rho.data, shape, [perm; perm .+ length(perm)])
    SparseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data)
end

operators.identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:SparseOperator} = SparseOperator(b1, b2, sparse(ComplexF64(1)*I, length(b1), length(b2)))
operators.identityoperator(b1::Basis, b2::Basis) = identityoperator(SparseOperator, b1, b2)
operators.identityoperator(b::Basis) = identityoperator(b, b)

"""
    diagonaloperator(b::Basis)

Create a diagonal operator of type [`SparseOperator`](@ref).
"""
function diagonaloperator(b::Basis, diag::Vector{T}) where T <: Number
  @assert 1 <= length(diag) <= prod(b.shape)
  SparseOperator(b, sparse(Diagonal(convert(Vector{ComplexF64}, diag))))
end

# Fast in-place multiplication implementations
operators.gemm!(alpha, M::SparseOperator{B1,B2}, b::DenseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = sparsematrix.gemm!(convert(ComplexF64, alpha), M.data, b.data, convert(ComplexF64, beta), result.data)
operators.gemm!(alpha, a::DenseOperator{B1,B2}, M::SparseOperator{B2,B3}, beta, result::DenseOperator{B1,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis} = sparsematrix.gemm!(convert(ComplexF64, alpha), a.data, M.data, convert(ComplexF64, beta), result.data)
operators.gemv!(alpha, M::SparseOperator{B1,B2}, b::Ket{B2}, beta, result::Ket{B1}) where {B1<:Basis,B2<:Basis} = sparsematrix.gemv!(convert(ComplexF64, alpha), M.data, b.data, convert(ComplexF64, beta), result.data)
operators.gemv!(alpha, b::Bra{B1}, M::SparseOperator{B1,B2}, beta, result::Bra{B2}) where {B1<:Basis,B2<:Basis} = sparsematrix.gemv!(convert(ComplexF64, alpha), b.data, M.data, convert(ComplexF64, beta), result.data)

# Broadcasting
struct SparseOperatorStyle{BL<:Basis,BR<:Basis} <: DataOperatorStyle{BL,BR} end
Broadcast.BroadcastStyle(::Type{<:SparseOperator{BL,BR}}) where {BL<:Basis,BR<:Basis} = SparseOperatorStyle{BL,BR}()
Broadcast.BroadcastStyle(::DenseOperatorStyle{B1,B2}, ::SparseOperatorStyle{B1,B2}) where {B1<:Basis,B2<:Basis} = DenseOperatorStyle{B1,B2}()
Broadcast.BroadcastStyle(::DenseOperatorStyle{B1,B2}, ::SparseOperatorStyle{B3,B4}) where {B1<:Basis,B2<:Basis,B3<:Basis,B4<:Basis} = throw(bases.IncompatibleBases())
Broadcast.BroadcastStyle(::SparseOperatorStyle{B1,B2}, ::SparseOperatorStyle{B3,B4}) where {B1<:Basis,B2<:Basis,B3<:Basis,B4<:Basis} = throw(bases.IncompatibleBases())

@inline function Base.copy(bc::Broadcast.Broadcasted{Style,Axes,F,Args}) where {BL<:Basis,BR<:Basis,Style<:SparseOperatorStyle{BL,BR},Axes,F,Args<:Tuple}
    bcf = Broadcast.flatten(bc)
    args_ = Tuple(a.data for a=bcf.args)
    bl,br = states.find_basis(bcf.args)
    bc_ = Broadcast.Broadcasted(bcf.f, args_, axes(bcf))
    return SparseOperator{BL,BR}(bl, br, copy(bc_))
end

end # module


text:
using RCall # before using PhyloNetworks and PhyloPlots
using PhyloPlots
using Test

using DataFrames
using PhyloNetworks # to read topologies
using Colors # is REQUIREd by Gadfly

@testset "PhyloPlots Tests" begin
  # each file should be its own testset, to run them all even if one has a failure
  include("test_phylonetworkPlots.jl")
  include("test_plotGadfly.jl") # uses Colors
  include("test_plotRCall.jl")  # uses DataFrames to test annotations
  include("test_substitutionmodels.jl")
  include("test_rexport.jl")    # uses RCall, but NOT library(ape)
end


text:
# ======== Extension of diff rules to complex functions ==========#

# preempt other function diffrules. The two entries correspond to d/dz and d/dz*
DiffRules.@define_diffrule Base.abs2(z) = :(conj($z)), :($z)
DiffRules.@define_diffrule Base.sqrt(z) = :(inv(2 * sqrt($z))), :(0)
DiffRules.@define_diffrule Base.log(z) = :(inv($z)), :(0)
DiffRules.@define_diffrule Base.:^(z,p) = :($p * ($z^($p - 1))), :(0)

DiffRules.@define_diffrule Base.abs(z) = :(0.5*conj($z)*inv(abs($z))), :(0.5*$z*inv(abs($z)))

# ======== Train functions to deal with complex duals ==========#

# Extend functions of a single complex argument to accept dual argument
macro extend_unary_dual_to_complex(fcn)
    M = :Base
    f = :($fcn)
    dfz, dfzstar = DiffRules.diffrule(M,f,:v)
    fdef = quote
        function $M.$f(z::ComplexDual{T}) where {T}
            v = value(z)
            dvz, dvzstar = dz_partials(z)

            # chain rule
            dfdz =     $dfz*dvz     + $dfzstar*conj(dvzstar)
            dfdzstar = $dfz*dvzstar + $dfzstar*conj(dvz)

            return ComplexDual{T}($M.$f(v),dfdz,dfdzstar)
        end
    end
    return esc(fdef)
end

# Extend functions of a single complex argument and another (non-dual) argument to accept
# dual argument in place of first argument
macro extend_binary_dual_to_complex(fcn)
    M = :Base
    f = :($fcn)
    dfz, dfzstar = DiffRules.diffrule(M,f,:v,:p)
    defs = quote end
    for R in AMBIGUOUS_TYPES
        expr = quote
            function $M.$f(z::ComplexDual{T},p::$R) where {T}
                v = value(z)
                dvz, dvzstar = dz_partials(z)

                # chain rule
                dfdz =     $dfz*dvz     + $dfzstar*conj(dvzstar)
                dfdzstar = $dfz*dvzstar + $dfzstar*conj(dvz)

                return ComplexDual{T}($M.$f(v,p),dfdz,dfdzstar)

            end
        end
        append!(defs.args, expr.args)
    end
    return esc(defs)
end

@extend_unary_dual_to_complex abs2
@extend_unary_dual_to_complex sqrt
@extend_unary_dual_to_complex log
@extend_unary_dual_to_complex abs

@extend_binary_dual_to_complex ^


text:
# tests Cairo and Pango

using Test
using Gadfly
using Cairo

@testset "Test exporting plot figure to PDF" begin

pl = Gadfly.plot(y=randn(10), Geom.line)

pl |> PDF("/tmp/test.pdf")

end


text:
##############################################################################
##
## Generalized reduced rank tests using the singularvalue decomposition
##
## See also the Stata command ranktest
## RANKTEST: Stata module to test the rank of a matrix using the Kleibergen-Paap rk statistic
## Authors: Frank Kleibergen, Mark E Schaffer
## More precisely, it corresponds to the Stata command:  ranktest  (X) (Z), wald full

## For Kleibergen-Paap first-stage F-statistics
## See "Generalized reduced rank tests using the singular value decomposition"
## Journal of Econometrics, 2006, 133:1
## Frank Kleibergen and Richard Paap
## doi:10.1016/j.jeconom.2005.02.011
##############################################################################

function ranktest!(X::Matrix{Float64},
                    Z::Matrix{Float64},
                    Pi::Matrix{Float64},
                    vcov_method::CovarianceEstimator,
                    dof_small::Int,
                    dof_fes::Int)

    if (size(X, 2) == 0) | (size(Z, 2) == 0)
        return NaN
    end

    # Compute theta
    Fmatrix = cholesky!(Symmetric(Z' * Z), check = false).U
    Gmatrix = cholesky!(Symmetric(X' * X), check = false).U
    theta = Fmatrix * (Gmatrix' \ Pi')'

    # compute lambda
    svddecomposition = svd(theta, full = true)
    u = svddecomposition.U
    vt = svddecomposition.Vt

    k = size(X, 2)
    l = size(Z, 2)

    u_sub = u[k:l, k:l]
    vt_sub = vt[k,k]

    # see p.102 of KP
    # the operation u_sup \ sqrt(u_sup * u_sub') may fail. For now, we handle the case k = l but not the more general one
    a_qq = iszero(u_sub) ? u[1:l, k:l] : u[1:l, k:l]  * (u_sub \ sqrt(u_sub * u_sub'))
    b_qq = iszero(vt_sub) ? vt[1:k, k]' : sqrt(vt_sub * vt_sub') * (vt_sub' \ vt[1:k, k]')

    kronv = kron(b_qq, a_qq')
    lambda = kronv * vec(theta)

    # compute vhat
    if vcov_method isa Vcov.SimpleCovariance
        vlab = cholesky!(Hermitian((kronv * kronv') ./ size(X, 1)), check = false)
    else
        K = kron(Gmatrix, Fmatrix)'
        vcovmodel = Vcov.VcovData(Z, K, X, size(Z, 1) - dof_small - dof_fes)
        matrix_vcov2 = Vcov.S_hat(vcovmodel, vcov_method)
        vhat = K \ (K \ matrix_vcov2)'
        vlab = cholesky!(Hermitian(kronv * vhat * kronv'), check = false)
    end
    r_kp = lambda' * (vlab \ lambda)
    return r_kp[1]
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# bring vsel to a normalized shape, i.e.:
#
# AbstractVector{<:Integer} -> AbstractVector{<:Integer}
# AbstractVector{AbstractUnitRange} -> AbstractVector{AbstractUnitRange}
#
# Expr -> Vector{Union{Symbol,Expr})
#
# Symbol -> Vector{Union{Symbol,Expr})

function _normalize_vsel_unshaped(vsel::Union{AbstractVector{T}, AbstractVector{UnitRange{T}}}) where T <: Integer
    return vsel
end

function _normalize_vsel_unshaped(vsel::Integer) where T <: Integer
    return Int[vsel]
end

function _normalize_vsel_unshaped(vsel::AbstractVector)
    @argcheck all(broadcast(x -> x isa Union{T, AbstractVector{T}, AbstractVector{UnitRange{T}}} where T <: Integer, vsel)) throw(ArgumentError("Samples or distribution are unshaped, please use Integer, Vector of Integers, or UnitRange for indexing."))
    return vsel
end

function _normalize_vsel_shaped(vs::AbstractValueShape, vsel::AbstractVector{UnitRange})
    @argcheck length(vsel) == 1 throw(ArgumentError("Can't use multidimensional cartesian index for shaped samples or distribution, please use keys instead"))
    return [key for key in keys(vs)[vsel[1]]]
end

function _normalize_vsel_shaped(vs::AbstractValueShape, vsel::Integer)
    return [keys(vs)[vsel]]
end

function _normalize_vsel_shaped(vs::AbstractValueShape, vsel::AbstractVector{Integer})
    return [keys(vs)[i] for i in vsel]
end

function _normalize_vsel_shaped(vs::AbstractValueShape, vsel::Expr)
    if isexpr(vsel, :vect) || isexpr(vsel, :call)
        if @capture(vsel, [s_:e_]) || @capture(vsel, (s_:e_))
            return [keys(vs)[i] for i in s:e]

        elseif @capture(vsel, [dims__]) && all(broadcast(x -> x isa Integer, dims))
            return [keys(vs)[i] for i in dims]
        else
            vsel_norm = [_normalize_vsel_shaped(vs, el isa QuoteNode ? el.value : el) for el in vsel.args]
            return vcat(vsel_norm...) # for user convenience, does :([a[1,2], b, c]) -> [:(a[1,2]), :b, :c]
        end
    end
    return [vsel]
end

function _normalize_vsel_shaped(vs::AbstractValueShape, vsel::Symbol)
    return [vsel]
end

function _normalize_vsel_shaped(vs::AbstractValueShape, vsel::AbstractVector)
    vsel_norm = [_normalize_vsel_shaped(vs, el isa QuoteNode ? el.value : el) for el in vsel]
    return vcat(vsel_norm...)
end

# Retrieve flat indexes of the selected dims of samples
# If shaped, build new shape

function marg_idxs_unshaped(samples::DensitySampleVector,
                            vsel::AbstractVector
)
    flat_v = ValueShapes.flatview(samples.v)
    return (firstindex(flat_v) - 1) .+ vsel
end

function marg_idxs_shaped(samples::DensitySampleVector,
                         vsel::AbstractVector
)
    shapes = []
    sel_idxs = Int[]

    vs = varshape(samples)
    flat_idxs = axes(ValueShapes.flatview(unshaped.(samples).v), 1)
    shaped_idxs = vs(flat_idxs)

    for el in vsel

        # :a
        if el isa Symbol
            shape_acc = getproperty(vs, el)
            idx_acc = getproperty(shaped_idxs, el)

            shape = shape_acc.shape
            append!(sel_idxs, idx_acc isa Integer ? Int64[idx_acc] : vec(idx_acc))

        # :(a[1]) or :(a[1:4]) or :(a[1,2,3]) or :(a[1,2,:])
        elseif @capture(el, a_[args__])
            shape_acc = getproperty(vs, a)
            idx_acc = getproperty(shaped_idxs, a)

            if @capture(el, a_[s_:e_])
                shape = ArrayShape{Real}(length(s:e))
                append!(sel_idxs, vec(idx_acc[s:e]))
            elseif @capture(el, a_[dim_])
                shape = ScalarShape{Real}()
                append!(sel_idxs, idx_acc[dim])
            else
                slice_dims = Int[]

                for (i,arg) in enumerate(args)
                    if @capture(arg, (s_:e_))

                        args[i] = s:e
                        append!(slice_dims, length(s:e))
                    else

                        append!(slice_dims, 0)
                    end
                end

                sliced_dims = slice_dims .> 0

                whole_dims = args .== :(:)
                args[whole_dims] .= Colon()

                new_dims = shape_acc.shape.dims .* whole_dims .+ slice_dims

                new_axes = new_dims[whole_dims .| sliced_dims]

                shape = length(new_axes) > 0 ? ArrayShape{Real}(new_axes...) : ScalarShape{Real}()

                idxs = idx_acc[args...]
                append!(sel_idxs, ndims(idxs) > 0 ? vec(idxs) : idxs)

            end
        end
        el = encode_name(el)
        append!(shapes, (el = shape,))
    end
    vsel = encode_name.(vsel)
    new_shape = NamedTupleShape(NamedTuple(vsel .=> shapes))

    return sel_idxs, new_shape
end

# Encode dim names to preserve indexing information
#
# unicode encodings to preserve indexing information:
# :a⌞2⌟     = : a \llcorner 2 \lrcorner
# :a⌞2ː4⌟   = : a \llcorner 2 \lmrk 4 \lrcorner
# :a⌞1ˌ3ˌ5⌟ = : a \llcorner 1 \verti 3 \verti 5 \lrcorner

function encode_name(name::Expr)
    code = replace(string(name)," " => "", "[" => "⌞", ":" => "ː", "," => "ˌ", "]" => "⌟")
    return Symbol(code)
end

function encode_name(name::Symbol)
    return name
end
#=
"""
    bat_marginalize(
        samples::DensitySampleVector,
        vsel::Union{E, S, I, U, AbstractVector{Union{E, S, I, U}}} where E<:Expr where S<:Symbol where I<:Integer where U<:UnitRange
    )

Create a new `DensitySampleVector` containing the selected
dimensions `vsel` of the input `samples`.

Returns a `DensitySampleVector`. If the input is shaped, the new samples
will have the shape of the selected dimensions.

`samples` may be shaped or unshaped, vsel may be an `Int`, `UnitRange`,
`Symbol`, `Expression`, or a Vector of these types.

To select certain elements or whole slices of a multidimensional parameter `a` of
`samples` do:

`vsel` = :(a[1]) or `vsel` = :(a[1,2]) or `vsel` = :(a[1,:,:])

with the brackets containing the cartesian index of the selected values.

To preserve the indexing information of the selection, `vsel` is encoded using
unicode characters. The new names in the marginalized samples are then `symbol`s
of these codes:

:(a[2])     => :a⌞2⌟
:(a[2:4])   => :a⌞2ː4⌟
:(a[1,2,:]) => :a⌞1ˌ2ˌː⌟

If `samples` is shaped, and `vsel` is an `Int` `i` or a vector of `Int`s `[i,j,k]`,
the i-th, j-th, and k-th parameters of `samples` are selected automatically.

If `samples` is unshaped, `vsel` may only be an ``Int` `i` or a vector
of `Int`s `[i,j,k]`. Then the i-th, j-th, and k-th columns of `samples` are selected.

For convenience, `vsel` can be of the shape

:([a[1,2], a[5:8], c]) instead of [:(a[1,2]), :(a[5:8]), :c]

The logd of the new samples is set to NaN, the info and aux fields of `samples` are copied.

Example:

```julia
marg_shaped_samples = bat_marginalize(orig_shaped_samples, :(a[1,2,:]))

varshape(marg_samples) == NamedTupleShape{(:a⌞1ˌ2ˌː⌟,), Tuple{ValueAccessor{ArrayShape{Real, 1}}}}((a⌞1ˌ2ˌː⌟ = ValueAccessor{ArrayShape{Real, 1}}(ArrayShape{Real, 1}((3,)), 0, 3),), 3)
```

```julia
marg_unshaped_samples = bat_marginalize(orig_unshaped_smaples, [1,2,3])
'''
"""
=#
function bat_marginalize(samples::DensitySampleVector,
                         vsel
)
    shaped = isshaped(samples)
    vs = varshape(samples)
    vsel = vsel isa Tuple ? [vsel...] : vsel

    if shaped
       vsel = _normalize_vsel_shaped(vs, vsel)
    else
       vsel = _normalize_vsel_unshaped(vsel)
    end

    if shaped
        flat_idxs, new_shape = marg_idxs_shaped(samples, vsel)
    else
        flat_idxs = marg_idxs_unshaped(samples, vsel)
    end

    flat_orig_data = ValueShapes.flatview(unshaped.(samples).v)
    marg_data = flat_orig_data[flat_idxs, :]
    marg_data = ndims(marg_data) > 1 ? nestedview(marg_data) : marg_data
    logd = fill(NaN, length(ndims(marg_data) > 1 ? axes(marg_data, 2) : marg_data))
    info = samples.info
    aux = samples.aux

    marg_samples = shaped ? new_shape.(DensitySampleVector(marg_data, logd, info = info, aux = aux)) : DensitySampleVector(marg_data, logd, info = info, aux = aux)
    return marg_samples
end


text:
module matrices
using LinearAlgebra
using ..layers
export ScatterMatrix,matrix_layer,matrix_ref,matrix_tra,concatenate
struct ScatterMatrix
    S11::Array{Complex{Float64},2}
    S12::Array{Complex{Float64},2}
    S21::Array{Complex{Float64},2}
    S22::Array{Complex{Float64},2}
end

function matrix_layer(V,W,X,V0)
    A=(Matrix(W)\I+(Matrix(V)\I)*V0)#*.5
    B=(Matrix(W)\I-(Matrix(V)\I)*V0)#*.5
    Ai=I/A
    S11=S22=(A-X*B*Ai*X*B)\(X*B*Ai*X*A-B)
    S12=S21=(A-X*B*Ai*X*B)\X*(A-B*Ai*B)
    return ScatterMatrix(S11,S12,S21,S22)
end

function matrix_layer(l::Layer,V0)
    return matrix_layer(l.V,l.W,l.X,V0)
end

function matrix_ref(V,W,V0)
    A=W+V0\Matrix(V)
    B=W-V0\Matrix(V)
    Ai=I/A
    S11=-Ai*B
    S12=2*Ai
    S21=.5*(A-B*Ai*B)
    S22=B*Ai
    return ScatterMatrix(S11,S12,S21,S22)
end

function matrix_ref(l::Halfspace,V0)
    return matrix_ref(l.V,l.W,V0)
end

function matrix_tra(V,W,V0)
    A=W+V0\Matrix(V)
    B=W-V0\Matrix(V)
    Ai=I/A
    S11=B*Ai
    S12=.5*(A-B*Ai*B)
    S21=2*Ai
    S22=-Ai*B
    return ScatterMatrix(S11,S12,S21,S22)
end

function matrix_tra(l::Halfspace,V0)
    return matrix_tra(l.V,l.W,V0)
end

function concatenate(S11a,S12a,S21a,S22a,S11b,S12b,S21b,S22b)
    S11=S11a+(S12a/(I-S11b*S22a))*S11b*S21a
    S12=(S12a/(I-S11b*S22a))*S12b
    S21=(S21b/(I-S22a*S11b))*S21a
    S22=S22b+(S21b/(I-S22a*S11b))*S22a*S12b
    return S11,S12,S21,S22
end

function concatenate(S1,S2)
    S11,S12,S21,S22=concatenate(S1.S11,S1.S12,S1.S21,S1.S22,S2.S11,S2.S12,S2.S21,S2.S22)
    return ScatterMatrix(S11,S12,S21,S22)
end

function concatenate(Sin::Array{ScatterMatrix,1})
    Sout=Sin[1]
    for i=2:length(Sin)
        Sout=concatenate(Sout,Sin[i])
    end
    return Sout
end

end


text:
# An ordered set of dicts which we traverse to find the parameter value.
# Returns nothing if no param setting is found in any of the dicts.
type Parameters
  dicts  # First hash is searched first and then in order until the last

  Parameters(dicts...) = begin
    # Ensure we have at least one hash there if none are given
    if length(dicts) == 0
      dicts = [Dict{Any, Any}()]
    end

    new(dicts)
  end
end

function getindex(p::Parameters, key)
  stringkey = string(key)
  symbolkey = symbol(key)
  for d in p.dicts
    if haskey(d, symbolkey)
      return getindex(d, symbolkey)
    elseif haskey(d, stringkey)
      return getindex(d, stringkey)
    end
  end
  return nothing
end

function setindex!(p::Parameters, value, key)
  setindex!(first(p.dicts), value, key)
end

import Base.haskey, Base.get, Base.merge, Base.delete!

function haskey(p::Parameters, key)
  for d in p.dicts
    if haskey(d, key)
      return true
    end
  end
  return false
end

function get(p::Parameters, key, default = nothing)
  value = getindex(p, key)
  return (value == nothing ? default : value)
end

# In a merge the last parameter should be prioritized since this is the way
# the normal Julia merge of dicts works.
merge(p1::Union(Dict, Parameters), p2::Union(Dict, Parameters)) = Parameters(p2, p1)

function delete!(p::Parameters, key)
  for d in p.dicts
    delete!(d, key)
  end
end


text:
function Base.getproperty(m::dnn_Layer, s::Symbol)
    if s==:blobs
        return cpp_to_julia(jlopencv_Layer_get_blobs(m))
    end
    if s==:name
        return cpp_to_julia(jlopencv_Layer_get_name(m))
    end
    if s==:type
        return cpp_to_julia(jlopencv_Layer_get_type(m))
    end
    if s==:preferableTarget
        return cpp_to_julia(jlopencv_Layer_get_preferableTarget(m))
    end
    return Base.getfield(m, s)
end

function Base.setproperty!(m::dnn_Layer, s::Symbol, v)
    if s==:name
        jlopencv_Layer_set_name(m, julia_to_cpp(v))
    end
    if s==:type
        jlopencv_Layer_set_type(m, julia_to_cpp(v))
    end
    if s==:preferableTarget
        jlopencv_Layer_set_preferableTarget(m, julia_to_cpp(v))
    end
    return Base.setfield!(m, s, v)
end

function finalize(cobj::cv_Ptr{T}, inputs::Array{InputArray, 1}, outputs::Array{InputArray, 1}) where {T <: dnn_Layer}
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Layer_cv_dnn_Layer_finalize(julia_to_cpp(cobj),julia_to_cpp(inputs),julia_to_cpp(outputs)))
end
finalize(cobj::cv_Ptr{T}, inputs::Array{InputArray, 1}; outputs::Array{InputArray, 1} = (Array{InputArray, 1}())) where {T <: dnn_Layer} = finalize(cobj, inputs, outputs)

# function run(cobj::cv_Ptr{T}, inputs::Array{InputArray, 1}, internals::Array{InputArray, 1}, outputs::Array{InputArray, 1}) where {T <: dnn_Layer}
# 	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Layer_cv_dnn_Layer_run(julia_to_cpp(cobj),julia_to_cpp(inputs),julia_to_cpp(internals),julia_to_cpp(outputs)))
# end
# run(cobj::cv_Ptr{T}, inputs::Array{InputArray, 1}, internals::Array{InputArray, 1}; outputs::Array{InputArray, 1} = (Array{InputArray, 1}())) where {T <: dnn_Layer} = run(cobj, inputs, internals, outputs)

function outputNameToIndex(cobj::cv_Ptr{T}, outputName::String) where {T <: dnn_Layer}
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Layer_cv_dnn_Layer_outputNameToIndex(julia_to_cpp(cobj),julia_to_cpp(outputName)))
end
function Base.getproperty(m::dnn_Net, s::Symbol)
    return Base.getfield(m, s)
end
function Base.setproperty!(m::dnn_Net, s::Symbol, v)
    return Base.setfield!(m, s, v)
end

function empty(cobj::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_empty(julia_to_cpp(cobj)))
end

function dump(cobj::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_dump(julia_to_cpp(cobj)))
end

function dumpToFile(cobj::dnn_Net, path::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_dumpToFile(julia_to_cpp(cobj),julia_to_cpp(path)))
end

function getLayerId(cobj::dnn_Net, layer::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getLayerId(julia_to_cpp(cobj),julia_to_cpp(layer)))
end

function getLayerNames(cobj::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getLayerNames(julia_to_cpp(cobj)))
end

function getLayer(cobj::dnn_Net, layerId::dnn_LayerId)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getLayer(julia_to_cpp(cobj),julia_to_cpp(layerId)))
end

function connect(cobj::dnn_Net, outPin::String, inpPin::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_connect(julia_to_cpp(cobj),julia_to_cpp(outPin),julia_to_cpp(inpPin)))
end

function setInputsNames(cobj::dnn_Net, inputBlobNames::Array{String, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_setInputsNames(julia_to_cpp(cobj),julia_to_cpp(inputBlobNames)))
end

function setInputShape(cobj::dnn_Net, inputName::String, shape::Array{Int32, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_setInputShape(julia_to_cpp(cobj),julia_to_cpp(inputName),julia_to_cpp(shape)))
end

function forward(cobj::dnn_Net, outputName::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_forward(julia_to_cpp(cobj),julia_to_cpp(outputName)))
end
forward(cobj::dnn_Net; outputName::String = (String())) = forward(cobj, outputName)

function forward(cobj::dnn_Net, outputBlobs::Array{InputArray, 1}, outputName::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_forward(julia_to_cpp(cobj),julia_to_cpp(outputBlobs),julia_to_cpp(outputName)))
end

function forward(cobj::dnn_Net, outBlobNames::Array{String, 1}, outputBlobs::Array{InputArray, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_forward(julia_to_cpp(cobj),julia_to_cpp(outBlobNames),julia_to_cpp(outputBlobs)))
end
forward(cobj::dnn_Net, outBlobNames::Array{String, 1}; outputBlobs::Array{InputArray, 1} = (Array{InputArray, 1}())) = forward(cobj, outBlobNames, outputBlobs)

function forwardAsync(cobj::dnn_Net, outputName::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_forwardAsync(julia_to_cpp(cobj),julia_to_cpp(outputName)))
end
forwardAsync(cobj::dnn_Net; outputName::String = (String())) = forwardAsync(cobj, outputName)

function setHalideScheduler(cobj::dnn_Net, scheduler::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_setHalideScheduler(julia_to_cpp(cobj),julia_to_cpp(scheduler)))
end

function setPreferableBackend(cobj::dnn_Net, backendId::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_setPreferableBackend(julia_to_cpp(cobj),julia_to_cpp(backendId)))
end

function setPreferableTarget(cobj::dnn_Net, targetId::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_setPreferableTarget(julia_to_cpp(cobj),julia_to_cpp(targetId)))
end

function setInput(cobj::dnn_Net, blob::InputArray, name::String, scalefactor::Float64, mean::Scalar)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_setInput(julia_to_cpp(cobj),julia_to_cpp(blob),julia_to_cpp(name),julia_to_cpp(scalefactor),julia_to_cpp(mean)))
end
setInput(cobj::dnn_Net, blob::InputArray; name::String = (""), scalefactor::Float64 = Float64(1.0), mean::Scalar = ()) = setInput(cobj, blob, name, scalefactor, mean)

function setParam(cobj::dnn_Net, layer::dnn_LayerId, numParam::Int32, blob::InputArray)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_setParam(julia_to_cpp(cobj),julia_to_cpp(layer),julia_to_cpp(numParam),julia_to_cpp(blob)))
end

function getParam(cobj::dnn_Net, layer::dnn_LayerId, numParam::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getParam(julia_to_cpp(cobj),julia_to_cpp(layer),julia_to_cpp(numParam)))
end
getParam(cobj::dnn_Net, layer::dnn_LayerId; numParam::Int32 = Int32(0)) = getParam(cobj, layer, numParam)

function getUnconnectedOutLayers(cobj::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getUnconnectedOutLayers(julia_to_cpp(cobj)))
end

function getUnconnectedOutLayersNames(cobj::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getUnconnectedOutLayersNames(julia_to_cpp(cobj)))
end

function getLayersShapes(cobj::dnn_Net, netInputShapes::Array{Array{Int32, 1}, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getLayersShapes(julia_to_cpp(cobj),julia_to_cpp(netInputShapes)))
end

function getLayersShapes(cobj::dnn_Net, netInputShape::Array{Int32, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getLayersShapes(julia_to_cpp(cobj),julia_to_cpp(netInputShape)))
end

function getFLOPS(cobj::dnn_Net, netInputShapes::Array{Array{Int32, 1}, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getFLOPS(julia_to_cpp(cobj),julia_to_cpp(netInputShapes)))
end

function getFLOPS(cobj::dnn_Net, netInputShape::Array{Int32, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getFLOPS(julia_to_cpp(cobj),julia_to_cpp(netInputShape)))
end

function getFLOPS(cobj::dnn_Net, layerId::Int32, netInputShapes::Array{Array{Int32, 1}, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getFLOPS(julia_to_cpp(cobj),julia_to_cpp(layerId),julia_to_cpp(netInputShapes)))
end

function getFLOPS(cobj::dnn_Net, layerId::Int32, netInputShape::Array{Int32, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getFLOPS(julia_to_cpp(cobj),julia_to_cpp(layerId),julia_to_cpp(netInputShape)))
end

function getLayerTypes(cobj::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getLayerTypes(julia_to_cpp(cobj)))
end

function getLayersCount(cobj::dnn_Net, layerType::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getLayersCount(julia_to_cpp(cobj),julia_to_cpp(layerType)))
end

function getMemoryConsumption(cobj::dnn_Net, netInputShape::Array{Int32, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getMemoryConsumption(julia_to_cpp(cobj),julia_to_cpp(netInputShape)))
end

function getMemoryConsumption(cobj::dnn_Net, layerId::Int32, netInputShapes::Array{Array{Int32, 1}, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getMemoryConsumption(julia_to_cpp(cobj),julia_to_cpp(layerId),julia_to_cpp(netInputShapes)))
end

function getMemoryConsumption(cobj::dnn_Net, layerId::Int32, netInputShape::Array{Int32, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getMemoryConsumption(julia_to_cpp(cobj),julia_to_cpp(layerId),julia_to_cpp(netInputShape)))
end

function enableFusion(cobj::dnn_Net, fusion::Bool)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_enableFusion(julia_to_cpp(cobj),julia_to_cpp(fusion)))
end

function getPerfProfile(cobj::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_cv_dnn_Net_getPerfProfile(julia_to_cpp(cobj)))
end

function Base.getproperty(m::dnn_Model, s::Symbol)
    return Base.getfield(m, s)
end
function Base.setproperty!(m::dnn_Model, s::Symbol, v)
    return Base.setfield!(m, s, v)
end

function setInputSize(cobj::dnn_Model, size::Size)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_setInputSize(julia_to_cpp(cobj),julia_to_cpp(size)))
end

function setInputSize(cobj::dnn_Model, width::Int32, height::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_setInputSize(julia_to_cpp(cobj),julia_to_cpp(width),julia_to_cpp(height)))
end

function setInputMean(cobj::dnn_Model, mean::Scalar)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_setInputMean(julia_to_cpp(cobj),julia_to_cpp(mean)))
end

function setInputScale(cobj::dnn_Model, scale::Float64)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_setInputScale(julia_to_cpp(cobj),julia_to_cpp(scale)))
end

function setInputCrop(cobj::dnn_Model, crop::Bool)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_setInputCrop(julia_to_cpp(cobj),julia_to_cpp(crop)))
end

function setInputSwapRB(cobj::dnn_Model, swapRB::Bool)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_setInputSwapRB(julia_to_cpp(cobj),julia_to_cpp(swapRB)))
end

function setInputParams(cobj::dnn_Model, scale::Float64, size::Size, mean::Scalar, swapRB::Bool, crop::Bool)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_setInputParams(julia_to_cpp(cobj),julia_to_cpp(scale),julia_to_cpp(size),julia_to_cpp(mean),julia_to_cpp(swapRB),julia_to_cpp(crop)))
end
setInputParams(cobj::dnn_Model; scale::Float64 = Float64(1.0), size::Size = (SizeOP()), mean::Scalar = (), swapRB::Bool = (false), crop::Bool = (false)) = setInputParams(cobj, scale, size, mean, swapRB, crop)

function predict(cobj::dnn_Model, frame::InputArray, outs::Array{InputArray, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_predict(julia_to_cpp(cobj),julia_to_cpp(frame),julia_to_cpp(outs)))
end
predict(cobj::dnn_Model, frame::InputArray; outs::Array{InputArray, 1} = (Array{InputArray, 1}())) = predict(cobj, frame, outs)

function dnn_Model(model::String, config::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_Model(julia_to_cpp(model),julia_to_cpp(config)))
end
dnn_Model(model::String; config::String = ("")) = dnn_Model(model, config)

function dnn_Model(network::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Model_cv_dnn_Model_Model(julia_to_cpp(network)))
end
function Base.getproperty(m::dnn_ClassificationModel, s::Symbol)
    return Base.getfield(m, s)
end
function Base.setproperty!(m::dnn_ClassificationModel, s::Symbol, v)
    return Base.setfield!(m, s, v)
end

function classify(cobj::dnn_ClassificationModel, frame::InputArray)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_ClassificationModel_cv_dnn_ClassificationModel_classify(julia_to_cpp(cobj),julia_to_cpp(frame)))
end

function dnn_ClassificationModel(model::String, config::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_ClassificationModel_cv_dnn_ClassificationModel_ClassificationModel(julia_to_cpp(model),julia_to_cpp(config)))
end
dnn_ClassificationModel(model::String; config::String = ("")) = dnn_ClassificationModel(model, config)

function dnn_ClassificationModel(network::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_ClassificationModel_cv_dnn_ClassificationModel_ClassificationModel(julia_to_cpp(network)))
end
function Base.getproperty(m::dnn_KeypointsModel, s::Symbol)
    return Base.getfield(m, s)
end
function Base.setproperty!(m::dnn_KeypointsModel, s::Symbol, v)
    return Base.setfield!(m, s, v)
end

function estimate(cobj::dnn_KeypointsModel, frame::InputArray, thresh::Float32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_KeypointsModel_cv_dnn_KeypointsModel_estimate(julia_to_cpp(cobj),julia_to_cpp(frame),julia_to_cpp(thresh)))
end
estimate(cobj::dnn_KeypointsModel, frame::InputArray; thresh::Float32 = Float32(0.5)) = estimate(cobj, frame, thresh)

function dnn_KeypointsModel(model::String, config::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_KeypointsModel_cv_dnn_KeypointsModel_KeypointsModel(julia_to_cpp(model),julia_to_cpp(config)))
end
dnn_KeypointsModel(model::String; config::String = ("")) = dnn_KeypointsModel(model, config)

function dnn_KeypointsModel(network::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_KeypointsModel_cv_dnn_KeypointsModel_KeypointsModel(julia_to_cpp(network)))
end
function Base.getproperty(m::dnn_SegmentationModel, s::Symbol)
    return Base.getfield(m, s)
end
function Base.setproperty!(m::dnn_SegmentationModel, s::Symbol, v)
    return Base.setfield!(m, s, v)
end

function segment(cobj::dnn_SegmentationModel, frame::InputArray, mask::InputArray)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_SegmentationModel_cv_dnn_SegmentationModel_segment(julia_to_cpp(cobj),julia_to_cpp(frame),julia_to_cpp(mask)))
end
segment(cobj::dnn_SegmentationModel, frame::InputArray; mask::InputArray = (CxxMat())) = segment(cobj, frame, mask)

function dnn_SegmentationModel(model::String, config::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_SegmentationModel_cv_dnn_SegmentationModel_SegmentationModel(julia_to_cpp(model),julia_to_cpp(config)))
end
dnn_SegmentationModel(model::String; config::String = ("")) = dnn_SegmentationModel(model, config)

function dnn_SegmentationModel(network::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_SegmentationModel_cv_dnn_SegmentationModel_SegmentationModel(julia_to_cpp(network)))
end
function Base.getproperty(m::dnn_DetectionModel, s::Symbol)
    return Base.getfield(m, s)
end
function Base.setproperty!(m::dnn_DetectionModel, s::Symbol, v)
    return Base.setfield!(m, s, v)
end

function detect(cobj::dnn_DetectionModel, frame::InputArray, confThreshold::Float32, nmsThreshold::Float32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_DetectionModel_cv_dnn_DetectionModel_detect(julia_to_cpp(cobj),julia_to_cpp(frame),julia_to_cpp(confThreshold),julia_to_cpp(nmsThreshold)))
end
detect(cobj::dnn_DetectionModel, frame::InputArray; confThreshold::Float32 = Float32(0.5), nmsThreshold::Float32 = Float32(0.0)) = detect(cobj, frame, confThreshold, nmsThreshold)

function dnn_DetectionModel(model::String, config::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_DetectionModel_cv_dnn_DetectionModel_DetectionModel(julia_to_cpp(model),julia_to_cpp(config)))
end
dnn_DetectionModel(model::String; config::String = ("")) = dnn_DetectionModel(model, config)

function dnn_DetectionModel(network::dnn_Net)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_DetectionModel_cv_dnn_DetectionModel_DetectionModel(julia_to_cpp(network)))
end

# function getAvailableTargets(be::dnn_Backend)
# 	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_getAvailableTargets(julia_to_cpp(be)))
# end

function Net_readFromModelOptimizer(xml::String, bin::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_readFromModelOptimizer(julia_to_cpp(xml),julia_to_cpp(bin)))
end

function Net_readFromModelOptimizer(bufferModelConfig::Array{UInt8, 1}, bufferWeights::Array{UInt8, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_Net_readFromModelOptimizer(julia_to_cpp(bufferModelConfig),julia_to_cpp(bufferWeights)))
end

function readNetFromDarknet(cfgFile::String, darknetModel::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromDarknet(julia_to_cpp(cfgFile),julia_to_cpp(darknetModel)))
end
readNetFromDarknet(cfgFile::String; darknetModel::String = (String())) = readNetFromDarknet(cfgFile, darknetModel)

function readNetFromDarknet(bufferCfg::Array{UInt8, 1}, bufferModel::Array{UInt8, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromDarknet(julia_to_cpp(bufferCfg),julia_to_cpp(bufferModel)))
end
readNetFromDarknet(bufferCfg::Array{UInt8, 1}; bufferModel::Array{UInt8, 1} = (stdggvectoriUInt8kOP())) = readNetFromDarknet(bufferCfg, bufferModel)

function readNetFromCaffe(prototxt::String, caffeModel::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromCaffe(julia_to_cpp(prototxt),julia_to_cpp(caffeModel)))
end
readNetFromCaffe(prototxt::String; caffeModel::String = (String())) = readNetFromCaffe(prototxt, caffeModel)

function readNetFromCaffe(bufferProto::Array{UInt8, 1}, bufferModel::Array{UInt8, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromCaffe(julia_to_cpp(bufferProto),julia_to_cpp(bufferModel)))
end
readNetFromCaffe(bufferProto::Array{UInt8, 1}; bufferModel::Array{UInt8, 1} = (stdggvectoriUInt8kOP())) = readNetFromCaffe(bufferProto, bufferModel)

function readNetFromTensorflow(model::String, config::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromTensorflow(julia_to_cpp(model),julia_to_cpp(config)))
end
readNetFromTensorflow(model::String; config::String = (String())) = readNetFromTensorflow(model, config)

function readNetFromTensorflow(bufferModel::Array{UInt8, 1}, bufferConfig::Array{UInt8, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromTensorflow(julia_to_cpp(bufferModel),julia_to_cpp(bufferConfig)))
end
readNetFromTensorflow(bufferModel::Array{UInt8, 1}; bufferConfig::Array{UInt8, 1} = (stdggvectoriUInt8kOP())) = readNetFromTensorflow(bufferModel, bufferConfig)

function readNetFromTorch(model::String, isBinary::Bool, evaluate::Bool)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromTorch(julia_to_cpp(model),julia_to_cpp(isBinary),julia_to_cpp(evaluate)))
end
readNetFromTorch(model::String; isBinary::Bool = (true), evaluate::Bool = (true)) = readNetFromTorch(model, isBinary, evaluate)

function readNet(model::String, config::String, framework::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNet(julia_to_cpp(model),julia_to_cpp(config),julia_to_cpp(framework)))
end
readNet(model::String; config::String = (""), framework::String = ("")) = readNet(model, config, framework)

function readNet(framework::String, bufferModel::Array{UInt8, 1}, bufferConfig::Array{UInt8, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNet(julia_to_cpp(framework),julia_to_cpp(bufferModel),julia_to_cpp(bufferConfig)))
end
readNet(framework::String, bufferModel::Array{UInt8, 1}; bufferConfig::Array{UInt8, 1} = (stdggvectoriUInt8kOP())) = readNet(framework, bufferModel, bufferConfig)

function readTorchBlob(filename::String, isBinary::Bool)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readTorchBlob(julia_to_cpp(filename),julia_to_cpp(isBinary)))
end
readTorchBlob(filename::String; isBinary::Bool = (true)) = readTorchBlob(filename, isBinary)

function readNetFromModelOptimizer(xml::String, bin::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromModelOptimizer(julia_to_cpp(xml),julia_to_cpp(bin)))
end

function readNetFromModelOptimizer(bufferModelConfig::Array{UInt8, 1}, bufferWeights::Array{UInt8, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromModelOptimizer(julia_to_cpp(bufferModelConfig),julia_to_cpp(bufferWeights)))
end

function readNetFromONNX(onnxFile::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromONNX(julia_to_cpp(onnxFile)))
end

function readNetFromONNX(buffer::Array{UInt8, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readNetFromONNX(julia_to_cpp(buffer)))
end

function readTensorFromONNX(path::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_readTensorFromONNX(julia_to_cpp(path)))
end

function blobFromImage(image::InputArray, scalefactor::Float64, size::Size, mean::Scalar, swapRB::Bool, crop::Bool, ddepth::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_blobFromImage(julia_to_cpp(image),julia_to_cpp(scalefactor),julia_to_cpp(size),julia_to_cpp(mean),julia_to_cpp(swapRB),julia_to_cpp(crop),julia_to_cpp(ddepth)))
end
blobFromImage(image::InputArray; scalefactor::Float64 = Float64(1.0), size::Size = (SizeOP()), mean::Scalar = (), swapRB::Bool = (false), crop::Bool = (false), ddepth::Int32 = Int32(CV_32F)) = blobFromImage(image, scalefactor, size, mean, swapRB, crop, ddepth)

function blobFromImages(images::Array{InputArray, 1}, scalefactor::Float64, size::Size, mean::Scalar, swapRB::Bool, crop::Bool, ddepth::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_blobFromImages(julia_to_cpp(images),julia_to_cpp(scalefactor),julia_to_cpp(size),julia_to_cpp(mean),julia_to_cpp(swapRB),julia_to_cpp(crop),julia_to_cpp(ddepth)))
end
blobFromImages(images::Array{InputArray, 1}; scalefactor::Float64 = Float64(1.0), size::Size = (SizeOP()), mean::Scalar = (), swapRB::Bool = (false), crop::Bool = (false), ddepth::Int32 = Int32(CV_32F)) = blobFromImages(images, scalefactor, size, mean, swapRB, crop, ddepth)

function imagesFromBlob(blob_::InputArray, images_::Array{InputArray, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_imagesFromBlob(julia_to_cpp(blob_),julia_to_cpp(images_)))
end
imagesFromBlob(blob_::InputArray; images_::Array{InputArray, 1} = (Array{InputArray, 1}())) = imagesFromBlob(blob_, images_)

function shrinkCaffeModel(src::String, dst::String, layersTypes::Array{String, 1})
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_shrinkCaffeModel(julia_to_cpp(src),julia_to_cpp(dst),julia_to_cpp(layersTypes)))
end
shrinkCaffeModel(src::String, dst::String; layersTypes::Array{String, 1} = (stdggvectoriStringkOP())) = shrinkCaffeModel(src, dst, layersTypes)

function writeTextGraph(model::String, output::String)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_writeTextGraph(julia_to_cpp(model),julia_to_cpp(output)))
end

function NMSBoxes(bboxes::Array{Rect{Float64}, 1}, scores::Array{Float32, 1}, score_threshold::Float32, nms_threshold::Float32, eta::Float32, top_k::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_NMSBoxes(julia_to_cpp(bboxes),julia_to_cpp(scores),julia_to_cpp(score_threshold),julia_to_cpp(nms_threshold),julia_to_cpp(eta),julia_to_cpp(top_k)))
end
NMSBoxes(bboxes::Array{Rect{Float64}, 1}, scores::Array{Float32, 1}, score_threshold::Float32, nms_threshold::Float32; eta::Float32 = Float32(1.0), top_k::Int32 = Int32(0)) = NMSBoxes(bboxes, scores, score_threshold, nms_threshold, eta, top_k)

function NMSBoxesRotated(bboxes::Array{RotatedRect, 1}, scores::Array{Float32, 1}, score_threshold::Float32, nms_threshold::Float32, eta::Float32, top_k::Int32)
	return cpp_to_julia(jlopencv_cv_dnn_cv_dnn_NMSBoxes(julia_to_cpp(bboxes),julia_to_cpp(scores),julia_to_cpp(score_threshold),julia_to_cpp(nms_threshold),julia_to_cpp(eta),julia_to_cpp(top_k)))
end
NMSBoxesRotated(bboxes::Array{RotatedRect, 1}, scores::Array{Float32, 1}, score_threshold::Float32, nms_threshold::Float32; eta::Float32 = Float32(1.0), top_k::Int32 = Int32(0)) = NMSBoxesRotated(bboxes, scores, score_threshold, nms_threshold, eta, top_k)


text:
"""
    gammaHSIC(X::Array{T, 2}, Y::Array{T, 2})

"""
function gammaHSIC(X::Array{T, 2}, Y::Array{T, 2}; α = 0.1, randomSubSet = 100, kernelSize = -1) where T <: Real

	M = size(X)[1]
    #@assert M > 1 "This implementation does not support univariate cases."

	# get kernel sizes
	if kernelSize == -1
		sigx = estimateKernelSize(X, sampleSize = randomSubSet)
		sigy = estimateKernelSize(Y, sampleSize = randomSubSet)
	else
		sigx = kernelSize
		sigy = kernelSize
	end

	bone = ones(T, M, 1)
    H = eye(M) - (1/M)*ones(T, M,M)

	K = rbfDotProduct(X, X, sigx)
	L = rbfDotProduct(Y, Y, sigy)

	# NOTE: these are slightly biased estimates of centred Gram matrices
	Kc = H*K*H
	Lc = H*L*H

	# NOTE: we fit Gamma to testStat*M
	testStat = 1/M * sum(Kc'.*Lc)

	varHSIC = (1/6 * Kc.*Lc).^2
	varHSIC = 1/M/(M-1)* (  sum(varHSIC) - sum(diag(varHSIC))  ) # second subtracted term is bias correction
	varHSIC = 72*(M-4)*(M-5)/M/(M-1)/(M-2)/(M-3) * varHSIC

	K = K-diagm(0 => diag(K))
	L = L-diagm(0 => diag(L))

    muX = (1/M/(M-1)) *bone'*(K*bone)
    muY = (1/M/(M-1)) *bone'*(L*bone)

	mHSIC  = 1/M * ( 1 .+ muX*muY - muX - muY )

	al = mHSIC^2 / varHSIC
	bet = varHSIC*M ./ mHSIC

	return (testStat, gammainvcdf(al[1], bet[1], 1 - α))
end


text:
"""
	$(SIGNATURES)

Shocks destroy a fraction (college specific) of what has been learned this year.
"""
Base.@kwdef mutable struct HcShocksSwitchesLoseLearn <: HcShockSwitches
    nColleges :: CollInt
    fracLostV :: Vector{Double}
end

mutable struct HcShockSetLoseLearn <: HcShockSet
    objId :: ObjectId
    switches :: HcShocksSwitchesLoseLearn
    shockProbs :: BoundedVector{Double}
end

struct HcShockLoseLearn <: HcShock
    fracLost :: Double
    prob :: Double
end

function Base.show(io :: IO, dh :: HcShocksSwitchesLoseLearn)
    print(io, typeof(dh));
end

# Set
function Base.show(io :: IO, dh :: HcShockSetLoseLearn)
    print(io, typeof(dh));
end

Base.show(io :: IO, dh :: HcShockLoseLearn) = print(io, typeof(dh));

StructLH.describe(dh :: HcShockSetLoseLearn) = StructLH.describe(dh.switches);
function StructLH.describe(dh :: HcShocksSwitchesLoseLearn)
    m = [
        "H shocks"  "Lose a fraction of what has been learned";
    ];
    return m
end

function settings_table(dh :: HcShocksSwitchesLoseLearn)
    return StructLH.describe(dh)
end

## -------------  Access

has_pvector(dh :: HcShockSetLoseLearn) = false;
# get_pvector(dh :: HcShockSetLoseLearn) = dh.switches.pvec;

# Shock prob by college
shock_probs(dh :: HcShockSetLoseLearn) = ModelParams.values(dh.shockProbs);
shock_probs(dh :: HcShockSetLoseLearn, idx) = shock_probs(dh)[idx];

shock_prob(dh :: HcShockLoseLearn) = dh.prob;

## -------------  Construction

function make_hshock_switches_lose_learn(nColleges, fracLostV)
    switches = HcShocksSwitchesLoseLearn(nColleges, fracLostV);
    @assert validate_hshocks(switches)
    return switches
end

function validate_hshocks(switches :: HcShocksSwitchesLoseLearn)
    isValid = true;
    isValid = isValid  &&  all(0.0 .<= switches.fracLostV .< 1.0);
    isValid = isValid  &&  (size(switches.fracLostV) == (switches.nColleges, ));
    return isValid
end

function make_h_shock_set(objId, switches :: HcShocksSwitchesLoseLearn)
    pProbId = make_child_id(objId, :shockProbs, "H shock probabilities");
    pProb = init_shock_probs(pProbId, switches);
    return HcShockSetLoseLearn(objId, switches, pProb)
end

# update description ++++++
function init_shock_probs(objId :: ObjectId, switches :: HcShocksSwitchesLoseLearn)
    nc = n_colleges(switches);
    dxV = fill(0.2, nc);
    pProb = BoundedVector{Double}(objId, ParamVector(objId),
        :decreasing, 0.01, 0.5, dxV);
    set_pvector!(pProb; description = ldescription(:hShockGrad),
        symbol = lsymbol(:hShockGrad), isCalibrated = true);
    return pProb
end

function validate_hshocks(dh :: HcShockSetLoseLearn)
    isValid = all_at_least(shock_probs(dh), 0.0)  .&
        all_at_most(shock_probs(dh), 1.0)
    return isValid
end

function make_h_shock(dh :: HcShockSetLoseLearn, iCollege :: Integer)
    return HcShockLoseLearn(dh.switches.fracLostV[iCollege],
        shock_probs(dh, iCollege));
end

function validate_hshocks(dh :: HcShockLoseLearn)
    return (shock_prob(dh) >= 0.0)  &&  (shock_prob(dh) <= 1.0)
end

# ----------------  Individual shock object for one college

"""
	$(SIGNATURES)

Simulate h shocks. Returns `h` next period, given `h` at the end of this period (before shocks are realized).

Important to draw random variables in the same order each time.

# Arguments
- `dh`: shock object
- `hEndIdxM`: h indices at end of period
"""
function sim_h_shocks(dh :: HcShockLoseLearn, hStartM, hEndIdxM :: Array{T1},
    hGridNextV, rng :: AbstractRNG) where T1 <: Integer

    @assert size(hStartM) == size(hEndIdxM);
    @assert all_at_least(hEndIdxM, 1);
    @assert all_at_most(hEndIdxM, length(hGridNextV));

    hOutM = similar(hEndIdxM);
    for (j, hEndIdx) in enumerate(hEndIdxM)
        hEnd = hGridNextV[hEndIdx];
        withShock = (rand(rng) < shock_prob(dh));
        hOutM[j] = hprime_idx(dh, hStartM[j], hEnd, hGridNextV; withShock);
    end
    return hOutM
end

function hprime_idx(dh :: HcShockLoseLearn, hStart, hEnd, hGridNextV; withShock :: Bool)
    if withShock
        fracKeep = 1.0 - dh.fracLost;
    else
        fracKeep = 1.0;
    end
    learn = fracKeep .* (hEnd .- hStart);
    return round_to_grid(hStart .+ learn, hGridNextV)
end

# Copied from MatrixLH

"""
round_to_grid()

Round a matrix to the nearest points on a grid. Returns the index, not the value.
"""
function round_to_grid(xM :: Array{T}, gridV :: AbstractVector{T}) where T <: Real
    idxM = similar(xM, Int);
    for (i1, x) in enumerate(xM)
        idxM[i1] = round_to_grid(x, gridV);
    end
    return idxM
end

function round_to_grid(x :: F1, gridV :: AbstractVector{F2}) where {F1 <: Real, F2 <: Real}
    # _, idx = findmin(abs.(gridV .- x));
    _, idx = fndmin(z -> abs(z - x),  gridV)
    return idx
end

# findmin for functions.
function fndmin(f, xV)
    fMin = f(first(xV));
    jMin = 1;
    for (j, x) in enumerate(xV)
        fVal = f(x);
        if fVal < fMin
            fMin = fVal;
            jMin = j;
        end
    end
    return fMin, jMin
end

"""
	$(SIGNATURES)

Prob(h on tomorrow's grid | h grid point today). Mainly for testing.
"""
function pr_hprime(dh :: HcShockLoseLearn, hEndIdx :: I1,
    hStart :: Real, hEnd :: Real, hGridNextV) where I1 <: Integer

    prV = (shock_prob(dh), 1.0 - shock_prob(dh));
    idxV = (hprime_idx(dh, hStart, hEnd, hGridNextV; withShock = true),
        hprime_idx(dh, hStart, hEnd, hGridNextV; withShock = false));
        # for hasShock in (true, false)];
    return idxV, prV
end

# function time_pr_hprime

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


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9]);F = [5162*x1^2*x3^2+638*x1^2*x3*x6+6241*x1^2*x4^2-1580*x1^2*x4*x5+100*x1^2*x5^2+4589*x1^2*x6^2+6724*x1^2*x7^2-11644*x1^2*x7*x8+5041*x1^2*x8^2-10676*x1*x2*x3*x6+8532*x1*x2*x4^2-1080*x1*x2*x4*x5+4104*x1*x2*x5*x6-2128*x1*x2*x6*x9+2272*x1*x3^2*x8-1034*x1*x3^2*x9+3496*x1*x3*x4*x6-1716*x1*x3*x4*x7-13206*x1*x3*x4*x9+1136*x1*x3*x5^2-4416*x1*x3*x5*x6-7776*x1*x3*x6^2-2726*x1*x3*x6*x9-12070*x1*x3*x7*x8-1782*x1*x3*x9^2+5904*x1*x4^2*x7-5112*x1*x4^2*x8-5672*x1*x4*x6*x7+994*x1*x4*x6*x8+13386*x1*x4*x6*x9-316*x1*x4*x7*x8-5214*x1*x4*x8^2+36*x1*x5*x6*x9+40*x1*x5*x7*x8+492*x1*x5*x7*x9+660*x1*x5*x8^2-426*x1*x5*x8*x9-6714*x1*x6*x9^2+3844*x2^2*x3^2-6696*x2^2*x3*x5+3472*x2^2*x3*x9+2916*x2^2*x4^2+2916*x2^2*x5^2-3024*x2^2*x5*x9+3700*x2^2*x6^2+5632*x2^2*x6*x9+7089*x2^2*x9^2-5704*x2*x3^2*x4+4968*x2*x3*x4*x5-2576*x2*x3*x4*x9-1344*x2*x3*x6*x8-8084*x2*x4^2*x9+15912*x2*x4*x6*x9-216*x2*x4*x7*x8-3564*x2*x4*x8^2+7332*x2*x4*x9^2-672*x2*x5^2*x6+4042*x2*x5*x9^2+616*x2*x6^2*x7+1408*x2*x6^2*x9+7140*x2*x6*x7*x8+5120*x2*x6*x7*x9+2048*x2*x6*x9^2+6144*x2*x7*x9^2+5828*x2*x9^3+2116*x3^2*x4^2+2116*x3^2*x5^2+7452*x3^2*x5*x6+6561*x3^2*x6^2+256*x3^2*x8^2+2209*x3^2*x9^2-1472*x3*x4*x5*x9-2592*x3*x4*x6*x9+7332*x3*x4*x7*x9-2976*x3*x4*x8*x9+256*x3*x5^2*x8-1472*x3*x5^2*x9-2592*x3*x5*x6*x9+1932*x3*x5*x9^2+3402*x3*x6*x9^2-2720*x3*x7*x8^2+7614*x3*x9^3+8692*x4^4-504*x4^3*x6-13416*x4^3*x9-7180*x4^2*x5*x9+49*x4^2*x6^2+6084*x4^2*x7^2+4325*x4^2*x9^2-1488*x4*x5^2*x9-42*x4*x5*x6*x9+7220*x4*x5*x9^2+15810*x4*x7*x8*x9+12636*x4*x7*x9^2+9000*x4*x9^3+64*x5^4-1360*x5^2*x7*x8+2114*x5^2*x9^2+4660*x5*x9^3+49*x6^2*x7^2+224*x6^2*x7*x9+5881*x6^2*x9^2+672*x6*x7^2*x9-300*x6*x7*x8*x9+1536*x6*x7*x9^2-4950*x6*x8^2*x9+7229*x7^2*x8^2+2304*x7^2*x9^2+132*x7*x8^3+1089*x8^4+10846*x9^4+2068*x1*x2*x9-14432*x1*x3*x7+12496*x1*x3*x8+792*x1*x3*x9-3784*x1*x4^2+3432*x1*x4*x9+1892*x1*x5*x9-6840*x1*x6*x7+2088*x1*x6*x9+2728*x1*x9^2+11160*x2*x3*x7-3520*x2*x5*x6-9720*x2*x5*x7-5120*x2*x5*x9+5040*x2*x7*x9-6336*x3*x4^2+1232*x3*x4*x6-8280*x3*x4*x7-528*x3*x5*x9-3384*x3*x9^2-5616*x4*x7*x9-560*x5*x6*x7-1280*x5*x6*x9-3840*x5*x7*x9-5832*x9^3+484*x1^2+7744*x3^2+1600*x5^2+8100*x7^2+1296*x9^2, -10676*x1*x3*x6+8532*x1*x4^2-1080*x1*x4*x5+4104*x1*x5*x6-2128*x1*x6*x9+7688*x2*x3^2-13392*x2*x3*x5+6944*x2*x3*x9+5832*x2*x4^2+5832*x2*x5^2-6048*x2*x5*x9+7400*x2*x6^2+11264*x2*x6*x9+14178*x2*x9^2-5704*x3^2*x4+4968*x3*x4*x5-2576*x3*x4*x9-1344*x3*x6*x8-8084*x4^2*x9+15912*x4*x6*x9-216*x4*x7*x8-3564*x4*x8^2+7332*x4*x9^2-672*x5^2*x6+4042*x5*x9^2+616*x6^2*x7+1408*x6^2*x9+7140*x6*x7*x8+5120*x6*x7*x9+2048*x6*x9^2+6144*x7*x9^2+5828*x9^3+2068*x1*x9+11160*x3*x7-3520*x5*x6-9720*x5*x7-5120*x5*x9+5040*x7*x9, 10324*x1^2*x3+638*x1^2*x6-10676*x1*x2*x6+4544*x1*x3*x8-2068*x1*x3*x9+3496*x1*x4*x6-1716*x1*x4*x7-13206*x1*x4*x9+1136*x1*x5^2-4416*x1*x5*x6-7776*x1*x6^2-2726*x1*x6*x9-12070*x1*x7*x8-1782*x1*x9^2+7688*x2^2*x3-6696*x2^2*x5+3472*x2^2*x9-11408*x2*x3*x4+4968*x2*x4*x5-2576*x2*x4*x9-1344*x2*x6*x8+4232*x3*x4^2+4232*x3*x5^2+14904*x3*x5*x6+13122*x3*x6^2+512*x3*x8^2+4418*x3*x9^2-1472*x4*x5*x9-2592*x4*x6*x9+7332*x4*x7*x9-2976*x4*x8*x9+256*x5^2*x8-1472*x5^2*x9-2592*x5*x6*x9+1932*x5*x9^2+3402*x6*x9^2-2720*x7*x8^2+7614*x9^3-14432*x1*x7+12496*x1*x8+792*x1*x9+11160*x2*x7-6336*x4^2+1232*x4*x6-8280*x4*x7-528*x5*x9-3384*x9^2+15488*x3, 12482*x1^2*x4-1580*x1^2*x5+17064*x1*x2*x4-1080*x1*x2*x5+3496*x1*x3*x6-1716*x1*x3*x7-13206*x1*x3*x9+11808*x1*x4*x7-10224*x1*x4*x8-5672*x1*x6*x7+994*x1*x6*x8+13386*x1*x6*x9-316*x1*x7*x8-5214*x1*x8^2+5832*x2^2*x4-5704*x2*x3^2+4968*x2*x3*x5-2576*x2*x3*x9-16168*x2*x4*x9+15912*x2*x6*x9-216*x2*x7*x8-3564*x2*x8^2+7332*x2*x9^2+4232*x3^2*x4-1472*x3*x5*x9-2592*x3*x6*x9+7332*x3*x7*x9-2976*x3*x8*x9+34768*x4^3-1512*x4^2*x6-40248*x4^2*x9-14360*x4*x5*x9+98*x4*x6^2+12168*x4*x7^2+8650*x4*x9^2-1488*x5^2*x9-42*x5*x6*x9+7220*x5*x9^2+15810*x7*x8*x9+12636*x7*x9^2+9000*x9^3-7568*x1*x4+3432*x1*x9-12672*x3*x4+1232*x3*x6-8280*x3*x7-5616*x7*x9, -1580*x1^2*x4+200*x1^2*x5-1080*x1*x2*x4+4104*x1*x2*x6+2272*x1*x3*x5-4416*x1*x3*x6+36*x1*x6*x9+40*x1*x7*x8+492*x1*x7*x9+660*x1*x8^2-426*x1*x8*x9-6696*x2^2*x3+5832*x2^2*x5-3024*x2^2*x9+4968*x2*x3*x4-1344*x2*x5*x6+4042*x2*x9^2+4232*x3^2*x5+7452*x3^2*x6-1472*x3*x4*x9+512*x3*x5*x8-2944*x3*x5*x9-2592*x3*x6*x9+1932*x3*x9^2-7180*x4^2*x9-2976*x4*x5*x9-42*x4*x6*x9+7220*x4*x9^2+256*x5^3-2720*x5*x7*x8+4228*x5*x9^2+4660*x9^3+1892*x1*x9-3520*x2*x6-9720*x2*x7-5120*x2*x9-528*x3*x9-560*x6*x7-1280*x6*x9-3840*x7*x9+3200*x5, 638*x1^2*x3+9178*x1^2*x6-10676*x1*x2*x3+4104*x1*x2*x5-2128*x1*x2*x9+3496*x1*x3*x4-4416*x1*x3*x5-15552*x1*x3*x6-2726*x1*x3*x9-5672*x1*x4*x7+994*x1*x4*x8+13386*x1*x4*x9+36*x1*x5*x9-6714*x1*x9^2+7400*x2^2*x6+5632*x2^2*x9-1344*x2*x3*x8+15912*x2*x4*x9-672*x2*x5^2+1232*x2*x6*x7+2816*x2*x6*x9+7140*x2*x7*x8+5120*x2*x7*x9+2048*x2*x9^2+7452*x3^2*x5+13122*x3^2*x6-2592*x3*x4*x9-2592*x3*x5*x9+3402*x3*x9^2-504*x4^3+98*x4^2*x6-42*x4*x5*x9+98*x6*x7^2+448*x6*x7*x9+11762*x6*x9^2+672*x7^2*x9-300*x7*x8*x9+1536*x7*x9^2-4950*x8^2*x9-6840*x1*x7+2088*x1*x9-3520*x2*x5+1232*x3*x4-560*x5*x7-1280*x5*x9, 13448*x1^2*x7-11644*x1^2*x8-1716*x1*x3*x4-12070*x1*x3*x8+5904*x1*x4^2-5672*x1*x4*x6-316*x1*x4*x8+40*x1*x5*x8+492*x1*x5*x9-216*x2*x4*x8+616*x2*x6^2+7140*x2*x6*x8+5120*x2*x6*x9+6144*x2*x9^2+7332*x3*x4*x9-2720*x3*x8^2+12168*x4^2*x7+15810*x4*x8*x9+12636*x4*x9^2-1360*x5^2*x8+98*x6^2*x7+224*x6^2*x9+1344*x6*x7*x9-300*x6*x8*x9+1536*x6*x9^2+14458*x7*x8^2+4608*x7*x9^2+132*x8^3-14432*x1*x3-6840*x1*x6+11160*x2*x3-9720*x2*x5+5040*x2*x9-8280*x3*x4-5616*x4*x9-560*x5*x6-3840*x5*x9+16200*x7, -11644*x1^2*x7+10082*x1^2*x8+2272*x1*x3^2-12070*x1*x3*x7-5112*x1*x4^2+994*x1*x4*x6-316*x1*x4*x7-10428*x1*x4*x8+40*x1*x5*x7+1320*x1*x5*x8-426*x1*x5*x9-1344*x2*x3*x6-216*x2*x4*x7-7128*x2*x4*x8+7140*x2*x6*x7+512*x3^2*x8-2976*x3*x4*x9+256*x3*x5^2-5440*x3*x7*x8+15810*x4*x7*x9-1360*x5^2*x7-300*x6*x7*x9-9900*x6*x8*x9+14458*x7^2*x8+396*x7*x8^2+4356*x8^3+12496*x1*x3, -2128*x1*x2*x6-1034*x1*x3^2-13206*x1*x3*x4-2726*x1*x3*x6-3564*x1*x3*x9+13386*x1*x4*x6+36*x1*x5*x6+492*x1*x5*x7-426*x1*x5*x8-13428*x1*x6*x9+3472*x2^2*x3-3024*x2^2*x5+5632*x2^2*x6+14178*x2^2*x9-2576*x2*x3*x4-8084*x2*x4^2+15912*x2*x4*x6+14664*x2*x4*x9+8084*x2*x5*x9+1408*x2*x6^2+5120*x2*x6*x7+4096*x2*x6*x9+12288*x2*x7*x9+17484*x2*x9^2+4418*x3^2*x9-1472*x3*x4*x5-2592*x3*x4*x6+7332*x3*x4*x7-2976*x3*x4*x8-1472*x3*x5^2-2592*x3*x5*x6+3864*x3*x5*x9+6804*x3*x6*x9+22842*x3*x9^2-13416*x4^3-7180*x4^2*x5+8650*x4^2*x9-1488*x4*x5^2-42*x4*x5*x6+14440*x4*x5*x9+15810*x4*x7*x8+25272*x4*x7*x9+27000*x4*x9^2+4228*x5^2*x9+13980*x5*x9^2+224*x6^2*x7+11762*x6^2*x9+672*x6*x7^2-300*x6*x7*x8+3072*x6*x7*x9-4950*x6*x8^2+4608*x7^2*x9+43384*x9^3+2068*x1*x2+792*x1*x3+3432*x1*x4+1892*x1*x5+2088*x1*x6+5456*x1*x9-5120*x2*x5+5040*x2*x7-528*x3*x5-6768*x3*x9-5616*x4*x7-1280*x5*x6-3840*x5*x7-17496*x9^2+2592*x9];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
"""
Generates various Genie files.
"""
module Generator

using Revise, Genie, Genie.Loggers, Genie.FileTemplates, Genie.Inflector, Genie.Configuration #, Genie.App

"""
    new_controller(cmd_args::Dict{String,Any}) :: Nothing

Generates a new Genie controller file and persists it to the resources folder.
"""
function new_controller(cmd_args::Dict{String,Any}) :: Nothing
  resource_name = cmd_args["controller:new"]
  Genie.Inflector.is_singular(resource_name) && (resource_name = Inflector.to_plural(resource_name) |> Base.get)
  resource_name = uppercasefirst(resource_name)

  resource_path = setup_resource_path(resource_name)
  cfn = controller_file_name(resource_name)
  write_resource_file(resource_path, cfn, resource_name, :controller) &&
    log("New controller created at $(joinpath(resource_path, cfn))")

  nothing
end

"""
    new_channel(cmd_args::Dict{String,Any}) :: Nothing

Generates a new Genie channel file and persists it to the resources folder.
"""
function new_channel(cmd_args::Dict{String,Any}) :: Nothing
  resource_name = uppercasefirst(cmd_args["channel:new"])
  if Genie.Inflector.is_singular(resource_name)
    resource_name = Genie.Inflector.to_plural(resource_name) |> Base.get
  end

  resource_path = setup_resource_path(resource_name)
  cfn = channel_file_name(resource_name)
  write_resource_file(resource_path, cfn, resource_name, :channel) &&
    log("New channel created at $(joinpath(resource_path, cfn))")

  nothing
end

"""
    new_resource(cmd_args::Dict{String,Any}, config::Settings) :: Nothing
    new_resource(resource_name::Union{String,Symbol}) :: Nothing

Generates all the files associated with a new resource and persists them to the resources folder.
"""
function new_resource(cmd_args::Dict{String,Any}) :: Nothing
  resource_name = uppercasefirst(cmd_args["resource:new"])

  try
    SearchLight.Generator.new_resource(resource_name)
  catch ex
    log("Skipping SearchLight", :warn)
  end

  if Genie.Inflector.is_singular(resource_name)
    resource_name = Genie.Inflector.to_plural(resource_name) |> Base.get
  end

  resource_path = setup_resource_path(resource_name)
  for (resource_file, resource_type) in [(controller_file_name(resource_name), :controller), (channel_file_name(resource_name), :channel), (Genie.GENIE_AUTHORIZATOR_FILE_NAME, :authorizer)]
    write_resource_file(resource_path, resource_file, resource_name, resource_type) &&
      log("New $resource_file created at $(joinpath(resource_path, resource_file))")
  end

  views_path = joinpath(resource_path, "views")
  ! isdir(views_path) && mkpath(views_path)

  ! isdir(Genie.TEST_PATH_UNIT) && mkpath(Genie.TEST_PATH_UNIT)
  test_file = resource_name * Genie.TEST_FILE_IDENTIFIER |> lowercase
  write_resource_file(Genie.TEST_PATH_UNIT, test_file, resource_name, :test) &&
    log("New $test_file created at $(joinpath(Genie.TEST_PATH_UNIT, test_file))")

  nothing
end
function new_resource(resource_name::Union{String,Symbol}) :: Nothing
  new_resource(Dict{String,Any}("resource:new" => string(resource_name)))
end

"""
    setup_resource_path(resource_name::String) :: String

Computes and creates the directories structure needed to persist a new resource.
"""
function setup_resource_path(resource_name::String) :: String
  resource_path = joinpath(Genie.RESOURCES_PATH, lowercase(resource_name))

  if ! isdir(resource_path)
    mkpath(resource_path)
    push!(LOAD_PATH, resource_path)
  end

  resource_path
end

"""
    write_resource_file(resource_path::String, file_name::String, resource_name::String) :: Bool

Generates all resouce files and persists them to disk.
"""
function write_resource_file(resource_path::String, file_name::String, resource_name::String, resource_type::Symbol) :: Bool
  resource_name = Base.get(Inflector.to_plural(resource_name)) |> Inflector.from_underscores

  try
    if resource_type == :controller
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, Genie.FileTemplates.new_controller(resource_name))
      end

    elseif resource_type == :authorizer
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, Genie.FileTemplates.new_authorizer())
      end

    elseif resource_type == :test
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, Genie.FileTemplates.new_test(resource_name, Base.get(Inflector.to_singular(resource_name)) ))
      end

    elseif resource_type == :channel
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, Genie.FileTemplates.new_channel(resource_name))
      end

    else
      error("Not supported, $file_name")
    end
  catch ex
    log(ex, :warn)
  end

  try
    Main.load_resources()
  catch ex
    log("Not in app, skipping autoload", :warn)
  end

  true
end

"""
"""
function resource_does_not_exist(resource_path::String, file_name::String) :: Bool
  if isfile(joinpath(resource_path, file_name))
    log("File already exists, $(joinpath(resource_path, file_name)) - skipping", :warn)
    return false
  end

  true
end

"""
"""
function controller_file_name(resource_name::Union{String,Symbol})
  string(resource_name) * Genie.GENIE_CONTROLLER_FILE_POSTFIX
end

"""
"""
function channel_file_name(resource_name::Union{String,Symbol})
  string(resource_name) * Genie.GENIE_CHANNEL_FILE_POSTFIX
end

end


text:
"""
Assumes x,y points go around the polygon in one direction.
Credits: https://rosettacode.org/wiki/Shoelace_formula_for_polygonal_area#Julia
"""
shoelacearea(x, y) =
    abs(sum(i * j for (i, j) in zip(x, append!(y[2:end], y[1]))) -
        sum(i * j for (i, j) in zip(append!(x[2:end], x[1]), y))) / 2


text:
using Query
using DataFrames
using Test

@testset "Standalone Syntax" begin

@testset "@take operator" begin
    df = DataFrame(a=[1,2,3], b=[3.,2.,1.], c=["a", "b", "c"])

    df2 = df |> @take(2) |> DataFrame

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [1,2]
    @test df2[:b] == [3.,2.]
    @test df2[:c] == ["a", "b"]

    df2 = DataFrame(@take(df, 2))

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [1,2]
    @test df2[:b] == [3.,2.]
    @test df2[:c] == ["a", "b"]
end

@testset "@drop operator" begin
    df = DataFrame(a=[1,2,3], b=[3.,2.,1.], c=["a", "b", "c"])

    df2 = df |> @drop(1) |> DataFrame

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [2,3]
    @test df2[:b] == [2.,1.]
    @test df2[:c] == ["b","c"]

    df2 = DataFrame(@drop(df, 1))

    @test df2 isa DataFrame
    @test size(df2) == (2,3)
    @test df2[:a] == [2,3]
    @test df2[:b] == [2.,1.]
    @test df2[:c] == ["b","c"]
end

@testset "@unique operator" begin
    df = DataFrame(a=[1,2,1], b=[3.,3.,3.])

    @test df |> @unique() |> collect == [(a=1,b=3.), (a=2,b=3.)]
    @test df |> @unique(_.b) |> collect == [(a=1,b=3.)]
end

end


text:
#####
##### This is the header for a region of memory in the heap.
#####

@enum Pool::UInt Local=0 Remote=1
struct PoolType{T} end

# It is imlemented as just a pointer, which points to the beginning of a region of memory.
# However, we reserve the first 64 bytes (cache line) for use by the header.
#
# These fields are accessible using `getproperty` and `setproperty` overloading.
#
# The layout looks like this
#
# bytes 7..0
#    The top 58 bits describe the size of the memory for this block
#    The lower 6 bits are reserved for misc metadata.
#       We can save the lower 6 bits because we assume sizes are always multiples of 64 bytes.
#       TODO: We can actually do better since sizes will be multiples of MIN_ALLOCATION.
#
#       bit 0: 1 if the block is free, 0 if it is not free.
#       bit 1-2: Pool ID: 00 if Local, 01 if Remote. These bits correspond to the `Pool` enum.
#       bit 3: dirty bit: 0 if clean, 1 if dirty
#       bit 4: evicting bit: 0 if normal free, 1 if being evicted
#       bit 5: queued bit: 1 if queued for freeing, 0 otherwise
#       bit 6: orphaned bit: 1 is freed while evicted.
#
#
#    Size can be obtained with `block.size`.
#    The free state is obtained with `block.free`.
#    All of the lower 6 bits are obtained with `block.bitmasks`.
#
# bytes 15..8
#    if the block is free, this is a pointer to the next block in the freelist.
#    Access with `block.next`.
#
#    If the block is NOT free, this is the ID, of the owner of the block.
#    Access this with `block.id`
#
# bytes 23..16
#    if the block is free, this is a pointer to the previous block in the freelist.
#    access with `block.previous`.
#
#    if block is not free, this is a pointer to the sibling array.
#    access with 'block.sibling'.
#    If 'isnull(block.sibling)', then the sibling does not exist.
#    A sibling is a copy of the array that lives in another memory pool.
#
# bytes 31..24
#    the size of the block before this block.
#    used to walk backward through the heap/

struct Block
    ptr::Ptr{Nothing}
    # Inner constructor for ambiguity resolution
    Block(ptr::Ptr{Nothing}) = new(ptr)
end

Block() = Block(Ptr{Nothing}(0))
Block(ptr::Ptr) = Block(convert(Ptr{Nothing}, ptr))
Block(address::UInt) = Block(Ptr{Nothing}(address))

Base.length(block::Block) = block.size - headersize()
Base.sizeof(block::Block) = block.size

Base.pointer(x::Block) = getfield(x, :ptr)
datapointer(x::Block) = pointer(x) + headersize()
isnull(x::Block) = isnull(x.ptr)
Base.isless(a::Block, b::Block) = a.ptr < b.ptr
address(x::Block) = convert(UInt, pointer(x))

Base.hash(x::Block, u::UInt = UInt(0x178e9aaf)) = hash(pointer(x), u)
Base.:(==)(a::Block, b::Block) = (pointer(a) == pointer(b))

unsafe_block(ptr::Ptr) = Block(convert(Ptr{Nothing}, ptr) - headersize())

struct BlockMeta end
metadata(x, ::BlockMeta) = unsafe_block(unsafe_pointer(x))

# Reserved room in for each allocation (bytes).
headersize() = 64

function zerometa!(block::Block)
    ptr = Ptr{UInt8}(pointer(block))
    for i in Base.OneTo(headersize())
        @_ivdep_meta
        unsafe_store!(ptr, zero(UInt8), i)
    end
    return nothing
end

mask(x) = one(x) << x
mask(x, y...) = mask(x) | mask(y...)

macro getbits(typ, ptr, x::Integer...)
    # Precompute the mask and minimum value
    msk = mask(x...)
    min = minimum(x)
    return :((unsafe_load(convert(Ptr{$typ}, $(esc(ptr)))) & $msk) >> $min)
end

macro setbits!(typ, ptr, v, x::Integer...)
    load_clear_mask = ~mask(x...)
    min = minimum(x)
    set_clear_mask = mask((x .- min)...)
    return quote
        cptr = convert(Ptr{$typ}, $(esc(ptr)))
        y = unsafe_load(cptr) & $load_clear_mask
        y |= ($(esc(v)) & $set_clear_mask) << $min
        unsafe_store!(cptr, y)
    end
end

const LOWER_BIT_MASK = UInt(0x7F)

# Accessing "x.size" automatically clears the lower bits used for masking metadata.
# Here, we just store the size directly back to clear the lower bits.
clearbits!(x::Block) = unsafe_store!(convert(Ptr{UInt64}, pointer(x)), x.size)

function Base.getproperty(x::Block, name::Symbol)
    if name == :ptr
        return pointer(x)
    elseif name == :size
        # Load the first 8 bytes which encodes the size as a UInt64
        # Mask out the lower 6 bits since those are reserved for other metadata.
        return unsafe_load(convert(Ptr{UInt64}, pointer(x))) & ~LOWER_BIT_MASK

    # Bitmask metadata
    elseif name == :bitmasks
        return unsafe_load(convert(Ptr{UInt64}, pointer(x))) & LOWER_BIT_MASK
    elseif name == :free
        return Bool(@getbits(UInt, pointer(x), 0))
    elseif name == :pool
        return Pool(@getbits(UInt, pointer(x), 1, 2))
    elseif name == :dirty
        return Bool(@getbits(UInt, pointer(x), 3))
    elseif name == :evicting
        return Bool(@getbits(UInt, pointer(x), 4))
    elseif name == :queued
        return Bool(@getbits(UInt, pointer(x), 5))
    elseif name == :orphaned
        return Bool(@getbits(UInt, pointer(x), 6))

    # Bytes 15..8
    elseif name == :next
        return unsafe_load(convert(Ptr{Block}, pointer(x) + 8))

    elseif name == :id
        return unsafe_load(convert(Ptr{UInt64}, pointer(x) + 8))

    # Bytes 23..16
    elseif name == :previous || name == :sibling
        return unsafe_load(convert(Ptr{Block}, pointer(x) + 16))

    # Bytes 31..24
    elseif name == :backsize
        return unsafe_load(convert(Ptr{UInt64}, pointer(x) + 24))

    else
        error()
    end
end

function Base.setproperty!(x::Block, name::Symbol, v)
    # Bytes 7..0
    if name == :size
        sz = convert(UInt64, v) & ~LOWER_BIT_MASK
        unsafe_store!(convert(Ptr{UInt64}, pointer(x)), sz | x.bitmasks)
    # Bit masks
    elseif name == :free
        @setbits!(UInt8, pointer(x), UInt8(v::Bool), 0)
    elseif name == :pool
        @setbits!(UInt8, pointer(x), UInt8(v::Pool), 1, 2)
    elseif name == :dirty
        @setbits!(UInt8, pointer(x), UInt8(v::Bool), 3)
    elseif name == :evicting
        @setbits!(UInt8, pointer(x), UInt8(v::Bool), 4)
    elseif name == :queued
        @setbits!(UInt8, pointer(x), UInt8(v::Bool), 5)
    elseif name == :orphaned
        @setbits!(UInt8, pointer(x), UInt8(v::Bool), 6)

    # bytes 15..8
    elseif name == :next
        unsafe_store!(convert(Ptr{Block}, pointer(x) + 8), v::Block)

    elseif name == :id
        unsafe_store!(convert(Ptr{UInt64}, pointer(x) + 8), convert(UInt64, v))

    # bytes 23..16
    elseif name == :previous || name == :sibling
        unsafe_store!(convert(Ptr{Block}, pointer(x) + 16), v::Block)

    # bytes 31..24
    elseif name == :backsize
        unsafe_store!(convert(Ptr{UInt64}, pointer(x) + 24), convert(UInt64, v))

    else
        error("Unknown field: $name")
    end
end

isfree(x::Block) = x.free

#####
##### Metadata API
#####

setdirty!(x::Block, flag::Bool = true) = (x.dirty = flag)
isdirty(x::Block) = x.dirty
getid(x::Block) = x.id
getpool(x::Block) = x.pool

function getsibling(x::Block)
    sibling = x.sibling
    return isnull(sibling) ? nothing : sibling
end

setsibling!(x::Block, y::Block) = (x.sibling = y)
setsibling!(x::Block, ::Nothing) = (x.sibling = Block())
link!(x::Block, y::Block) = (x.sibling, y.sibling) = (y, x)
unlink!(x::Block, y::Block) = (x.sibling, y.sibling) = (Block(), Block())

@inline markqueued!(x::Block) = (x.queued = true)
@inline isqueued(x::Block) = x.queued

#####
##### Display
#####

function Base.show(io::IO, block::Block)
    if isnull(block)
        println(io, "Null Block")
        return nothing
    end

    println(io, "Block")
    println(io, "    Address: 0x", string(address(block); base = 16))
    println(io, "    Size: ", block.size)
    println(io, "    Backsize: ", block.backsize)
    print(io, "    Availability: ")

    # Display the free attributes of the block
    if isfree(block)
        printstyled(io, "Free\n"; color = :green)
        println(io, "    Next: 0x", string(address(block.next); base = 16))
        println(io, "    Previous: 0x", string(address(block.previous); base = 16))
    else
        printstyled(io, "Taken\n"; color = :red)
        println(io, "    Pool: ", block.pool)
        println(io, "    ID: ", block.id)
        println(io, "    Dirty: ", block.dirty)
        println(io, "    Evicting: ", block.evicting)
        println(io, "    Orphaned: ", block.orphaned)
        println(io, "    Sibling: 0x", string(address(block.sibling); base = 16))
    end
end


text:
using PyPlot
include("/Users/ericagol/Students/Luger/Starry/starry/julia/s2_stable.jl")

r0 = [0.01,100.]
nb = 1000
fig,axes = subplots(1,2)

epsilon = 1e-12; delta = 1e-3
for i=1:2
#  r=r0[i]
#  ax = axes[i]
#  b = [linspace(r-1,r-1+epsilon,nb); linspace(r-1+epsilon,r-1+1e-3,nb); linspace(r-1+1e-3,r-1e-3,nb); linspace(r-1e-3,r-epsilon,nb); linspace(r-epsilon,r+epsilon,nb); linspace(r+epsilon,r+1e-3,nb);  linspace(r+1e-3,r+1-1e-3,nb); linspace(r+1-1e-3,r+1-epsilon,nb); linspace(r+1-epsilon,r+1,nb)]
  ax=axes[i]
  r=r0[i]
  if r < 1.0
    b = [linspace(1e-15,epsilon,nb); linspace(epsilon,delta,nb); linspace(delta,r-delta,nb);
     r-logspace(log10(delta),log10(epsilon),nb); linspace(r-epsilon,r+epsilon,nb); r+logspace(log10(epsilon),log10(delta),nb);
     linspace(r+delta,1-r-delta,nb); 1-r-logspace(log10(delta),log10(epsilon),nb); linspace(1-r-epsilon,1-r+epsilon,nb);
     1-r+logspace(log10(epsilon),log10(delta),nb); linspace(1-r+delta,1+r-delta,nb); 1+r-logspace(log10(delta),log10(epsilon),nb);linspace(1+r-epsilon,1+r-1e-15,nb)]
     nticks = 14
     xticknames=[L"$10^{-15}$",L"$10^{-12}$",L"$10^{-3}$",L"$r-10^{-3}$",L"$r-10^{-12}$",L"$r+10^{-12}$",L"$r+10^{-3}$",
     L"$1-r-10^{-3}$",L"$1-r-10^{-12}$",L"$1-r+10^{-12}$",L"$1-r+10^{-3}$",L"$1+r-10^{-3}$",L"$1+r-10^{-12}$",L"$1+r-10^{-15}$"]
  else
    b = [r-1+logspace(log10(epsilon),log10(delta),nb); linspace(r-1+delta,r-delta,nb);
     r-logspace(log10(delta),log10(epsilon),nb); linspace(r-epsilon,r+epsilon,nb); r+logspace(log10(epsilon),log10(delta),nb);
     linspace(r+delta,r+1-delta,nb); r+1-logspace(log10(delta),log10(epsilon),nb)]
     nticks = 8
     xticknames=[L"$r-1+10^{-12}$",L"$r-1+10^{-3}$",L"$r-10^{-3}$",L"$r-10^{-12}$",L"$r+10^{-12}$",L"$r+10^{-3}$",
     L"$r+1-10^{-3}$",L"$r+1-10^{-12}$"]
  end

  s2_grid = s2.(r,b)
  s2_big = convert(Array{Float64,1},s2.(big(r),big.(b)))
  diff = s2_grid-s2_big
#  ax[:plot](b,diff)
  ax[:plot](diff)
  ax[:set_xlabel]("b")
  ax[:set_ylabel]("s2-s2(big)")
#  ax[:axis]([0,length(b),1e-16,1])
  ax[:set_xticks](nb*linspace(0,nticks-1,nticks))
  ax[:set_xticklabels](xticknames,rotation=45)
  if i==1
    ax[:set_title]("r = 0.01")
  else
    ax[:set_title]("r = 100")
  end
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark