text:
export get_dot_graph, plot_dot_graph, numericalgrad, gradcheck

function topsort(top::Var)
    sorted = Var[]
    visited = Set{Var}()
    function visit(x)
        x in visited && return
        push!(visited, x)
        x.creator isa Union{NullFunc,ParamCreator} && return push!(sorted, x)
        for v in x.creator.args
            v isa Var && visit(v)
        end
        push!(sorted, x)
    end
    visit(top)
    reverse(sorted)
end

function _dot_var(v::Var, verbose=false)
    name = getname(v)
    if verbose && isdatadefined(v)
        if name != "" name *= ": " end
        name *= string(size(v)) * " " * string(eltype(v))
    end
    return "$(objectid(v)) [label=\"$name\", color=orange, style=filled]\n"
end

function _dot_func(f::Func)
    func_name = split(string(typeof(f)), ".")[end]
    txt = "$(objectid(f)) [label=\"$(func_name)\", color=lightblue, style=filled, shape=box]\n"
    for x in f.args
        txt *= "$(objectid(x)) -> $(objectid(f))\n"
    end
    for y in f.outputs
        txt *= "$(objectid(f)) -> $(objectid(y.value))\n"
    end
    return txt
end

function get_dot_graph(output::Var; verbose=true)
    txt = ""
    variables = topsort(output)
    txt *= _dot_var(output, verbose)

    for v in variables
        f = v.creator
        f isa SymbolFunc && continue
        txt *= _dot_func(f)
        for x in f.args
            txt *= _dot_var(x, verbose)
        end
    end
    return "digraph g {\n" * txt * "}"
end

function plot_dot_graph(output::Var; verbose=false, file="graph.png")
    dot_graph = get_dot_graph(output, verbose=verbose)
    graph_path = tempname() * ".dot"
    open(graph_path, "w") do f
        write(f, dot_graph)
    end

    extension = split(file, ".")[end]
    cmd = `dot $graph_path -T $extension -o $file`
    run(cmd)
end

function allclose(a, b; rtol=1e-4, atol=1e-5)
    return all(@. abs(a - b) <= (atol + rtol * abs(b)))
end

function numericalgrad(f, x, args...; eps=1e-4)
    x = x isa Var ? x.data : x
    grad = zeros(Float64, size(x))
    for idx in 1:length(x)
        tmp = copy(x[idx])
        x[idx] = tmp + eps
        fxh1 = f(x, args...)
        x[idx] = tmp - eps
        fxh2 = f(x, args...)
        if fxh1 isa Var grad[idx] = sum(fxh1.data - fxh2.data) / (2 * eps)
        else grad[idx] = sum(fxh1 - fxh2) / (2 * eps) end
        x[idx] = tmp
    end
    return grad
end

function gradcheck(f, x, args...; rtol=1e-4, atol=1e-5)
    x = asvar(x)
    num_grad = numericalgrad(f, x, args...)
    y = f(x, args...)
    gradient!(y)
    bp_grad = x.grad.data
    size(num_grad) == size(bp_grad) || return false
    allclose(num_grad, bp_grad) || return false
    return true
end

function _sumto(x, shape)
    target_dim = length(shape)
    lead = ndims(x) - target_dim
    lead_dims = Tuple(target_dim + 1:ndims(x))
    dims = ()
    for i in 1:target_dim
        if shape[i] == 1
            dims = tuple(dims..., i + lead)
        end
    end
    y = sum(x, dims=(lead_dims..., dims...))
    lead > 0 && (y = dropdims(y, dims=lead_dims))
    return y
end


text:
using LinearAlgebra

using Yao
using Yao.Blocks
using QuAlgorithmZoo

"""
Quantum GAN.

Reference:
    Benedetti, M., Grant, E., Wossnig, L., & Severini, S. (2018).
    Adversarial quantum circuit learning for pure state approximation, 1–14.
"""
struct QuGAN{N}
    target::DefaultRegister
    generator::MatrixBlock{N}
    discriminator::MatrixBlock
    reg0::DefaultRegister
    witness_op::MatrixBlock
    circuit::AbstractBlock
    gdiffs
    ddiffs

    function QuGAN(target::DefaultRegister, gen::MatrixBlock, dis::MatrixBlock)
        N = nqubits(target)
        c = sequence(gen, addbit(1), dis)
        witness_op = put(N+1, (N+1)=>P0)
        gdiffs = collect(gen, AbstractDiff)
        ddiffs = collect(dis, AbstractDiff)
        new{N}(target, gen, dis, zero_state(N), witness_op, c, gdiffs, ddiffs)
    end
end

"""loss function"""
#loss(qcg::QuGAN) = p0t(qcg) - p0g(qcg)
"""probability to get evidense qubit 0 on generation set."""
#p0g(qg::QuGAN) = expect(qg.witness_op, psi_discgen(qg)) |> real
"""probability to get evidense qubit 0 on target set."""
#p0t(qg::QuGAN) = expect(qg.witness_op, psi_disctarget(qg)) |> real
"""generated wave function"""
psi(qg::QuGAN) = copy(qg.reg0) |> qg.generator
"""input |> generator |> discriminator"""
psi_discgen(qg::QuGAN) = copy(qg.reg0) |> qg.circuit
"""target |> discriminator"""
psi_disctarget(qg::QuGAN) = copy(qg.target) |> qg.circuit[2:end]
"""tracedistance between target and generated wave function"""
distance(qg::QuGAN) = tracedist(qg.target, psi(qg))[]

"""obtain the gradient"""
function grad(qcg::QuGAN)
    ggrad_g = opdiff.(()->psi_discgen(qcg), qcg.gdiffs, Ref(qcg.witness_op))
    dgrad_g = opdiff.(()->psi_discgen(qcg), qcg.ddiffs, Ref(qcg.witness_op))
    dgrad_t = opdiff.(()->psi_disctarget(qcg), qcg.ddiffs, Ref(qcg.witness_op))
    [-ggrad_g; dgrad_t - dgrad_g]
end

"""the training process"""
function train(qcg::QuGAN{N}, g_learning_rate::Real,
                             d_learning_rate::Real, niter::Int) where N
    ng = length(qcg.gdiffs)
    for i in 1:niter
        g = grad(qcg)
        dispatch!(+, qcg.generator, -g[1:ng]*g_learning_rate)
        dispatch!(-, qcg.discriminator, -g[ng+1:end]*d_learning_rate)
        (i*20)%niter==0 && println("Step = $i, Trance Distance = $(distance(qcg))")
    end
end

nbit = 3
target = rand_state(nbit)
gen = dispatch!(random_diff_circuit(nbit, 2,
                       pair_ring(nbit)), :random) |> autodiff(:QC)
discriminator = dispatch!(random_diff_circuit(nbit+1,
                        5, pair_ring(nbit+1)), :random) |> autodiff(:QC)
qcg = QuGAN(target, gen, discriminator)
train(qcg, 0.1, 0.2, 1000)

ext: jl


text:
@generator SWCSimpleGen begin
    start() = mult(a())
	a() = choose(Bool) ? "u" : "v"
	a() = choose(Bool) ? "x" : "y"
end

@testset "Switch choice model" begin

	@testset "constructors" begin

		gn = SWCSimpleGen()
		cm1 = GEChoiceModel(4)
		cm2 = SamplerChoiceModel(choicepointinfo(gn))

		cps = sort(collect(keys(choicepointinfo(gn))))

		cm = SwitchChoiceModel([(cps[[2,3]], cm1),], cm2)
		@test typeof(cm) <: SwitchChoiceModel

	end

	@testset "set/get parameters and ranges" begin

		gn = SWCSimpleGen()
		cm1 = GEChoiceModel(4)
		cm2 = SamplerChoiceModel(choicepointinfo(gn))
		cm3 = GEChoiceModel([1,2,3])

		cps = sort(collect(keys(choicepointinfo(gn))))
		cm = SwitchChoiceModel([(cps[[1,]], cm1), (cps[[3,4,]], cm3),], cm2)

		@test numparams(cm) == numparams(cm1) + numparams(cm3) + numparams(cm2)
		@test paramranges(cm) == [paramranges(cm1); paramranges(cm3); paramranges(cm2); ]
		@test getparams(cm) == [getparams(cm1); getparams(cm3); getparams(cm2); ]
		setparams!(cm, [5.0, 3.0, 2.0, 1.0, 10.0, 11.0, 12.0, 0.2, 0.7, 0.55, 0.4, 0.6])
		# Assumes choice point ordering of (specifically that categorical cp are tbe last two params of cm2)
	 	# 	Dict{Symbol,Any}(Pair{Symbol,Any}(:datatype,Bool),Pair{Symbol,Any}(:type,:value),Pair{Symbol,Any}(:max,true),Pair{Symbol,Any}(:min,false))
	 	#   Dict{Symbol,Any}(Pair{Symbol,Any}(:datatype,Bool),Pair{Symbol,Any}(:type,:value),Pair{Symbol,Any}(:max,true),Pair{Symbol,Any}(:min,false))
	 	#   Dict{Symbol,Any}(Pair{Symbol,Any}(:type,:sequence),Pair{Symbol,Any}(:max,9223372036854775807),Pair{Symbol,Any}(:min,0))
	 	#   Dict{Symbol,Any}(Pair{Symbol,Any}(:type,:rule),Pair{Symbol,Any}(:rulename,:a),Pair{Symbol,Any}(:max,2),Pair{Symbol,Any}(:min,1))
		@test getparams(cm) == [5.0, 3.0, 2.0, 1.0, 10.0, 11.0, 12.0, 0.2, 0.7, 0.55, 0.4, 0.6]

	end

	@testset "sampling" begin

		gn = SWCSimpleGen()
		cm1 = GEChoiceModel([3,2,1])
		cm2 = SamplerChoiceModel(choicepointinfo(gn))
		cm3 = GEChoiceModel([8,7,6])

		cps = sort(collect(keys(choicepointinfo(gn))))
		# Assumes choice point ordering of:
	 	# 	Dict{Symbol,Any}(Pair{Symbol,Any}(:datatype,Bool),Pair{Symbol,Any}(:type,:value),Pair{Symbol,Any}(:max,true),Pair{Symbol,Any}(:min,false))
	 	#   Dict{Symbol,Any}(Pair{Symbol,Any}(:datatype,Bool),Pair{Symbol,Any}(:type,:value),Pair{Symbol,Any}(:max,true),Pair{Symbol,Any}(:min,false))
	 	#   Dict{Symbol,Any}(Pair{Symbol,Any}(:type,:sequence),Pair{Symbol,Any}(:max,9223372036854775807),Pair{Symbol,Any}(:min,0))
	 	#   Dict{Symbol,Any}(Pair{Symbol,Any}(:type,:rule),Pair{Symbol,Any}(:rulename,:a),Pair{Symbol,Any}(:max,2),Pair{Symbol,Any}(:min,1))
		cm = SwitchChoiceModel([(cps[[3,]], cm1), (cps[[2,4,]], cm3),], cm2)
		setchoicemodel!(gn, cm)

		@mtestset "check sampling from all choice models and reset" begin
			d = join(choose(gn))
			@test ismatch(r"[u|v]y[u|v]", d) # multi implicitly checks reset
			@mtest_values_vary d # checks that "default" sampler choice model is applied to u/v choice
		end

	end

end


text:
# A `Table` presents itself as an `AbstractArray` of `NamedTuples`s

"""
    FlexTable(name1 = array1, ...)

Create a column-storage-based `FlexTable` with column names `name1`, etc, from arrays `array1`,
etc. The input arrays `array1`, etc, must share the same dimensionality and indices.

`FlexTable` itself is an `AbstractArray` whose elements are `NamedTuple`s of the form
`(name1 = first(array1), ...)`, etc. Rows of the table are obtained via standard array
indexing `table[i]`, and columns via `table.name`.

`FlexTable` differs from `Table` in that the columns are mutable - you may add, remove,
rename and replace entire columns of a `FlexTable`, but not a `Table`. However, `Table` can
access and iterate rows in local scope with fast, fully type-inferred code while `FlexTable`
will be more efficient with a higher-order interface.
"""
mutable struct FlexTable{N} <: AbstractArray{NamedTuple, N}
    data::NamedTuple{<:Any, <:Tuple{Vararg{AbstractArray{<:Any, N}}}}

    # Inner constructor, to compare axes?
end

FlexTable(;kwargs...) = FlexTable(kwargs.data)
FlexTable(nt::NamedTuple) = FlexTable{_ndims(nt)}(nt)

function FlexTable(x)
    cols = columns(x)
    if cols isa NamedTuple{<:Any, <:Tuple{Vararg{AbstractArray{N}} where N}}
        return FlexTable(cols)
    else
        return FlexTable(columntable(cols))
    end
end
FlexTable(t::Table) = FlexTable(columns(t))
Table(t::FlexTable) = Table(columns(t))

FlexTable{N}(t::Table{<:Any, N}) where {N} = FlexTable{N}(columns(t))

Tables.AccessStyle(::FlexTable) = Tables.ColumnAccess()
Tables.schema(t::FlexTable) = _eltypes(columns(t))

"""
    columns(dataframe::FlexTable)

Convert a `FlexTable` into a `NamedTuple` of it's columns.
"""
@inline columns(t::FlexTable) = getfield(t, :data)

@inline rows(t::FlexTable) = Table(columns(t))

# Simple column access via `table.columnname`
@inline Base.getproperty(t::FlexTable, name::Symbol) = getproperty(columns(t), name)

# Private fields are never exposed since they can conflict with column names
Base.propertynames(t::FlexTable, private::Bool=false) = columnnames(t)

function Base.setproperty!(t::FlexTable, name::Symbol, a::AbstractArray)
    setfield!(t, :data, merge(columns(t), NamedTuple{(name,)}((a,))))
    return t
end

"""
    columnnames(table)

Return a tuple of the column names of a `Table`.
"""
columnnames(t::FlexTable) = keys(columns(t))

# Basic AbstractArray interface

@inline Base.size(t::FlexTable{N}) where {N} = size(first(columns(t)))::NTuple{N, Integer}
@inline Base.axes(t::FlexTable{N}) where {N} = axes(first(columns(t)))::NTuple{N, Any}
@inline Base.IndexStyle(t::FlexTable) = IndexStyle(first(columns(t)))

function Base.checkbounds(::Type{Bool}, t::FlexTable, i...)
    # Make sure we are in bounds for *every* column. Only safe to do
    # here because each column might mutate size independently of the others!
    all(col -> checkbounds(Bool, col, i...), columns(t))
end

@inline function Base.getindex(t::FlexTable, i::Int)
    @boundscheck checkbounds(t, i)
    map(col -> @inbounds(getindex(col, i)), columns(t))::NamedTuple
end

@inline function Base.getindex(t::FlexTable, i::Int...)
    @boundscheck checkbounds(t, i...)
    map(col -> @inbounds(getindex(col, i...)), columns(t))::NamedTuple
end

@inline function Base.setindex!(t::FlexTable, v::NamedTuple, i::Int)
    @boundscheck begin
        checkbounds(t, i)
        @assert keys(v) === keys(columns(t))
    end
    map((val, col) -> @inbounds(setindex!(col, val, i)), v, columns(t))
    return t
end

@inline function Base.setindex!(t::FlexTable, v::NamedTuple, i::Int...)
    @boundscheck begin
        checkbounds(t, i)
        @assert keys(v) === keys(columns(t))
    end
    map((val, col) -> @inbounds(setindex!(col, val, i...)), v, columns(t))
    return t
end

# similar
@inline Base.similar(t::FlexTable{N}) where {N} = FlexTable{N}(similar(Table(t)))::FlexTable{N}
@inline Base.similar(t::FlexTable{N}, ::Type{NamedTuple}) where {N} = FlexTable{N}(similar(Table(t)))
@inline Base.similar(t::FlexTable{N}, ::Type{NamedTuple{names,T}}) where {N, names, T} = FlexTable{N}(similar(Table(t, NamedTuple{names, T})))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims) = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names,T}}, dims) where {names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))

# Ambiguities...
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims::Union{Integer, AbstractUnitRange}...) = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims::Tuple{Vararg{Int64,N}}) where {N} = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims::Tuple{Union{Integer, OneTo},Vararg{Union{Integer, OneTo}}}) = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))

@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names, T}}, dims::Union{Integer, AbstractUnitRange}...) where {names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names, T}}, dims::Tuple{Vararg{Int64,N}}) where {N, names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names, T}}, dims::Tuple{Union{Integer, OneTo},Vararg{Union{Integer, OneTo}}}) where {names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))

@inline __ndims(::Integer) = 1
@inline __ndims(::AbstractUnitRange) = 1
@inline __ndims(::NTuple{N, Any}) where {N} = N

# empty
Base.empty(t::FlexTable) = empty(t, _eltypes(columns(t)))
function Base.empty(t::FlexTable, ::Type{NamedTuple{names, T}}) where {names, T <: Tuple}
    FlexTable(empty(Table(t)))
end

# Support Vector / deque interface (mutable-length vectors)

function Base.empty!(t::FlexTable)
    map(empty!, columns(t))
    return t
end

function Base.pop!(t::FlexTable)
    return map(pop!, columns(t))::NamedTuple
end

function Base.push!(t::FlexTable, v::NamedTuple)
    map(push!, columns(t), v)
    return t
end

function Base.append!(t::Union{FlexTable, Table}, t2::Union{FlexTable, Table})
    map(append!, columns(t), columns(t2))
    return t
end

function Base.popfirst!(t::FlexTable)
    return map(popfirst!, columns(t))::NamedTuple
end

function Base.pushfirst!(t::FlexTable, v::NamedTuple)
    map(pushfirst!, columns(t), v)
    return t
end

function Base.prepend!(t::Union{FlexTable, Table}, t2::Union{FlexTable, Table})
    map(prepend!, columns(t), columns(t2))
    return t
end

function Base.deleteat!(t::FlexTable, i)
    map(col -> deleteat!(col, i), columns(t))
    return t
end

function Base.insert!(t::FlexTable, i::Integer, v::NamedTuple)
    map((col, val) -> insert!(col, i, val), columns(t), v)
    return t
end

function Base.splice!(t::FlexTable, inds::Integer)
    return map(col -> splice!(col, inds), columns(t))::NamedTuple
end

function Base.splice!(t::FlexTable, inds::AbstractVector)
    return FlexTable{1}(map(col -> splice!(col, inds), columns(t)))
end

function Base.splice!(t::FlexTable, inds::Integer, ins::NamedTuple)
    return map((col, vals) -> splice!(col, inds, vals), columns(t), ins)::NamedTuple
end

function Base.splice!(t::FlexTable, inds::AbstractVector, ins::NamedTuple)
    return FlexTable{1}(map((col, vals) -> splice!(col, inds, vals), columns(t), ins))
end

function Base.splice!(t::Union{FlexTable, Table}, inds::Integer, ins::Union{FlexTable, Table})
    return map((col, vals) -> splice!(col, inds, vals), columns(t), columns(ins))::NamedTuple
end

function Base.splice!(t::Union{FlexTable, Table}, inds::AbstractVector, ins::Union{FlexTable, Table})
    return FlexTable{1}(map((col, vals) -> splice!(col, inds, vals), columns(t), columns(ins)))
end

# TODO splicing in an `AbstractArray{<:NamedTuple}` should be possible...

# Speedups for column-based storage

function Base.getindex(t::FlexTable, inds::Union{AbstractArray, Colon}...)
    return FlexTable{_getindex_dims(inds)}(map(col -> getindex(col, inds...), columns(t)))
end

function Base.view(t::FlexTable, inds::Union{AbstractArray, Colon}...)
    return FlexTable{_getindex_dims(inds)}(map(col -> view(col, inds...), columns(t)))
end

@inline _getindex_dims(inds) = __getindex_dims(0, inds...)
@inline __getindex_dims(n::Int) = n
@inline __getindex_dims(n::Int, ::Int, inds...) = __getindex_dims(n, inds...)
@inline __getindex_dims(n::Int, ::AbstractArray{<:Any, m}, inds...) where {m} = __getindex_dims(n + m, inds...)
@inline __getindex_dims(n::Int, ::Colon, inds...) = __getindex_dims(n + 1, inds...)

# Deprecated for .= syntax (via Base.Broadcast.materialize!)
# It seems `Ref` might be the new cool here. Could also consider `AbstractArray{<:NamedTuple, 0}`?
#function Base.setindex!(t::Table{<:NamedTuple{names}}, v::NamedTuple{names}, inds::Union{AbstractArray, Colon}...) where {names}
#    map((col, val) -> setindex!(col, val, inds...), columns(t), v)
#    return t
#end

function Base.setindex!(t::FlexTable, t2::Union{FlexTable, Table}, inds::Union{AbstractArray, Colon}...)
    map((col, col2) -> setindex!(col, col2, inds...), columns(t), columns(t2))
    return t
end

function Base.vcat(t::Union{FlexTable, Table}, t2::Union{FlexTable, Table})
    return FlexTable{_vcat_ndims(ndims(t), ndims(t2))}(map(vcat, columns(t), columns(t2)))
end

function Base.hcat(t::Union{FlexTable, Table}, t2::Union{FlexTable, Table})
    return FlexTable{_hcat_ndims(ndims(t), ndims(t2))}(map(hcat, columns(t), columns(t2)))
end

function Base.hvcat(rows::Tuple{Vararg{Int}}, ts::Union{FlexTable, Table}...)
    return FlexTable(map((cols...,) -> hvcat(rows, cols...), map(columns, ts)...))
end

@pure function _vcat_ndims(i::Int, j::Int)
    max(i, j, 1)
end

@pure function _hcat_ndims(i::Int, j::Int)
    max(i, j, 2)
end

function Base.resize!(t::FlexTable, len::Int)
    for c in columns(t)
        resize!(c, len)
    end
    t
end


text:
using Rheologies ; const R = Rheologies

include(joinpath(@__DIR__,"sample_circ_inclusion.jl"))

output_path = @__DIR__

## SPATIAL DOMAIN ##
dim = 2 # spatial dimensions

    Lx = 1 # same as inclusion.msh mesh
    Ly = 2 # same as inclusion.msh mesh
    radius = Lx/50
    el_geom = Triangle # Quadrilateral or Triangle, equivalent to Cell{dim,nnodes,nfaces}

    # generate the grid and sets:
    meshfile = joinpath(output_path,"inclusion.msh")

    generate_mesh_inclusion(; Lx, Ly, radius ,lowres=0.04, highres=0.003, file=meshfile)
    grid = generate_grid(meshfile)
    #grid = generate_grid(el_geom, (nx, ny), corner1, Vec(corner1[1]+Lx, corner1[2]+Ly)) # can take 2 or 4 corners

    # Add facesets / nodesets / cellsets on which bc will be applied (except for the top, bottom, left and right boundary that are implemented by default)
    #addnodeset!(grid, "clamped", x -> (x[1] == corner_u[1] && 24.5<x[2]<=25.5));
    #addnodeset!(grid, "clamped", x -> ( (abs(x[1]-Lx/2)<0.01) & any(abs.(x[2] .- (0.0,Ly)).< 0.01) ) );
    addnodeset!(grid, "clamped", x -> ( (x[1] ≈ Lx/2) & any(x[2] .≈ (0.0,Ly)) ) );

    ## VARIABLES AND INTERPOLATIONS
    variables = PrimitiveVariables{1}((:u,), (2,), el_geom)
    #variables = PrimitiveVariables{2}((:u,:p), (2,1), el_geom) # not good !!!!

    quad_order = 2
    quad_type = :legendre # or :lobatto

    ## BOUNDARY CONDITIONS ##
    # DIRICHLET :
    # KEY : set on which to apply bc
    # VALUES :
    # 1) constrained variable
    # 2) constraint function of (x,t)
    # 3) constrained components
    # NEUMANN :
    # KEY : set on which to apply bc
    # VALUES : traction as a function of (x,t)
    v_in = 1e-5
    bc = BoundaryConditions(dirichlet = Dict("top" => [:u, (x,t)-> -v_in*t, 2],
                                             "bottom" => [:u, (x,t)-> v_in*t, 2],
                                             "clamped" => [:u, (x,t)-> 0.0, 1]),
                            neumann   = nothing)#Dict("left" => (x,t)->Vec(1e3,0.0),"right" => (x,t)->Vec(-1e3,0.0)))

    bf = BodyForces([0.0,0.0])#-9.81*2700000

    ## MATERIAL RHEOLOGY :
    seed(x,x0,r,val_in::T,val_out::T) where {T<:Real} = (sqrt((x[1]-x0[1])^2 + (x[2]-x0[2])^2) <= r ?
                                       val_in : val_out)

    # Define seed location accordingly to the mesh seed
    x0 = Vec(Lx/2,Ly/2)

    elas = Elasticity(E = x -> seed(x,x0,radius,30e9,70e9),
                      ν = x -> seed(x,x0,radius,0.45,0.3) )

    Δσ = 5e4 ; ϵ̇ᵖ_estimate = 2v_in/(Ly/40) # to

    plas  = DruckerPrager(ϕ = 30.0, # friction angle
                          C = x -> seed(x,x0,radius,1e6,2e6), # cohesion
                          H = 0.0, # hardening modulus
                          ηᵛᵖ = Δσ/ϵ̇ᵖ_estimate) # viscosity of the visco-plastic element

    rheology = Rheology(elasticity = elas,
                        plasticity = plas )

    ### CLOCK ###
    clock = Clock(tspan = (0.0,15.0), Δt = 1.0, Δt_max = 10.0, Δt_fact_up = 1.2)

    ### SOLVER ###
    linsolver = BackslashSolver()
    #linsolver = MUMPS()
    # linsolver = ConjugateGradient(package = :IterativeSolvers,
    #                              (max_iter = 1000,))
    #linsolver = MINRES(package = :IterativeSolvers)
    nlsolver = NewtonRaphson(max_iter_number = 15, atol = 1e-5, linear_solver = linsolver)

    ### MODEL ###
    model = Model( grid = grid,
                   variables = variables,
                   quad_order = quad_order,
                   quad_type = quad_type,
                   bc = bc,
                   body_forces = bf,
                   rheology = rheology,
                   initial_state = nothing,
                   clock = clock,
                   solver = nlsolver )

    ### OUTPUT ###
    outputs = Dict(:σxx     => (r,s)-> s.σ[1,1],
                   :σyy     => (r,s)-> s.σ[2,2],
                   :σzz     => (r,s)-> s.σ[3,3],
                   :ϵxx     => (r,s)-> s.ϵ[1,1],
                   :ϵyy     => (r,s)-> s.ϵ[2,2],
                   :ϵzz     => (r,s)-> s.ϵ[3,3],
                   :ϵ_vol   => (r,s)-> tr(s.ϵ),
                   :gamma   => (r,s)-> sqrt(2* dev(s.ϵ) ⊡ dev(s.ϵ)),
                   :E       => (r,s)-> r.elasticity.E,
                   :nu      => (r,s)-> r.elasticity.ν,
                   :norm_ep => (r,s)-> norm(s.ϵᵖ),
                   :τ       => (r,s)-> R.get_τ(dev(s.σ),r.plasticity),
                   :p       => (r,s)-> -1/3 * tr(s.σ),
                   :τoverp  => (r,s)-> R.get_τ(dev(s.σ),r.plasticity)/(-1/3 * tr(s.σ))
                   )
    # mat_outputs = Dict(:sigma_xx     => (r,s)-> s.σ[1,1],
    #               :sigma_yy     => (r,s)-> s.σ[2,2],
    #               :sigma_zz     => (r,s)-> s.σ[3,3],
    #               :eps_xx     => (r,s)-> s.ϵ[1,1],
    #               :eps_yy     => (r,s)-> s.ϵ[2,2],
    #               :eps_zz     => (r,s)-> s.ϵ[3,3],
    #               :norm_ep    => (r,s)-> norm(s.ϵᵖ),
    #               :tau        => (r,s)-> R.get_τ(dev(s.σ),r.plasticity),
    #               :p          => (r,s)-> -1/3 * tr(s.σ),
    #               :tauoverp  => (r,s)-> R.get_τ(dev(s.σ),r.plasticity)/-(1/3 * tr(s.σ))
    #               )

    VTK_ow = VTKOutputWriter(model, joinpath(output_path,"inclusion_elastique_gmsh_EP_with_eta_vp_$(plas.ηᵛᵖ)"), outputs, interval = 1, force_path = true);
    #MAT_ow = MATOutputWriter(model, joinpath(path,"sawcut_EP_MAT"), mat_outputs, interval = 1, force_path = true);
    #JLD2_ow = JLD2OutputWriter(model, joinpath(path,"seed_EP_JLD2"), outputs, interval = 1, force_path = true);
    #ow = MixedOutputWriter(VTK_ow,JLD2_ow);
### SOLVE ###
    @time modelsol, u = solve(model, output_writer = VTK_ow ,log = true);


text:
module QuadraturesTests

using Gridap.Fields
using Gridap.Integration

coords = Point{2,Float64}[(0.0,0.0),(0.5,0.5),(1.0,1.0)]
weights = [0.5, 1.0, 0.5]
quad = GenericQuadrature(coords,weights)
test_quadrature(quad)

end # module


text:
################################################################################
#
#	CONCRETE TYPE
#
#   BrillouinZone{R}
#   --> R is the reciprocal unitcell type (<: AbstractReciprocalUnitcell{P,B})
#       --> P is the reciprocal site type (<: AbstractReciprocalPoint{D})
#           --> D is the dimension of embedding space
#       --> B is the bond type (<: AbstractBond{L,N})
#           --> N is the dimension of the Bravais lattice that the bond is located in
#           --> L is the label type of bonds
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface implementation
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
mutable struct BrillouinZone{R} <: AbstractBrillouinZone{R}

    # reciprocal unitcell
    reciprocal_unitcell :: R

    # corners
    corners :: Vector{Vector{Float64}}

    # faces
    faces :: Vector{Vector{Int64}}

end

# export the type
export BrillouinZone

################################################################################
#
#	INTERFACING / ACCESSING BRILLOUIN ZONES
#	(functions that have to be overwritten by concrete types)
#
################################################################################

# constructor
function newBrillouinZone(
            ::Type{BrillouinZone{RU}},
            reciprocal_unitcell :: RU,
            corners             :: Vector{<:Vector{<:Real}},
            faces               :: Vector{<:Vector{<:Integer}}
        ) :: BrillouinZone{RU} where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},RU<:AbstractReciprocalUnitcell{P,B}}

    # call the constructor
    return BrillouinZone{RU}(reciprocal_unitcell, corners, faces)
end

# export constructor
export newBrillouinZone

# accessing the corner points of the BZ
function corners(
            bz :: BrillouinZone{RU}
        ) :: Vector{Vector{Float64}} where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},RU<:AbstractReciprocalUnitcell{P,B}}

    # return the corners
    return bz.corners
end

# setting the corner points of the BZ
function corners!(
            bz      :: BrillouinZone{RU},
            corners :: Vector{<:Vector{<:Real}}
        ) where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},RU<:AbstractReciprocalUnitcell{P,B}}

    # set the corners
    bz.corners = corners
end

# export corner interface
export corners, corners!

# accessing the faces of the BZ
function faces(
            bz :: BrillouinZone{RU}
        ) :: Vector{Vector{Int64}} where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},RU<:AbstractReciprocalUnitcell{P,B}}

    # return the faces
    return bz.faces
end

# setting the faces of the BZ
function faces!(
            bz    :: BrillouinZone{RU},
            faces :: Vector{<:Vector{<:Integer}}
        ) where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},RU<:AbstractReciprocalUnitcell{P,B}}

    # set the faces
    bz.faces = faces
end

# export faces interface
export faces, faces!

# accessing the reciprocal unitcell of the BZ
function reciprocalUnitcell(
            bz :: BrillouinZone{RU}
        ) :: RU where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},RU<:AbstractReciprocalUnitcell{P,B}}

    # return the reciprocal unitcell
    return bz.reciprocal_unitcell
end

# setting the reciprocal unitcell of the BZ
function reciprocalUnitcell!(
            bz  :: BrillouinZone{RU},
            ruc :: RU
        ) where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},RU<:AbstractReciprocalUnitcell{P,B}}

    # set the reciprocal unitcell
    bz.reciprocal_unitcell = ruc
end

# export reciprocalUnitcell interface
export reciprocalUnitcell, reciprocalUnitcell!


text:
# Autogenerated wrapper script for SuiteSparse_LBT_jll for i686-w64-mingw32
export libamd, libbtf, libcamd, libccolamd, libcholmod, libcolamd, libklu, libldl, librbio, libspqr, libsuitesparse_wrapper, libsuitesparseconfig, libumfpack

using libblastrampoline_jll
JLLWrappers.@generate_wrapper_header("SuiteSparse_LBT")
JLLWrappers.@declare_library_product(libamd, "libamd.dll")
JLLWrappers.@declare_library_product(libbtf, "libbtf.dll")
JLLWrappers.@declare_library_product(libcamd, "libcamd.dll")
JLLWrappers.@declare_library_product(libccolamd, "libccolamd.dll")
JLLWrappers.@declare_library_product(libcholmod, "libcholmod.dll")
JLLWrappers.@declare_library_product(libcolamd, "libcolamd.dll")
JLLWrappers.@declare_library_product(libklu, "libklu.dll")
JLLWrappers.@declare_library_product(libldl, "libldl.dll")
JLLWrappers.@declare_library_product(librbio, "librbio.dll")
JLLWrappers.@declare_library_product(libspqr, "libspqr.dll")
JLLWrappers.@declare_library_product(libsuitesparse_wrapper, "libsuitesparse_wrapper.dll")
JLLWrappers.@declare_library_product(libsuitesparseconfig, "libsuitesparseconfig.dll")
JLLWrappers.@declare_library_product(libumfpack, "libumfpack.dll")
function __init__()
    JLLWrappers.@generate_init_header(libblastrampoline_jll)
    JLLWrappers.@init_library_product(
        libamd,
        "bin\\libamd.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbtf,
        "bin\\libbtf.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcamd,
        "bin\\libcamd.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libccolamd,
        "bin\\libccolamd.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcholmod,
        "bin\\libcholmod.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolamd,
        "bin\\libcolamd.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libklu,
        "bin\\libklu.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libldl,
        "bin\\libldl.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librbio,
        "bin\\librbio.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspqr,
        "bin\\libspqr.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsuitesparse_wrapper,
        "bin\\libsuitesparse_wrapper.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsuitesparseconfig,
        "bin\\libsuitesparseconfig.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libumfpack,
        "bin\\libumfpack.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const USE_GPU = false
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, Statistics
# pyplot()
do_viz = false

@parallel function compute_flux!(qx::Data.Array, qy::Data.Array, T::Data.Array, λ::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(qx) = -λ*@d_xi(T)/dx
    @all(qy) = -λ*@d_yi(T)/dy
    return
end

@parallel function update_T!(T::Data.Array, qx::Data.Array, qy::Data.Array, dt::Data.Number, ρCp::Data.Number, dx::Data.Number, dy::Data.Number)
    @inn(T) = @inn(T) - dt/ρCp*(@d_xa(qx)/dx + @d_ya(qy)/dy)
    return
end

@views function heat_2D_xpu()
    # Physics
    lx, ly = 10.0, 10.0
    λ      = 1.0
    ρCp    = 1.0
    nt     = 200
    # Numerics
    nx, ny = 128, 128
    nout   = 10
    # Derived numerics
    dx, dy = lx/nx, ly/ny
    dt     = min(dx^2,dy^2)/ρCp/λ/4.1
    xc     = LinRange(dx/2, lx-dx/2, nx)
    yc     = LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    qx     = @zeros(nx-1,ny-2)
    qy     = @zeros(nx-2,ny-1)
    # Initial condition
    T      = Data.Array( exp.(.-(xc.-lx./2.0).^2 .-(yc.-ly./2.0)'.^2) )
    # Time loop
    for it = 1:nt
        if (it==11) global t0 = Base.time() end
        @parallel compute_flux!(qx, qy, T, λ, dx, dy)
        @parallel update_T!(T, qx, qy, dt, ρCp, dx, dy)
        if mod(it,nout)==0 && do_viz
            display(heatmap(xc, yc, Array(T)', xlabel="lx", ylabel="ly", title="heat diffusion, it=$it", clims=(0.,1.)))
            # sleep(.01)
        end
    end
    time_s = (Base.time()-t0)
    @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", time_s, round((2/1e9*nx*ny*sizeof(lx))/(time_s/(nt-10)), sigdigits=2))
    return
end

@time heat_2D_xpu()

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text: include("zwarpedarray.jl")


text:
module RBC_3ways

greet() = print("Hello World!")

Perturbation() =  Execute RBC_Perturbation.jl

end # module


text:
using Asteroids: main

main(keepalive=!Base.isinteractive())


text:
using Test, CSV, ibm2ieee

# ibm32 requires two checks
check(ibm::UInt32, ieee::Float64) = t32(ibm, ieee) && t64(ibm, ieee)
t32(ibm::UInt32, ieee::Float64) = reinterpret(Float32, ibm2ieee.ibm32ieee32(ibm)) == Float32(ieee)

# we can't verify guarantee equality because ibm32=>ieee32 is lossy
# but we can validate against the lossage
t64(ibm::UInt32, ieee::Float64) = Float32(reinterpret(Float64, ibm2ieee.ibm32ieee64(ibm))) == Float32(ieee)

# ibm64 only requires one check.
check(ibm::UInt64, ieee::Float64) = reinterpret(Float64, ibm2ieee.ibm64ieee64(ibm)) == ieee

# examples from wikipedia and other random sources verified by Python version
@test check(0xc276a000, -118.625)
@test check(0xC2808000, -128.5)
@test check(0x40600000, 0.375)
@test check(0xBEC80000, -0.0030517578)
@test check(0x42808000, 128.5)
@test check(0x427B7333, 123.45)

@test check(0x427B733333333333, 123.45)
@test check(0x4280000000000000, 128.0)
@test check(0x4060000000000000, 0.375)

# now check 10,000 randomly chosen bit patterns
@test begin
    ibm = CSV.File("input.csv", type=UInt32).ibm
    ieee32 = CSV.File("output.csv").ieee
    ieee64 = CSV.File("output64.csv").ieee
    u32 = reinterpret.(Float32, ibm2ieee.ibm32ieee32.(ibm))
    u64 = reinterpret.(Float64, ibm2ieee.ibm32ieee64.(ibm))
    return all(isinf.(ieee32[isinf.(u32)])) && all(u32[.!isinf.(u32)] .== ieee32[.!isinf.(u32)]) && all(u64 == ieee64)
end


text:
function generate_switched_linear_systems(numMode,dim)
    A = [randn(dim,dim) for i=1:numMode]
    for i = 1:numMode
        A[i] = A[i]/maximum(abs.(eigvals(A[i])))#opnorm(A[i])
    end
    return A
end

function generate_trajectories(horizon,A,numTraj)

numMode = size(A)[1]
dim = size(A[1])[1]

x0 = zeros(dim,numTraj)
x = zeros(dim,numTraj)
for i=1:numTraj
    y0 = randn(dim)
    x0[:,i] = y0/norm(y0)
    x[:,i] = x0[:,i]
    b = rand((1:numMode),horizon)
    for j = 1:horizon
        x[:,i] = A[b[j]]*x[:,i]
    end
end
return x0, x
end

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
include("utils.jl")

""" Adaptive facility location model, with JuMP using robust counterpart,
    under budget uncertainty.
    Note: We restrict u and V to be elementwise positive, for tractability! """
function adaptive_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = Model(Gurobi.Optimizer)

    # VARIABLES
    @variable(model, x[1:n], Bin)           # Facility locations
    @variable(model, V[1:n, 1:m, 1:m] >= 0) # Affine policy w.r.t each z
    @variable(model, u[1:n, 1:m] >= 0)

    # CONSTRAINTS
    for i = 1:n # Positivity constraints
        for j = 1:m
            infdummy = @variable(model)
            normdummy = @variable(model, [1:m])
            y = @variable(model, [1:m])
            @constraint(model, [k=1:m], normdummy[k] >= y[k])
            @constraint(model, [k=1:m], normdummy[k] >= -y[k])
            @constraint(model, [k=1:m], infdummy >= -y[k] + V[i,j,k])
            @constraint(model, [k=1:m], infdummy >= y[k] - V[i,j,k])
            @constraint(model, u[i,j] >= rho*sum(normdummy) + Gamm * infdummy)
        end
    end
    for j=1:m # Demand constraints
        infdummy = @variable(model)
        normdummy = @variable(model, [1:m])
        y = @variable(model, [1:m])
        @constraint(model, [i=1:m], normdummy[i] >= y[i])
        @constraint(model, [i=1:m], normdummy[i] >= -y[i])
        dt = zeros(n, m);
        [dt[i,j] = 1 for i=1:n];
        @constraint(model, [i=1:m], infdummy >= -y[i] - P[i,j] + sum(V[:,:,i].*dt))
        @constraint(model, [i=1:m], infdummy >= y[i] + P[i,j] - sum(V[:,:,i].*dt))
        @constraint(model, sum(dt .* u) >= d[j] + rho*sum(normdummy) + Gamm * infdummy)
    end
    for i=1:n # Capacity constraints
        infdummy = @variable(model)
        normdummy = @variable(model, [1:m])
        y = @variable(model, [1:m])
        @constraint(model, [j=1:m], normdummy[j] >= y[j])
        @constraint(model, [j=1:m], normdummy[j] >= -y[j])
        dt = zeros(n, m);
        [dt[i,j] = 1 for j=1:m];
        @constraint(model, [j=1:m], infdummy >= -y[j] + sum(V[:,:,j].*dt))
        @constraint(model, [j=1:m], infdummy >= y[j] - sum(V[:,:,j].*dt))
        @constraint(model, sum(dt .* u) + rho*sum(normdummy) + Gamm * infdummy <=
                    s[i] * x[i]) # capacity constraint
    end
    # OBJECTIVE
    @variable(model, F)
    infdummy = @variable(model)
    normdummy = @variable(model, [1:m])
    y = @variable(model, [1:m])
    @constraint(model, [j=1:m], normdummy[j] >= y[j])
    @constraint(model, [j=1:m], normdummy[j] >= -y[j])
    @constraint(model, [j=1:m], infdummy >= - y[j] + sum(c.* V[:,:,j]))
    @constraint(model, [j=1:m], infdummy >= y[j] - sum(c.* V[:,:,j]))
    @constraint(model, F >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * u[i, j] for i=1:n, j=1:m) +
                                rho*sum(normdummy) + Gamm * infdummy)
    @objective(model, Min, F)
    return model, x, u, V
end

# Adaptive model through robust counterpart
rho = 1
Gamm = 5
model, x, u, V = adaptive_facility_model(c, f, rho, Gamm)
@time optimize!(model)
plt = plot_solution(model, x, u)
@show plt


text:
function demoBurglarSimple()
    # This demo is the same as demoBurglar.jl, except that here we do not make use of the available potential indexing methods and also a potential dictionary (as in demoBurglarDictPot).

    burglar,earthquake,alarm,radio=1,2,3,4 # variables must start from 1 and be contiguous
    yes, no= 1,2 # define states, starting from 1. These do need to be contiguous, starting from 1, since currently the probability tables are defined by arrays of contiguous states.

    # This defines the meaning of the variables and their states.
    # It is not necessary for the computation, but is useful for
    # printing the results in a meaningful form:
    DictVariable=Dict{Integer,DiscreteVariable}()
    DictVariable[burglar]=DiscreteVariable("burglar",["yes","no"])
    DictVariable[earthquake]=DiscreteVariable("earthquake",["yes","no"])
    DictVariable[alarm]=DiscreteVariable("alarm",["yes","no"])
    DictVariable[radio]=DiscreteVariable("radio",["yes","no"])

    # Define the discrete probability tables:
    PotBurglar=PotArray(burglar,[0.01 0.99])
    PotEarthquake=PotArray(earthquake,[0.000001 1-0.000001])

    table=zeros(2,2)
    table[yes,yes]=1
    table[no,yes]=0
    table[yes,no]=0
    table[no,no]=1
    PotRadio=PotArray([radio earthquake],table)

    table=zeros(2,2,2)
    table[yes,yes,yes]=0.9999
    table[yes,yes,no]=0.99
    table[yes,no,yes]=0.99
    table[yes,no,no]=0.0001
    table[no,:,:]=1-table[yes,:,:]
    PotAlarm=PotArray([alarm burglar earthquake],table)

    jointpot=prod([PotBurglar PotRadio PotAlarm PotEarthquake])

    println("p(burglar|alarm=yes):")
    show(condpot(setpot(jointpot,alarm,yes),burglar),DictVariable)

    println("p(burglar|alarm=yes,radio=yes):")
    show(condpot(setpot(jointpot,[alarm radio],[yes yes]),burglar),DictVariable)

end


text:
using HDF5
using ProgressMeter

@generated function write_solution_to_hdf5(op::CollisionOperator{M,N,ℳ,𝒩,RT,CT}, u::Matrix{RT}, n::Int, h5ϕ::HDF5.Dataset, h5ω::HDF5.Dataset) where {M,N,ℳ,𝒩,RT,CT}
    local û::Matrix{CT} = zeros(CT,ℳ,𝒩)
    local ϕ̂::Matrix{CT} = zeros(CT,ℳ,𝒩)
    local ϕ::Matrix{RT} = zeros(RT,M,N)

    quote
        prfft!(op.ft, u, $û)
        $ϕ̂ .= op.Δ⁻¹ .* $û
        irfft!(op.ft, $ϕ̂, $ϕ)

        h5ϕ[:,:,n] = $ϕ
        h5ω[:,:,n] = u
    end
end

"""
Run simulation for nt time steps Δt starting from u₀ and write output to HDF5 file <output>.
"""
function run_simulation(op::CollisionOperator{M,N,ℳ,𝒩,RT,CT}, nt::Int, Δt::RT, u₀::Matrix{RT}, output::String, nsave::Int=1) where {M,N,ℳ,𝒩,RT,CT}

    # create HDF5 output file
    h5 = h5open(output, "w")

    h5ϕ = create_dataset(h5, "ϕ", RT, ((M,N,nt+1), (M,N,-1)), chunk=(M,N,1))
    h5ω = create_dataset(h5, "ω", RT, ((M,N,nt+1), (M,N,-1)), chunk=(M,N,1))

    # TODO Save grid and time steps to HDF5.

    # write initial conditions to HDF5 file
    write_solution_to_hdf5(op, u₀, 1, h5ϕ, h5ω)

    # run for nt time steps
    u₁ = zero(u₀)

    @showprogress 1 for n in 1:nt
        timestep!(op, u₀, u₁, Δt)
        mod(n, nsave) == 0 || n == nt ? write_solution_to_hdf5(op, u₁, n+1, h5ϕ, h5ω) : nothing
        u₀ .= u₁

        mod(n, 10*nsave) == 0 ? flush(h5) : nothing
    end

    # close HDF5 output file
    close(h5)

end

"""
Initialise u₀ from uinit(x,y) function and run simulation.
"""
function run_simulation(op::CollisionOperator, nt::Int, Δt::Number, uinit::Function, output::String, nsave::Int=1)
    u₀ = get_field(op.grid)

    for j in 1:size(u₀,2)
        for i in 1:size(u₀,1)
             u₀[i,j]  = uinit(op.grid.x[i], op.grid.y[j])
        end
    end

    run_simulation(op, nt, Δt, u₀, output, nsave)
end

"""
Initialise u₀ from HDF5 file <input> and run simulation.
"""
function run_simulation(op::CollisionOperator, nt::Int, Δt::Number, input::String, output::String, nsave::Int=1)
    u₀ = get_field(op.grid)

    h5open(input, "r") do h5
        @assert size(h5.ω,1) == size(op.grid,1)
        @assert size(h5.ω,2) == size(op.grid,2)
        u₀ .= view(h5["ω"][:,:,end], :, :, 1)
    end

    run_simulation(op, nt, Δt, u₀, output, nsave)
end


text:
function plotreturn(x, title::String,  n::Int, col::String)

#  gnuplot_send("set style fill solid border -1") # Segmentation fault: 11

  histogram(x,
             "title", title,
             "bins", n,
             "color", col,
             "xlabel", "returns",
             "ylabel", "number of occurenaces",
             "legend", "legend"
             )
end
plotreturn(x) = plotreturn(x, "", 10, "blue")

function plotseries(x::DataFrame, col::String)

  plot([1:nrow(x)], float(x[:,col]))

end


text:
using CoordinateTransformations

@testset "Datum transforms" begin

@testset "GDA94 from ITRF" begin
    # Tests for ITRF -> GDA94 conversion.  Unfortunately there doesn't seem to be a
    # public testset for this.

    @testset "Basic example" begin
        # Basic worked example from Dawson and Woods
        x_ITRF  = ECEF(-4052052.3678, 4212836.0411, -2545105.1089)
        x_GDA94 = ECEF(-4052051.7615, 4212836.1945, -2545106.0145)
        @test norm(x_GDA94 - Geodesy.GDA94_from_ITRF_Dawson2010(2005, Date(2010,6,16))(x_ITRF)) < 2e-4

        # Test inverse
        epoch = Date(2010,6,16)
        trans = datum_shift_ECEF(ITRF{2008}(epoch), GDA94) ∘ datum_shift_ECEF(GDA94, ITRF{2008}(epoch))
        @test trans.m ≈ eye(3)
        @test isapprox(trans.v, zeros(3), atol=10*eps())
    end

    @testset "ITRS realizations" begin
        x_GDA94 = ECEF(-4052051.7615, 4212836.1945, -2545106.0145)
        epoch = Date(2010,6,16)

        # Test dispatch
        @test datum_shift_ECEF(GDA94(), ITRF{2008}(epoch))(x_GDA94) == Geodesy.GDA94_from_ITRF_Dawson2010(2008, epoch)(x_GDA94)
        @test datum_shift_ECEF(GDA94(), ITRF{1996}(epoch))(x_GDA94) == Geodesy.GDA94_from_ITRF_Dawson2010(1996, epoch)(x_GDA94)

        # Nonexistent ITRF realization
        @test_throws Exception datum_shift_ECEF(ITRF{10000}(epoch), GDA94())
    end

    @testset "GDA94 from ITRF" begin
        # The following test data was obtained by submitting a GPS processing
        # request to the Geosciences Australia AUSPOS service, extracting the
        # cartesian coordinates of the base stations used in the solution from
        # the pdf report.

        # Cartesian, GDA94
        stations_GDA94 = [
            # Station X (m)         Y (m)        Z (m)
            #=BDST=# ECEF(-5021920.612, 2559339.872, -2975290.670),
            #=BNDY=# ECEF(-5125976.801, 2688801.593, -2669891.529),
            #=CBLT=# ECEF(-5061144.438, 2584178.830, -2886586.871),
            #=CNDO=# ECEF(-4494145.811, 2901721.598, -3462001.509),
            #=CWRA=# ECEF(-4532143.250, 2755335.963, -3531098.620),
            #=DALB=# ECEF(-4979266.855, 2730160.233, -2895220.922),
            #=GATT=# ECEF(-5012218.430, 2628002.744, -2931853.186),
            #=GONG=# ECEF(-4601821.645, 2561462.987, -3585678.678),
            #=IPS2=# ECEF(-5028440.649, 2588779.894, -2938802.549),
            #=MSVL=# ECEF(-4571853.513, 2599983.302, -3597310.167),
            #=PTKL=# ECEF(-4599805.569, 2558779.075, -3590076.333),
            #=ROBI=# ECEF(-5034843.827, 2523322.876, -2984064.629),
            #=RSBY=# ECEF(-5121088.692, 2863243.621, -2493144.670),
            #=TOOW=# ECEF(-4994481.724, 2663618.280, -2931172.448),
            #=WARW=# ECEF(-4967998.840, 2638152.833, -2997613.670),
        ]

        # Cartesian, ITRF2008
        ITRF_instance = ITRF{2008}(Date(2012,05,08))
        stations_ITRF2008 = [
            # Station X (m)         Y (m)       Z (m)
            #=BDST=# ECEF(-5021921.133, 2559339.665, -2975289.754),
            #=BNDY=# ECEF(-5125977.300, 2688801.332, -2669890.587),
            #=CBLT=# ECEF(-5061144.951, 2584178.607, -2886585.948),
            #=CNDO=# ECEF(-4494146.434, 2901721.514, -3462000.602),
            #=CWRA=# ECEF(-4532143.863, 2755335.886, -3531097.728),
            #=DALB=# ECEF(-4979267.387, 2730160.019, -2895219.988),
            #=GATT=# ECEF(-5012218.954, 2628002.532, -2931852.262),
            #=GONG=# ECEF(-4601822.242, 2561462.912, -3585677.803),
            #=IPS2=# ECEF(-5028441.169, 2588779.681, -2938801.629),
            #=MSVL=# ECEF(-4571854.116, 2599983.231, -3597309.290),
            #=PTKL=# ECEF(-4599806.166, 2558779.000, -3590075.459),
            #=ROBI=# ECEF(-5034844.345, 2523322.669, -2984063.717),
            #=RSBY=# ECEF(-5121089.191, 2863243.334, -2493143.707),
            #=TOOW=# ECEF(-4994482.252, 2663618.069, -2931171.521),
            #=WARW=# ECEF(-4967999.373, 2638152.635, -2997612.749),
        ]

        to_GDA = datum_shift_ECEF(GDA94(), ITRF_instance)
        # Need a nonzero tolerance which is expected, but it should be half
        # this size if (1) GA correctly rounded the coordinates (2) we used
        # exactly the same time computation and (3) list of transformation
        # coefficients.
        tolerance = 1.1e-3

        for i = 1:size(stations_ITRF2008,1)
            x_ITRF = stations_ITRF2008[i]
            x_GDA  = stations_GDA94[i]
            @test norm(x_GDA - to_GDA(x_ITRF)) < tolerance
        end
    end

end

end # @testset


text:
function value(obj::UnivariateVectorFunction, arg0::jdouble)
    return jcall(obj, "value", Vector{jdouble}, (jdouble,), arg0)
end


text:
# # `interfaces.jl`
#
# Definitions of secondary interfaces.
#
# ## Inplace interface

encode!(buf, method, context, sample) = encode(method, context, sample)
encodeinput!(buf, method, context, input) = encodeinput(method, context, input)
encodetarget!(buf, method, context, target) = encodetarget(method, context, target)
decodey!(buf, method, context, y) = decodey(method, context, y)
decodeŷ!(buf, method, context, ŷ) = decodeŷ(method, context, ŷ)

# ## Interpretation interface

function interpretsample end
function interpretinput end
function interprettarget end
function interpretx end
function interpretŷ end
function interprety end

# ## Training interface

"""
    methodlossfn(method)

Default loss function to use when training models for `method`.
"""
function methodlossfn end

"""
    methodmodel(method, backbone)

Construct a model for `method` from a backbone architecture, for example
by attaching a method-specific head model.
"""
function methodmodel end

# ## Testing interface

"""
    mocksample(method)

Generate a random `sample` compatible with `method`.
"""
mocksample(method) = (mockinput(method), mocktarget(method))

"""
    mockinput(method)

Generate a random `input` compatible with `method`.
"""
function mockinput end

"""
    mocktarget(method)

Generate a random `target` compatible with `method`.
"""
function mocktarget end

"""
    mockmodel(method)

Generate a random `model` compatible with `method`.
"""
function mockmodel end


text:
using TailRiskScenGen
using Test

println("Running test_monotonic_utils.jl...")
include("test_monotonic_utils.jl")

println("Running test_RiskRegion.jl...")
include("test_RiskRegion.jl")

println("Running test_BruteForceRiskRegion.jl...")
include("test_BruteForceRiskRegion.jl")

println("Running test_scen_gen.jl...")
include("test_scen_gen.jl")

println("Done tests")

## # Parameters
## mean = [0.1, 0.05]
## cov = [[1.0, -0.5] [-0.5, 1.5]]
## num_scen = 100
## beta = 0.95

## dim = length(mean)
## mv_norm = MvNormal(mean, cov)
## scenarios = rand(mv_norm, num_scen)
## probs = fill(1.0/num_scen, num_scen)
## alpha = quantile(Normal(), 0.95)
## cone_A = eye(dim)

## print("Running aggregate_scenarios with test arguments...")
## agg_scen, agg_prob = aggregate_scenarios(scenarios, probs, mean, cov, cone_A, alpha)
## print(agg_scen)
## print(agg_prob)
## print("\n")


text:
module TestDataFrames

  using GBIF
  using DataFrames
  using Test

  df = DataFrame(occurrences())

  @test typeof(df) <: DataFrame

end


text:
__precompile__()

module JsonCall

using FunctionalData, JSON, HttpServer, Requests, HDF5

export jsoncall, jsonhandler, serve, savehdf5, loadhdf5

const HDF5EXT = ".hdf5"
const SINGLEVARIABLE = "jsoncall_singlevariable"

######################
##  parsejson

function parsejson(a::AbstractString)
    r = JSON.parse(a)
    @p parsejson_ r
end

parsejson_(a) = a

function parsejson_(a::AbstractString)
    r = endswith(a, HDF5EXT) ? loadhdf5(a) : a
end

function parsejson_(a::Dict)
    if length(collect(keys(a))) == 2 && haskey(a, "arraysize") && haskey(a,"arraydata")
        return typed(reshape(a["arraydata"], a["arraysize"]...))
    else
        @p mapvalues a parsejson_
    end
end

######################
##  makejson

makejson(a) = JSON.json(makejson_(a))
makejson_(a) = a
makejson_(a::Array) = try makejson_(typed(a)) catch a end
function makejson_{T<:Number,N}(a::Array{T,N})
    Dict("arraysize" => size(a), "arraydata" => vec(a))
end
makejson_(a::Dict) = @p mapvalues a makejson_

######################
##  jsonhandler

function jsonhandler(f)
    jsonhandler = HttpHandler() do req::Request, res::Response
        d = parsejson(ascii(req.data))
        r = f(d)
        Response(makejson(r))
    end
end

######################
##  jsoncall

function jsoncall(d::Dict, port = 8000, url = "http://localhost:$port")
    r = readall(post("http://localhost:$port"; data = makejson(d)))
    parsejson(r)
end

######################
##  jsoncall

serve(f, port = 8000) = @async @p jsonhandler f | Server | run port

######################
##  savehdf5

savehdf5(filename, data) = savehdf5(filename, Dict(SINGLEVARIABLE => data))
function savehdf5(filename, data::Dict)
    endswith(filename, HDF5EXT) || error("JsonCall: Filename must end in '$HDF5EXT'")
    h5open(filename, "w") do file
        for key in keys(data)
            write(file, key, data[key])
        end
    end
    filename
end

######################
##  loadhdf5

function loadhdf5(filename)
    d = Dict()
    h5open(filename, "r") do file
        for key in names(file)
            d[key] = read(file, key)
        end
    end
    isa(d, Dict) && haskey(d, SINGLEVARIABLE) ? d[SINGLEVARIABLE] : d
end

end # module


text:
using LinearAlgebra
using RecipesBase

struct Doggo{T}
    leg::DoggoLeg{T}
    T_leg::Vector{AffineMap{Quat{T}, Vector{Int}}}
    n_legs::Int
    part::NamedTuple
end

function Doggo(leg::DoggoLeg, T_Leg)
    n_legs = length(T_leg)
    part = create_partition((3,4,))

function leg_kinematics(doggo::Doggo, q, leg::Int)
    leg = doggo.leg
    q = BlockArray(q, doggo.part)
    q_leg = q[Symbol("l" * leg)]

    # Get 2D kinematics from the leg model
    r_2d = fk(leg, q_leg)

    # Convert 2D kinematics to 3D coordinates
    ri = leg2Dto3D_2.(r)

    # Move to body frame
    Tbi = [doggo.T_leg[leg] ∘ r for r in ri]

    # Create world transformation
    r0 = q.pos
    q0 = Quaternion(q.rot)

    [[q0*Tbij.translation + r0; SVector(q0*Quaternion(Tbij.linear))] for Tbij in Tbi]
end

struct DoggoLeg{T}
    L::Vector{T}
    l::Vector{T}
    M::Vector{Diagonal{T,Vector{T}}}
    mass::Vector{T}
    gravity::T
    dof::Int
end

function DoggoLeg(L_top, L_bot, m::Vector, J::Vector)
    L1 = L_top
    L2 = L1
    L3 = L_bot
    L4 = L3
    L = [L1,L2,L3,L4]
    l = [Lk/2 for Lk in L]
    M = [Diagonal([m[k]*ones(2); J[k]]) for k = 1:4]
    DoggoLeg(L, l, M, m, 9.81, 2)
end

num_links(model::DoggoLeg) = length(model.L)

function trigtheta(doggo::DoggoLeg,d,h)
    L1,L3 = doggo.L[[1,3]]
    st = h*d/(L1*L3)
    ct = (L1^2 + L3^2 - h^2)/(2*L1*L3)
    return st,ct
end

function theta(doggo::DoggoLeg,d,h)
    atan(2*h*d, doggo.L[1]^2 + doggo.L[3]^2 - h^2)
end

function leg_dims(doggo::DoggoLeg, q)
    L1,L3 = doggo.L[[1,3]]
    a,b = q

    # Calculate internal sizes
    d = L1*sin(0.5*(a-b))
    h1 = L1*cos(0.5*(a-b))
    h2 = sqrt(L3^2 -d^2)
    h = h1 + h2
    return d, h, h1, h2
end

function fk(doggo::DoggoLeg, q)
    dof = doggo.dof
    n_links = num_links(doggo)

    # Get pieces out of model
    L1,L2,L3,L4 = doggo.L
    l1,l2,l3,l4 = doggo.l
    a,b, = q

    d,h = leg_dims(doggo, q)
    t = theta(doggo, d, h)
    st, ct = trigtheta(doggo, d, h)
    r = [zeros(eltype(q), dof) for k = 1:n_links]
    r[1] = [l1*sin(a), -l1*cos(a), a]
    r[2] = [l2*sin(b), -l2*cos(b), b]
    r[3] = [(L1 - l3*ct)*sin(a) - l3*st*cos(a),
           (-L1 + l3*ct)*cos(a) - l3*st*sin(a),
             a + t - pi]
    r[4] = [(L2 - l4*ct)*sin(b) + l4*st*cos(b),
           (-L2 + l4*ct)*cos(b) + l4*st*sin(b),
             b - t + pi]
    return r
end

function fk_joints(doggo::DoggoLeg, q)
    dof = doggo.dof
    n_links = num_links(doggo)

    # Get pieces out of model
    L1,L2,L3,L4 = doggo.L
    l1,l2,l3,l4 = doggo.l
    a,b, = q

    d,h = leg_dims(doggo, q)
    t = theta(doggo, d, h)
    st, ct = trigtheta(doggo, d, h)
    r = [zeros(eltype(q), dof) for k = 1:n_links]
    r[1] = [0,0,0]
    r[2] = [L1*sin(a), -L1*cos(a), a]
    r[3] = [L2*sin(b), -L2*cos(b), b]
    r[4] = [(L1 - L3*ct)*sin(a) - L3*st*cos(a),
           (-L1 + L3*ct)*cos(a) - L3*st*sin(a),
             a + t - pi]
    return r
end

function jacobian(doggo::DoggoLeg, q)
    dof = doggo.dof
    n_links = num_links(doggo)

    # Get pieces out of model
    L1,L2,L3,L4 = doggo.L
    l1,l2,l3,l4 = doggo.l
    a,b = q

    # Calculate internal sizes
    d = L1*sin(0.5*(a-b))
    h1 = L1*cos(0.5*(a-b))
    h2 = sqrt(L3^2 -d^2)
    h = h1 + h2

    # Calculate the sines and cosines of θ
    phi = (a-b)/2
    st,ct = trigtheta(doggo,d,h)
    dsda = ((h-d^2/h2)*cos(phi) - d*sin(phi))/(2L3)
    dcda = (h*sin(phi) + h*d/h2*cos(phi))/(2L3)
    dsdb = -dsda
    dcdb = -dcda

    # Derivatives of theta
    L = L1^2 + L3^2
    dta = (-d*(L+h^2)*sin(phi) + (-d*(L+h^2)*d/h2 + h*(L-h^2))*cos(phi))/(4*L1*L3^2)
    dtb = -dta

    # Gradient of forward Kinematics
    grad = [zeros(eltype(q),3,dof) for k = 1:n_links]
    grad[1][1,1] = l1*cos(a)
    grad[1][1,2] = 0

    grad[1][2,1] = l1*sin(a)
    grad[1][2,2] = 0

    grad[1][3,1] = 1
    grad[1][3,2] = 0

    grad[2][1,1] = 0
    grad[2][1,2] = l2*cos(b)

    grad[2][2,1] = 0
    grad[2][2,2] = l2*sin(b)

    grad[2][3,1] = 0
    grad[2][3,2] = 1

    # dx3a = L1*cos(a) + l3*cos(a2)*da2a
    # dx3b = l3*cos(a2)*da2b
    grad[3][1,1] = (L1-l3*ct)*cos(a) + l3*st*sin(a) - l3*sin(a)*dcda - l3*cos(a)*dsda
    grad[3][1,2] = -l3*sin(a)*dcdb - l3*cos(a)*dsdb

    # dy3a = L1*sin(a) + l3*sin(a2)*da2a
    # dy3b = l3*sin(a2)*da2b
    grad[3][2,1] = (L1-l3*ct)*sin(a) - l3*st*cos(a) + l3*cos(a)*dcda - l3*sin(a)*dsda
    grad[3][2,2]= l3*cos(a)*dcdb - l3*sin(a)*dsdb

    grad[3][3,1] = 1 + dta
    grad[3][3,2] = dtb

    # dx4a = l4*cos(b2)*db2a
    # dx4b = L2*cos(b) + l4*cos(b2)*db2b
    grad[4][1,1] = -l4*sin(b)*dcda + l4*cos(b)*dsda
    grad[4][1,2] = (L2 - l4*ct)*cos(b) - l4*st*sin(b) - l4*sin(b)*dcdb + l4*cos(b)*dsdb

    # dy4a = l4*sin(b2)*db2a
    # dy4b = L2*sin(b) + l4*sin(b2)*db2b
    grad[4][2,1] = l4*cos(b)*dcda + l4*sin(b)*dsda
    grad[4][2,2] = (L2 - l4*ct)*sin(b) + l4*st*cos(b) + l4*cos(b)*dcdb + l4*sin(b)*dsdb

    grad[4][3,1] = -dta
    grad[4][3,2] = 1-dtb

    return grad
end

function fk_hess(doggo::DoggoLeg, q)
    dof = doggo.dof
    n_links = num_links(doggo)

    # Get pieces out of model
    L1,L2,L3,L4 = doggo.L
    l1,l2,l3,l4 = doggo.l
    a,b, = q

    # Calculate internal sizes
    d = L1*sin(0.5*(a-b))
    h1 = L1*cos(0.5*(a-b))
    h2 = sqrt(L3^2 -d^2)
    h = h1 + h2

    # Calculate the sines and cosines of θ
    phi = (a-b)/2
    st,ct = trigtheta(doggo,d,h)
    dsda = ((h-d^2/h2)*cos(phi) - d*sin(phi))/(2L3)
    dcda = (h*sin(phi) + h*d/h2*cos(phi))/(2L3)
    dsdb = -dsda
    dcdb = -dcda

    # Hessian of sines and cosines
    dsdaa = 1/(4L3)*((-(d^2/h2^2 + 3)*h1*d/h2)*cos(phi) - (h-d^2/h2)*sin(phi) - 2h1*sin(phi) - d*cos(phi))
    dcdaa = 1/4L3*(1/h2*(-d^2 + h*h1 + (h/h2 - 1)*d^2*h1/h2 + h*h2)*cos(phi) - d*(h/h2 + (1+h1/h2))*sin(phi))
    dsdbb = dsdaa
    dcdbb = dcdaa
    dsdab = -dsdaa
    dcdab = -dcdaa

    # Derivatives of theta
    L = L1^2 + L3^2
    dta = (-d*(L+h^2)*sin(phi) + (-d*(L+h^2)*d/h2 + h*(L-h^2))*cos(phi))/(4*L1*L3^2)
    dtb = -dta

    dtaa = (( 2d^2*h*(1+h1/h2)  + (d^2 - h1*h2)*(L+h^2)/h2 - h*(L-h^2) )*sin(phi)/2  -
        ( (2h1/h2 + 1 + d^2*h1/h2^3)*(L+h^2) + (L - 3h^2 - d^2*2h/h2)*(1+h1/h2) )*d*cos(phi)/2 ) /
        (4*L1*L3^2)
    dtab = -dtaa
    dtbb = dtaa

    # Gradient of forward Kinematics
    grad = [zeros(3,dof) for k = 1:n_links]
    grad[1][1,1] = l1*cos(a)
    grad[1][1,2] = 0

    grad[1][2,1] = l1*sin(a)
    grad[1][2,2] = 0

    grad[1][3,1] = 1
    grad[1][3,2] = 0

    grad[2][1,1] = 0
    grad[2][1,2] = l2*cos(b)

    grad[2][2,1] = 0
    grad[2][2,2] = l2*sin(b)

    grad[2][3,1] = 0
    grad[2][3,2] = 1

    # dx3a = L1*cos(a) + l3*cos(a2)*da2a
    # dx3b = l3*cos(a2)*da2b
    grad[3][1,1] = (L1-l3*ct)*cos(a) + l3*st*sin(a) - l3*sin(a)*dcda - l3*cos(a)*dsda
    grad[3][1,2] = -l3*sin(a)*dcdb - l3*cos(a)*dsdb

    # dy3a = L1*sin(a) + l3*sin(a2)*da2a
    # dy3b = l3*sin(a2)*da2b
    grad[3][2,1] = (L1-l3*ct)*sin(a) - l3*st*cos(a) + l3*cos(a)*dcda - l3*sin(a)*dsda
    grad[3][2,2]= l3*cos(a)*dcdb - l3*sin(a)*dsdb

    grad[3][3,1] = 1 + dta
    grad[3][3,2] = dtb

    # dx4a = l4*cos(b2)*db2a
    # dx4b = L2*cos(b) + l4*cos(b2)*db2b
    grad[4][1,1] = -l4*sin(b)*dcda + l4*cos(b)*dsda
    grad[4][1,2] = (L2 - l4*ct)*cos(b) - l4*st*sin(b) - l4*sin(b)*dcdb + l4*cos(b)*dsdb

    # dy4a = l4*sin(b2)*db2a
    # dy4b = L2*sin(b) + l4*sin(b2)*db2b
    grad[4][2,1] = l4*cos(b)*dcda + l4*sin(b)*dsda
    grad[4][2,2] = (L2 - l4*ct)*sin(b) + l4*st*cos(b) + l4*cos(b)*dcdb + l4*sin(b)*dsdb

    grad[4][3,1] = -dta
    grad[4][3,2] = 1-dtb

    # Hessian of forward Kinematics
    hess = [zeros(3*dof, dof) for k = 1:n_links]
    hess[1][1,1] = -l1*sin(a)         # x1/αα
    hess[1][1,2] = 0                  # x1/αβ
    hess[1][2,1] = l1*cos(a)          # y1/αα
    hess[1][2,2] = 0                  # y1/αβ
    hess[1][3,1] = 0                  # θ1/αα
    hess[1][3,2] = 0                  # θ1/αβ
    hess[1][4,1] = hess[1][1,2]       # x1/βα
    hess[1][4,2] = 0                  # x1/ββ
    hess[1][5,1] = hess[1][2,2]       # y1/βα
    hess[1][5,2] = 0                  # y1/ββ
    hess[1][6,1] = hess[1][3,2]       # θ1/βα
    hess[1][6,2] = 0                  # θ1/ββ

    hess[2][1,1] = 0                  # x2/αα
    hess[2][1,2] = 0                  # x2/αβ
    hess[2][2,1] = 0                  # y2/αα
    hess[2][2,2] = 0                  # y2/αβ
    hess[2][3,1] = 0                  # θ2/αα
    hess[2][3,2] = 0                  # θ2/αβ
    hess[2][4,1] = hess[1][1,2]       # x2/βα
    hess[2][4,2] = -l2*sin(b)         # x0/ββ
    hess[2][5,1] = hess[1][2,2]       # y2/βα
    hess[2][5,2] = l2*cos(b)          # y2/ββ
    hess[2][6,1] = hess[1][3,2]       # θ2/βα
    hess[2][6,2] = 0                  # θ2/ββ

    dx3a = (L1-l3*ct)*cos(a) + l3*st*sin(a) - l3*sin(a)*dcda - l3*cos(a)*dsda
    dx3b = -l3*sin(a)*dcdb - l3*cos(a)*dsdb
    dx3aa = -(L1-l3*ct)* sin(a) - l3*cos(a)*dcda +
             l3*st * cos(a) + l3*sin(a)*dsda -
             l3*dcda*cos(a) - l3*sin(a)*dcdaa +
             l3*dsda*sin(a) - l3*cos(a)*dsdaa

    dx3aa = -(L1-l3*ct)* sin(a) - l3*cos(a)*dcda +
             l3*st * cos(a) + l3*sin(a)*dsda -
             l3*dcda*cos(a) - l3*sin(a)*dcdaa +
             l3*dsda*sin(a) - l3*cos(a)*dsdaa
    dx3ab = -l3*cos(a)*dcdb  + l3*sin(a)*dsdb -
             l3*sin(a)*dcdab - l3*cos(a)*dsdab
    dx3bb = -l3*sin(a)*dcdbb - l3*cos(a)*dsdbb

    dy3a = (L1-l3*ct)*sin(a) - l3*st*cos(a) + l3*cos(a)*dcda - l3*sin(a)*dsda
    dy3b = l3*cos(a)*dcdb - l3*sin(a)*dsdb

    dy3aa = (L1-l3*ct)*cos(a)  - l3*sin(a)*dcda +
                l3*st*sin(a)   - l3*cos(a)*dsda -
                l3*sin(a)*dcda + l3*cos(a)*dcdaa -
                l3*cos(a)*dsda - l3*sin(a)*dsdaa
    dy3ab = -l3*sin(a)*dcdb - l3*cos(a)*dsdb +
             l3*cos(a)*dcdab - l3*sin(a)*dsdab
    dy3bb = l3*cos(a)*dcdbb -l3*sin(a)*dsdbb

    dx4a = -l4*sin(b)*dcda + l4*cos(b)*dsda
    dx4b = (L2 - l4*ct)*cos(b) - l4*st*sin(b) - l4*sin(b)*dcdb + l4*cos(b)*dsdb

    dx4aa = -l4*sin(b)*dcdaa + l4*cos(b)*dsdaa
    dx4ab = -l4*dcda*cos(b) - l4*sin(b)*dcdab +
            -l4*dsda*sin(b) + l4*cos(b)*dsdab
    dx4bb = -(L2-l3*ct)*sin(b) - l4*cos(b)*dcdb +
            -l4*st*cos(b)      - l4*sin(b)*dsdb +
            -l3*dcdb*cos(b)    - l3*sin(b)*dcdbb +
            -l4*dsdb*sin(b)    + l3*cos(b)*dsdbb

    dy4a = l4*cos(b)*dcda + l4*sin(b)*dsda
    dy4b = (L2 - l4*ct)*sin(b) + l4*st*cos(b) + l4*cos(b)*dcdb + l4*sin(b)*dsdb

    dy4aa = l4*cos(b)*dcdaa + l4*sin(b)*dsdaa
    dy4ab = -l4*dcda*sin(b) + l4*cos(b)*dcdab +
             l4*dsda*cos(b) + l4*sin(b)*dsdab
    dy4bb = (L2-l4*ct)*cos(b) - l4*sin(b)*dcdb +
            -l4*st*sin(b)     + l4*cos(b)*dsdb +
            -l4*dcdb*sin(b)   + l4*cos(b)*dcdbb +
             l4*dsdb*cos(b)   + l4*sin(b)*dsdbb

    hess[3][1,1] = dx3aa
    hess[3][1,2] = dx3ab
    hess[3][2,1] = dy3aa
    hess[3][2,2] = dy3ab
    hess[3][3,1] = dtaa
    hess[3][3,2] = dtab
    hess[3][4,1] = dx3ab
    hess[3][4,2] = dx3bb
    hess[3][5,1] = dy3ab
    hess[3][5,2] = dy3bb
    hess[3][6,1] = dtab
    hess[3][6,2] = dtbb

    hess[4][1,1] = dx4aa
    hess[4][1,2] = dx4ab
    hess[4][2,1] = dy4aa
    hess[4][2,2] = dy4ab
    hess[4][3,1] = -dtaa
    hess[4][3,2] = -dtab
    hess[4][4,1] = dx4ab
    hess[4][4,2] = dx4bb
    hess[4][5,1] = dy4ab
    hess[4][5,2] = dy4bb
    hess[4][6,1] = -dtab
    hess[4][6,2] = -dtbb

    return grad, hess
end

function mass_matrix(doggo::DoggoLeg, jac::Vector{Matrix{T}}) where T
    dof = doggo.dof
    M = doggo.M
    M̄ = zeros(T,dof,dof)

    for k = 1:num_links(doggo)
        M̄ += jac[k]'*M[k]*jac[k]
    end

    return M̄
end
mass_matrix(doggo::DoggoLeg{T}, q::AbstractVector) where T = mass_matrix(doggo, jacobian(doggo, q))

get_T(doggo::DoggoLeg, q, q̇) = 0.5*q̇'mass_matrix(doggo, q)*q̇

function get_V(doggo::DoggoLeg, q)
    pos = fk(doggo, q)
    n_links = num_links(doggo)
    m = doggo.mass
    g = doggo.gravity

    heights = [pos[k][2] for k = 1:n_links]
    V = 0.0
    for k = 1:n_links
        V += m[k]*g*heights[k]
    end
    return V
end

function get_∇V(doggo::DoggoLeg, jac)
    dof = doggo.dof
    g = doggo.gravity
    ∇V = zeros(dof)
    m = doggo.mass
    for k = 1:num_links(doggo)
        ∇V += m[k]*g*jac[k][2,:]
    end
    return ∇V
end

function lagrangian(doggo::DoggoLeg, q, q̇)
    T = get_T(doggo, q, q̇)
    V = get_V(doggo, q)
    return T - V
end

function grad_L(doggo::DoggoLeg, q, q̇)
    dof = doggo.dof
    grad = zeros(2dof)

    jac, hess = fk_hess(doggo, q)
    C = comm(3,dof)
    dMdq = sum([(kron(speye(dof), jac[k]'M[k])*hess[k] +
                 kron(jac[k]'M[k], speye(dof))*C*hess[k]) for k = 1:4])
    grad[1:dof] = (0.5*kron(q̇', q̇')*dMdq)' - get_∇V(doggo, jac)
    grad[dof+1:2dof] = mass_matrix(doggo, jac)*q̇
    return grad
end

function euler_lagrange(doggo::DoggoLeg, q, qd, qdd)
    dof = doggo.dof
    jac, hess = fk_hess(doggo, q)

    dMdq = calc_dMdq(doggo, jac, hess)
    Mdot = dMdq*qd
    return reshape(Mdot,dof,dof)*qd + mass_matrix(doggo, jac)*qdd - (0.5*kron(qd',qd')*dMdq)' + get_∇V(doggo, jac)
end

function calc_dMdq(doggo::DoggoLeg, jac, hess)
    dof = doggo.dof

    C = comm(3,dof)
    M = doggo.M
    dMdq = sum([(kron(speye(dof), jac[k]'M[k])*hess[k] +
                 kron(jac[k]'M[k], speye(dof))*C*hess[k]) for k = 1:4])
    return dMdq
end

function dynamics(doggo::DoggoLeg, q, qd, tau)
    dof = doggo.dof
    jac, hess = fk_hess(doggo, q)

    dMdq = calc_dMdq(doggo, jac, hess)
    Mdot = dMdq*qd

    nl_terms = reshape(Mdot,dof,dof)*qd - (0.5*kron(qd',qd')*dMdq)' + get_∇V(doggo, jac)
    qdd = mass_matrix(doggo, jac)\(tau - nl_terms)
end

function euler_lagrange_auto(doggo, q, qd, qdd)
    dof = doggo.dof
    x = [q; qd]
    # part = create_partition2((dof,dof),(:x,:v))
    part1 = (x=1:dof,v=dof+1:2dof)
    part = create_partition2((dof,dof))
    part = NamedTuple{(:xx,:xv,:vx,:vv)}(part)
    res = DiffResults.HessianResult(x)
    lag(x) = lagrangian(doggo, x[part1.x], x[part1.v])
    ForwardDiff.hessian!(res, lag, x)
    L_hess = BlockArray(DiffResults.hessian(res), part)
    L_grad = DiffResults.gradient(res)
    # return L_hess, part
    L_hess.vx*qd + L_hess.vv*qdd - L_grad[1:2]
end

function myfun(A)
    part = create_partition2((2,2))
    # part = ntuple(i->part[i],2)
    part2 = NamedTuple{(:xx,:xv,:vx,:vv)}(part)
    A[part2.xx...] + A[part2.vv...]
end

function mypart(len)
    n = length(len)
    ntuple(i->len[i]:len[i+1],n-1)
end

function RecipesBase.plot(doggo::DoggoLeg, q, lim=3)
    pos = fk_joints(doggo, q)
    α,β = q

    p = plot(title="Doggo Leg alpha=$(round(α,digits=2)),beta=$(round(β,digits=2))",xlim=(-lim,lim),ylim=(-lim,lim),aspectratio=:equal,label="")
    plot!([pos[1][1],pos[2][1]],[pos[1][2],pos[2][2]],color=:blue,linewidth=2,label="")
    plot!([pos[1][1],pos[3][1]],[pos[1][2],pos[3][2]],color=:green,linewidth=2,label="")
    plot!([pos[2][1],pos[4][1]],[pos[2][2],pos[4][2]],color=:blue,linewidth=2,label="")
    plot!([pos[3][1],pos[4][1]],[pos[3][2],pos[4][2]],color=:green,linewidth=2,label="")
end


text:
using Atum
using Atum.Advection

using Test
using Printf
using LinearAlgebra: norm
using StaticArrays: SVector

function square(law, x⃗)
  FT = eltype(law)
  ρ = @inbounds abs(x⃗[1]) < FT(1 // 2) ? FT(2) : FT(1)
  SVector(ρ)
end

function run(A, FT, N, K; volume_form=WeakForm())
  Nq = N + 1

  law = AdvectionLaw{FT, 1}()

  cell = LobattoCell{FT, A}(Nq)
  v1d = range(FT(-1), stop=FT(1), length=K+1)
  grid = brickgrid(cell, (v1d,); periodic=(true,))

  dg = DGSEM(; law, grid, volume_form,
               surface_numericalflux = CentralFlux())

  cfl = FT(1 // 4)
  dt = cfl * min_node_distance(grid) / norm(constants(law).u⃗)
  timeend = FT(5.0)

  q = square.(Ref(law), points(grid))
  η0 = entropyintegral(dg, q)

  @info @sprintf """Starting
  N           = %d
  K           = %d
  volume_form = %s
  norm(q)     = %.16e
  η0          = %.16e
  """ N K volume_form weightednorm(dg, q) η0

  odesolver = RLSRK54(dg, q, dt)

  solve!(q, timeend, odesolver)
  ηf = entropyintegral(dg, q)

  Δη = (ηf - η0) / abs(η0)

  @info @sprintf """Finished
  norm(q) = %.16e
  ηf      = %.16e
  Δη      = %.16e
  """ weightednorm(dg, q) ηf Δη

  Δη
end

let
  A = Array
  FT = Float64
  N = 4
  K = 5

  volume_form = FluxDifferencingForm(CentralFlux())
  Δη = run(A, FT, N, K; volume_form)

  @test abs(Δη) <= 10eps(FT)
end


text:
using Extents
using Documenter

makedocs(;
    modules=[Extents],
    authors="Rafael Schouten <rafaelschouten@gmail.com>",
    repo="https://github.com/rafaqz/Extents.jl/blob/{commit}{path}#{line}",
    sitename="Extents.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rafaqz.github.io/Extents.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/rafaqz/Extents.jl",
    devbranch="main",
)


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module TwoSpin

using ExtractMacro
using ..Interface

export GraphTwoSpin

import ..Interface: energy, delta_energy, neighbors, allΔE, getN

struct GraphTwoSpin <: DiscrGraph{Int}
end

@doc """
    GraphTwoSpin() <: DiscrGraph{Int}

A trivial `DiscrGraph` type with 2 spins inteacting ferromagnetically
(\$J=1\$), without fields.

Only useful for testing/debugging purposes.
""" -> GraphTwoSpin()

getN(::GraphTwoSpin) = 2

function energy(X::GraphTwoSpin, C::Config)
    @extract C : N s
    @assert N == 2
    return -(2s[1] - 1) * (2s[2] - 1)
end

function delta_energy(X::GraphTwoSpin, C::Config, move::Int)
    @extract C : N s
    @assert 1 ≤ move ≤ N
    @assert N == 2

    return 2 * (2s[1] - 1) * (2s[2] - 1)
end

neighbors(X::GraphTwoSpin, i::Int) = return (3 - i,)
allΔE(::Type{GraphTwoSpin}) = (2,)

end


text:
include("AbstractTypes/AbstractType.jl")

mutable struct SimulationSettings <: Parameters
    day2hour::Int64 # One day has 24 hours
    year2day::Int64 # representative days in one year
    maxUnitNumbers::Int64 #Maximum units of each technology that can be once invested
    downpaymentRatio::Float64 # Initial payment
    representativeData::Dict{Symbol,Matrix} # The weight of each representative day
    investmentSequence::Symbol # sequential or parallel
    investmentModelSelection::Symbol #EMLab,EMCAS,GEP
end

mutable struct SystemSettings <: Parameters
    GenCoNumbers::Int64 # Total number of generation companies
    totalHorizon::Int64 # Simulation length
    milestoneYear::Int64
    interestRate::Float64 # risk-free interest rate in the bank
    simulationSettings::SimulationSettings
end


text:
#=______________________________________________________________________________
______________________________________________________________________________=#
function bin2hdf5(bin_path, out_path)
  try
    word_vectors = readgooglebin(bin_path)
    model_keys=collect(keys(word_vectors))
    model_values=collect(values(word_vectors))
    info("writing model ... ")
    fid=h5open(out_path,"w")
    fid["keys"]= model_keys
    dset_label = d_create(fid, "vectors", datatype(Float32),
                 dataspace(1, size(model_values[1],1), size(model_keys, 1)))
    map(i -> dset_label[:,:,i] = model_values[i] ,collect(1:size(model_keys, 1)))
    close(fid)
  catch E
    warn("Failed: $E")
    return false
  end
  return true
end

#=______________________________________________________________________________
______________________________________________________________________________=#
function readgooglebin(bin_path)
  f=open(bin_path)
  word_vectors=Dict{String,Any}()
  info("reading bin file:")
  @time full = readbytes(f)
  i=13
  ncount = 0
  while i < size(full,1)
    per_complete = i/size(full,1) * 100
    ncount += 1
    ncount % 50000 == 0 ? info("progress: $(per_complete)%") : nothing
    stringbuf = ""
    while Char(full[i])!=' ' && Char(full[i])!='\n'
      stringbuf = "$stringbuf$(string(Char(full[i])))"
      i += 1
    end
    i += 1
    floatarray=reinterpret(Float32,full[i:(1200+i-1)])
    i += 1200
    floatarray = floatarray.*(1/(norm(floatarray)))
    word_vectors[stringbuf] = floatarray
  end
  close(f)
  return word_vectors
end

#=______________________________________________________________________________
______________________________________________________________________________=#
function init(model_path)
  try
    fid = h5open(model_path, "r")
    words = read(fid["keys"])
    vectors = read(fid["vectors"])
    close(fid)
    global vectors = vectors[1,:,:]
    global words_key = Dict(zip(words,1:size(words,1)))
    return true
  catch E
    warn("Failed: $E")
    return false
  end
end

#=______________________________________________________________________________
String,String-> getsimilarity(vec1::Array,vec2::Array)
similarity between two words
getsimilarity("king","queen")
______________________________________________________________________________=#
function getsimilarity(word1::String, word2::String)
  return getsimilarity(getvector(word1, false), getvector(word2, false))
end

#=______________________________________________________________________________
Array(Float32,1),Array(Float32,1)-> Float32
similarity between two vectors
getsimilarity(getvector("king"),getvector("queen"))
______________________________________________________________________________=#
function getsimilarity(vec1, vec2)
  return vec1 != false && vec2 != false ?
    (dot(vec1,vec2)/(norm(vec1)*norm(vec2))) : false
end

#=______________________________________________________________________________
ASCIIString-> Array(Float32,1)
vector representation of word
vec1 = getvector("king")
______________________________________________________________________________=#
function getvector(word::String, default_value)
  idx = get(words_key, word, false)
  return idx != false ? vectors[:,idx] : default_value
end

#=______________________________________________________________________________
______________________________________________________________________________=#
function  getvector(words, default_value)
  vectorset = Array(Float32, size(words, 1), 300)
  for i = 1:size(words,1)
    value = getvector(words[i], false)
    if value == false
      value = default_value
    end
    vectorset[i,:] = value
  end
  return vectorset
end

#=______________________________________________________________________________
Array,Int64,Array(ASCIIString,1)-> tuple
Nearest nTop words to the given list
wordNearest(wordvectors,nTop,words)
______________________________________________________________________________=#
function wordnearest(wordvectors::Array, nTop, words)
  pq = PriorityQueue(AbstractString,Float32)
  for word in setdiff(collect(keys(words_key)), words)
    dist = getsimilarity(wordvectors, vectors[:,get(words_key, word, false)])
    enqueue!(pq, word, dist)
    (length(pq) > nTop) && dequeue!(pq)
  end
  return sort(collect(pq), by = t -> t[2], rev=true)
end

#=______________________________________________________________________________
Array(String,1),Array(String,1),Int64-> wordnearest(wv::Array,nTop,words)
Nearest nTop words to the positve and negative list
wordnearest(wordvectors,nTop,words)
______________________________________________________________________________=#
function wordnearest(positive::Array, negative::Array, nTop)
  wordvectors = zeros(300)
  words = union(positive, negative)
  for word in words
    vec = getvector(word, false)
    if vec != false
      wordvectors += word in positive ? vec : vec.*(-1)
    else
      warn("$word does not exist")
      return false
    end
  end
  return wordnearest(wordvectors, nTop, words)
end

#=______________________________________________________________________________
ASCIIString,Int64-> tuple
Nearest nTop words to the positve and negative list
wordNearest(word,nTop)
______________________________________________________________________________=#
function wordnearest(word::String, nTop)
  gc()
  collection = wordnearest([word], [], nTop)
  return collection
end

#=______________________________________________________________________________
ASCIIString-> ASCIIString
odd word out from string with words separated by space
odd_one_out("breakfast sushi dinner lunch")
______________________________________________________________________________=#
function odd1out(wordsAsString)
  words=split(wordsAsString)
  words = intersect(collect(keys(words_key)), words)
  if isempty(words)
    warn("Cannot select from empty wordList.")
    return false
  end
  oddvectors = Array(Float32,size(words,1),300)
  for i =1:size(words,1)
    oddvectors[i,:] = vectors[:,words_key[words[i]]]
  end
  meanvec = unitvec(vec(mean(oddvectors,1)),"l2")
  dists = dotmat(oddvectors, meanvec)
  return sort(collect(zip(words,dists)), by = t-> t[2])[1][1]
end

#=______________________________________________________________________________
______________________________________________________________________________=#
function unitvec(arr::Array,norm_type)
  return norm_type == "l1" ? arr./norm(arr,1) : arr./norm(arr)
end

#=______________________________________________________________________________
______________________________________________________________________________=#
function dotmat(x::Array,y::Array)
  ans = Array(Float32,size(x,1))
  ans = map(i -> (ans[i] = dot(vec(x[i,:]),vec(y))),collect(1:size(x,1)))
  return vec(ans)
end

#=______________________________________________________________________________
list,list-> Float
Compute cosine similarity between two sets of words.
n_similarity(['sushi', 'shop'], ['japanese', 'restaurant'])
______________________________________________________________________________=#
function nsimilarity(wordset1, wordset2)
  vectorset1 = getvector(wordset1, false)
  vectorset2 = getvector(wordset2, false)
  return vectorset1 != false && vectorset2 != false ?
     dot(unitvec(vec(mean(vectorset1,1)),"l2"),
     unitvec(vec(mean(vectorset2,1)),"l2")) : false
end


text:
## Sockets ##
mutable struct Socket
    data::Ptr{Cvoid}
    pollfd::_FDWatcher

    function Socket(ctx::Context, typ::Integer)
        p = ccall((:zmq_socket, libzmq), Ptr{Cvoid}, (Ptr{Cvoid}, Cint), ctx, typ)
        if p == C_NULL
            throw(StateError(jl_zmq_error_str()))
        end
        socket = new(p)
        setfield!(socket, :pollfd, _FDWatcher(fd(socket), #=readable=#true, #=writable=#false))
        finalizer(close, socket)
        push!(getfield(ctx, :sockets), WeakRef(socket))
        return socket
    end
    Socket(typ::Integer) = Socket(context(), typ)
end

function Socket(f::Function, args...)
    socket = Socket(args...)
    try
        f(socket)
    finally
        close(socket)
    end
end

Base.unsafe_convert(::Type{Ptr{Cvoid}}, s::Socket) = getfield(s, :data)

Base.isopen(socket::Socket) = getfield(socket, :data) != C_NULL
function Base.close(socket::Socket)
    if isopen(socket)
        close(getfield(socket, :pollfd), #=readable=#true, #=writable=#false)
        rc = ccall((:zmq_close, libzmq), Cint,  (Ptr{Cvoid},), socket)
        setfield!(socket, :data, C_NULL)
        if rc != 0
            throw(StateError(jl_zmq_error_str()))
        end
    end
end

# Raw FD access
if Sys.isunix()
    Base.fd(socket::Socket) = RawFD(socket.fd)
end
if Sys.iswindows()
    using Base.Libc: WindowsRawSocket
    Base.fd(socket::Socket) = WindowsRawSocket(convert(Ptr{Cvoid}, socket.fd))
end

Base.wait(socket::Socket) = wait(getfield(socket, :pollfd), readable=true, writable=false)
Base.notify(socket::Socket) = @preserve socket uv_pollcb(getfield(socket, :pollfd).handle, Int32(0), Int32(UV_READABLE))

function Sockets.bind(socket::Socket, endpoint::AbstractString)
    rc = ccall((:zmq_bind, libzmq), Cint, (Ptr{Cvoid}, Ptr{UInt8}), socket, endpoint)
    if rc != 0
        throw(StateError(jl_zmq_error_str()))
    end
end

function Sockets.connect(socket::Socket, endpoint::AbstractString)
    rc=ccall((:zmq_connect, libzmq), Cint, (Ptr{Cvoid}, Ptr{UInt8}), socket, endpoint)
    if rc != 0
        throw(StateError(jl_zmq_error_str()))
    end
end


text:
################################################################################
#
#   BAND STRUCTURE PLOTTING
#
################################################################################

# function to plot a bandstructure
function plotBandstructure(
            bs :: BS
            ;
            new_figure :: Bool = true,
            figsize :: Tuple = (6,4),
            color :: Vector{<:Integer} = [100,120,255],
            kwargs...
        ) where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, BS <: AbstractBandstructure{P,H}}

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # configure plot environment
    #rc("font", family="serif")

    # create a new figure
    if new_figure
        fig = figure(figsize=figsize)
    else
        fig = gcf()
    end

    ###########################
    #   PLOT BANDS
    ###########################

    # collect the segment breaks
    k_point_indices = zeros(Int64, length(energies(bs))+1)
    # get the labels of these breaks
    k_point_labels = label.(path(bs))

    # plot the band structure
    for s in 1:length(energies(bs))
        # push the next k point as index into the array
        k_point_indices[s+1] = length(energies(bs)[s][1]) + k_point_indices[s] - 1
        # plotting segment s, collecting x values for bands
        xvals = range(k_point_indices[s], stop=k_point_indices[s+1], length=length(energies(bs)[s][1]))
        # plot all bands
        for band in energies(bs)[s]
            plot(xvals,band, color=color./255)
        end
    end

    ###########################
    #   SET ALL TICKS (POINTS)
    ###########################

    # get the current axis
    ax = gca()
    axx = ax.get_xaxis()

    # configure tick labels on x axis
    xticks(k_point_indices, k_point_labels)
    # momentum limits (x axis)
    xlim(0, k_point_indices[end]+1)

    # configure ticks on x axis
    axx.set_tick_params(which="both", direction="out")
    axx.set_tick_params(which="top", color="none")

    # configure ticks on x axis
    axy = ax.get_yaxis()
    axy.set_tick_params(which="both", direction="out")

    # plot vertical lines for each point
    for p in k_point_indices
        axvline(p,color=[0.6, 0.6, 0.6], linestyle="--")
    end

    ###########################
    #   CONFIGURE AXIS
    ###########################

    # label the axis
    xlabel("momentum")
    ylabel("energy")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # return the figure object
    return fig
end

# pass unknown arguments directly to construction of band structure
function plotBandstructure(
        args...
        ;
        kwargs...
    )

    # create and plot a bandstructure
    plotBandstructure(getBandstructure(args...); kwargs...)
end

# export the function
export plotBandstructure


text:
using Mill, Test, Flux
using Mill: partialeval, identity_model

@testset "partialeval" begin
    metadata = fill("metadata", 4)
    an1 = ArrayNode(rand(3,4))
    b = BagNode(an1, [1:4, 0:-1], metadata)
    an2 = ArrayNode(randn(5, 4))
    wb = WeightedBagNode(an2, [1:2,3:4], rand(1:4, 4), metadata)
    pn = ProductNode((b=b,wb=wb))
    an3 = ArrayNode(rand(10, 2))
    ds = ProductNode((pn, an3))
    # printtree(ds)
    m = reflectinmodel(ds, d -> Chain(Dense(d, 4, relu), Dense(4,3)), d-> SegmentedMeanMax(d))

    @test partialeval(m.ms[2], ds.data[2], an3)[1] === m.ms[2]
    @test partialeval(m.ms[2], ds.data[2], an3)[2] === ds.data[2]
    @test partialeval(m.ms[2], ds.data[2], an1)[1] === identity_model()
    @test partialeval(m.ms[2], ds.data[2], an1)[2].data ≈ m.ms[2](ds.data[2]).data
    tm, td = partialeval(m.ms[2], ds.data[2], an3)
    @test tm(td).data ≈ m.ms[2](ds.data[2]).data
    @test partialeval(m.ms[2], ds.data[2], an1)[2].data ≈ m.ms[2](ds.data[2]).data

    @test partialeval(m.ms[1][:b], ds.data[1][:b], an1)[1] === m.ms[1][:b]
    @test partialeval(m.ms[1][:b], ds.data[1][:b], an1)[2] === ds.data[1][:b]
    @test partialeval(m.ms[1][:b], ds.data[1][:b], b)[1] === m.ms[1][:b]
    @test partialeval(m.ms[1][:b], ds.data[1][:b], b)[2] === ds.data[1][:b]
    tm, td = partialeval(m.ms[1][:b], ds.data[1][:b], an3)
    @test tm(td).data ≈ m.ms[1][:b](ds.data[1][:b]).data
    @test partialeval(m.ms[1][:b], ds.data[1][:b], an3)[1] === identity_model()
    @test partialeval(m.ms[1][:b], ds.data[1][:b], an3)[2].data ≈ m.ms[1][:b](ds.data[1][:b]).data

    @test partialeval(m.ms[1][:wb], ds.data[1][:wb], an2)[1] === m.ms[1][:wb]
    @test partialeval(m.ms[1][:wb], ds.data[1][:wb], an2)[2] === ds.data[1][:wb]
    @test partialeval(m.ms[1][:wb], ds.data[1][:wb], wb)[1] === m.ms[1][:wb]
    @test partialeval(m.ms[1][:wb], ds.data[1][:wb], wb)[2] === ds.data[1][:wb]
    tm, td = partialeval(m.ms[1][:wb], ds.data[1][:wb], an3)
    @test tm(td).data ≈ m.ms[1][:wb](ds.data[1][:wb]).data
    @test partialeval(m.ms[1][:wb], ds.data[1][:wb], an3)[1] === identity_model()
    @test partialeval(m.ms[1][:wb], ds.data[1][:wb], an3)[2].data ≈ m.ms[1][:wb](ds.data[1][:wb]).data

    @test partialeval(m.ms[1], ds.data[1], an1)[1].ms[1] === m.ms[1].ms[1]
    @test partialeval(m.ms[1], ds.data[1], an1)[2].data[1] === ds.data[1].data[1]
    @test partialeval(m.ms[1], ds.data[1], an1)[1].ms[2] === identity_model()
    @test partialeval(m.ms[1], ds.data[1], an1)[2].data[2].data ≈ m.ms[1].ms[2](ds.data[1].data[2]).data
    @test partialeval(m.ms[1], ds.data[1], an3)[1] === identity_model()
    @test partialeval(m.ms[1], ds.data[1], an3)[2].data ≈ m.ms[1](ds.data[1]).data
    tm, td = partialeval(m.ms[1], ds.data[1], an1)
    @test tm(td).data ≈ m.ms[1](ds.data[1]).data
    @test partialeval(m.ms[2], ds.data[2], an1)[2].data ≈ m.ms[2](ds.data[2]).data

    @test partialeval(m, ds, 1)[1] === identity_model()
    @test partialeval(m, ds, 1)[2].data ≈ m(ds).data
    tm, td = partialeval(m, ds, an1)
    @test tm(td).data ≈ m(ds).data
end


text:
#=

# Vulkan functions

Functions in C behave differently that in Julia. In particular, they can't return multiple values and mutate pointer memory instead. Other patterns emerge from the use of pointers with a separately-provided length, where a length/size parameter can be queried, so that you build a pointer with the right size, and pass it in to the API to be filled with data.
All these patterns were automated, so that wrapper functions feel a lot more natural and straightforward for Julia users than the API functions.

## Implicit return values

Functions almost never directly return a value in Vulkan, and usually return either a return code or nothing. This is a limitation of C where only one value can be returned. Instead, they fill pointers with data, and it is your responsibility to initialize them before the call and dereference them afterwards. Here is an example:

=#

using Vulkan
using Vulkan.core

function example_create_instance()
    instance_ref = Ref{VkInstance}()
    ## we'll cheat a bit for the create info
    code = vkCreateInstance(
        InstanceCreateInfo([], []), # create info
        C_NULL,                     # allocator
        instance_ref,
    )

    @assert code == VK_SUCCESS

    instance_ref[]
end

example_create_instance()

#=

We did not create a `VkInstanceCreateInfo` to stay concise. Note that the create info structure can be used as is by the `vkCreateInstance`, even if it is a wrapper. Indeed, it implements `Base.cconvert` and `Base.unsafe_convert` to automatically interface with the C API.

All this setup code is now automated by the wrapper:

=#

instance = unwrap(create_instance(InstanceCreateInfo([], []); allocator=C_NULL))

#=

Note that the `create_instance` wrapper returns a `ResultTypes.Result` type, so the result must be unwrapped to be used. To know more about `unwrap` and the `Result` type see the [Error handling](@ref Error-handling) section.

When there are multiple implicit return values (i.e. multiple pointers being written to), they are returned as a tuple:

```julia
actual_data_size, data = unwrap(get_pipeline_cache_data(device, pipeline_cache, data_size))
```

## Queries

### Enumerated items

Sometimes, when enumerating objects or properties for example, a function may need to be called twice: a first time for returning the number of elements to be enumerated, then a second time with an initialized array of the right length to be filled with Vulkan objects:

=#

function example_enumerate_physical_devices(instance)
    pPhysicalDeviceCount = Ref{UInt32}(0)

    ## get the length in pPhysicalDeviceCount
    code = vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, C_NULL)

    @assert code == VK_SUCCESS
    ## initialize the array with the returned length
    pPhysicalDevices = Vector{VkPhysicalDevice}(undef, pPhysicalDeviceCount[])

    ## fill the array
    code = vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, pPhysicalDevices)
    @assert code == VK_SUCCESS

    pPhysicalDevices
end

example_enumerate_physical_devices(instance)

# The relevant enumeration functions are wrapped with this, so that only one call needs to be made, without worrying about creating intermediate arrays:

unwrap(enumerate_physical_devices(instance))

#=

### Incomplete retrieval

Some API commands such as `vkEnumerateInstanceLayerProperties` may return a `VK_INCOMPLETE` code indicating that some items could not be written to the provided array. This happens if the number of available items changes after that the length is obtained, making the array too small. In this case, it is recommended to simply query the length again, and provide a vector of the updated size, starting over if the number of items changes again. To avoid doing this by hand, this step is automated in a while loop. Here is what it may look like:

=#

function example_enumerate_physical_devices_2(instance)
    pPhysicalDeviceCount = Ref{UInt32}(0)

    @assert vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, C_NULL) == VK_SUCCESS
    pPhysicalDevices = Vector{VkPhysicalDevice}(undef, pPhysicalDeviceCount[])
    code = vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, pPhysicalDevices)

    while code == VK_INCOMPLETE
        @assert vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, C_NULL) == VK_SUCCESS
        pPhysicalDevices = Vector{VkPhysicalDevice}(undef, pPhysicalDeviceCount[])
        code = vkEnumeratePhysicalDevices(instance, pPhysicalDeviceCount, pPhysicalDevices)
    end

    pPhysicalDevices
end

example_enumerate_physical_devices_2(instance)

# Fortunately, this is automated for you and you can just call

unwrap(enumerate_physical_devices(instance))


text:
println(sqrt(1im))
println(6//9)
println(numerator(6//9))
println(denominator(6//9))

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
Pkg.develop(Pkg.PackageSpec(path=joinpath(@__DIR__, "ReferenceTests")))

using ReferenceTests

ReferenceTests.record_tests()
# needs GITHUB_TOKEN to be set:
# ReferenceTests.upload_reference_images()
# Needs a backend to actually have something recoreded:
# ReferenceTests.reference_tests(recorded)


text:
using Documenter, DeterminantalPointProcesses

makedocs(;
    modules=[DeterminantalPointProcesses],
    format=Documenter.Writers.HTMLWriter.HTML(;
        assets=["assets/icon.ico"], analytics="UA-129106538-2"
    ),
    sitename="DeterminantalPointProcesses",
    authors="Theo Galy-Fajou, Maruan Al-Shedivat",
    pages=["Home" => "index.md"],
)

deploydocs(;
    deps=Deps.pip("mkdocs", "python-markdown-math"),
    repo="github.com/theogf/DeterminantalPointProcesses.jl.git",
    target="build",
)


text:
using MyFiles
using Test

@testset "MyFiles.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
@doc Markdown.doc"""
    rref(M::SMat{T}; truncate = false) where {T <: FieldElement} -> (Int, SMat{T})

Return a tuple $(r, A)$ consisting of the rank $r$ of $M$ and a reduced row echelon
form $A$ of $M$.
If the function is called with `truncate = true`, the result will not contain zero
rows, so `nrows(A) == rank(M)`.
"""
function rref(A::SMat{T}; truncate::Bool = false) where {T <: FieldElement}
  B = deepcopy(A)
  r = rref!(B, truncate = truncate)
  return r, B
end

# This does not really work in place, but it certainly changes A
function rref!(A::SMat{T}; truncate::Bool = false) where {T <: FieldElement}
  B = sparse_matrix(base_ring(A))
  B.c = A.c
  number_of_rows = A.r

  # Remove empty rows, so they don't get into the way when we sort
  i = 1
  while i <= length(A.rows)
    if iszero(A.rows[i])
      deleteat!(A.rows, i)
    else
      i += 1
    end
  end

  # Prefer sparse rows and, if the number of non-zero entries is equal, rows
  # with more zeros in front. (Appears to be a good heuristic in practice.)
  rows = sort!(A.rows, lt = (x, y) -> length(x) < length(y) || (length(x) == length(y) && x.pos[1] > y.pos[1]))

  for r in rows
    b = _add_row_to_rref!(B, r)
    if nrows(B) == ncols(B)
      break
    end
  end

  A.nnz = B.nnz
  A.rows = B.rows
  rankA = B.r
  if !truncate
    while length(A.rows) < number_of_rows
      push!(A.rows, sparse_row(base_ring(A)))
    end
  else
    A.r = B.r
  end
  return rankA
end

# Reduce v by M and if the result is not zero add it as a row (and then reduce
# M to maintain the rref).
# Return true iff v is not in the span of the rows of M.
# M is supposed to be in rref and both M and v are changed in place.
function _add_row_to_rref!(M::SMat{T}, v::SRow{T}) where { T <: FieldElem }
  if iszero(v)
    return false
  end

  pivot_found = false
  s = one(base_ring(M))
  i = 1
  new_row = 1
  while i <= length(v)
    c = v.pos[i]
    r = find_row_starting_with(M, c)
    if r > nrows(M) || M.rows[r].pos[1] > c
      # We found an entry in a column of v, where no other row of M has the pivot.
      @assert !iszero(v.values[i])
      i += 1
      if pivot_found
        # We already found a pivot
        continue
      end

      @assert i == 2 # after incrementing
      pivot_found = true
      new_row = r
      continue
    end

    # Reduce the entries of v by M.rows[r]
    t = -v.values[i] # we assume M.rows[r].pos[1] == 1 (it is the pivot)
    v = add_scaled_row!(M.rows[r], v, t)
    # Don't increase i, we deleted the entry
  end
  if !pivot_found
    return false
  end

  # Multiply v by inv(v.values[1])
  if !isone(v.values[1])
    t = inv(v.values[1])
    for j = 2:length(v)
      v.values[j] = mul!(v.values[j], v.values[j], t)
    end
    v.values[1] = one(base_ring(M))
  end
  insert!(M, new_row, v)

  # Reduce the rows above the newly inserted one
  for i = 1:new_row - 1
    c = M.rows[new_row].pos[1]
    j = searchsortedfirst(M.rows[i].pos, c)
    if j > length(M.rows[i].pos) || M.rows[i].pos[j] != c
      continue
    end

    t = -M.rows[i].values[j]
    l = length(M.rows[i])
    M.rows[i] = add_scaled_row!(M.rows[new_row], M.rows[i], t)
    if length(M.rows[i]) != l
      M.nnz += length(M.rows[i]) - l
    end
  end
  return true
end

###############################################################################
#
#   Kernel
#
###############################################################################

@doc Markdown.doc"""
    nullspace(M::SMat{T}) where {T <: FieldElement}

Return a tuple $(\nu, N)$ consisting of the nullity $\nu$ of $M$ and
a basis $N$ (consisting of column vectors) for the right nullspace of $M$,
i.e. such that $MN$ is the zero matrix. If $M$ is an $m\times n$ matrix
$N$ will be a $n\times \nu$ matrix in dense representation. The columns of $N$
are in upper-right reduced echelon form.
"""
function nullspace(M::SMat{T}) where {T <: FieldElement}
  rank, A = rref(M, truncate = true)
  nullity = ncols(M) - rank
  K = base_ring(M)
  X = zero_matrix(K, ncols(M), nullity)
  if rank == 0
    for i = 1:nullity
      X[i, i] = one(K)
    end
  elseif nullity != 0
    r = 1
    k = 1
    for c = 1:ncols(A)
      if r <= rank && A.rows[r].pos[1] == c
        r += 1
        continue
      end

      for i = 1:r - 1
        j = searchsortedfirst(A.rows[i].pos, c)
        if j > length(A.rows[i].pos) || A.rows[i].pos[j] != c
          continue
        end
        X[A.rows[i].pos[1], k] = -A.rows[i].values[j]
      end
      X[c, k] = one(K)
      k += 1
    end
  end
  return nullity, X
end

@doc Markdown.doc"""
    left_kernel(M::SMat{T}) where {T <: FieldElement}

Return a tuple $\nu, N$ where $N$ is a matrix whose rows generate the
left kernel of $M$, i.e. $NM = 0$ and $\nu$ is the rank of the kernel.
If $M$ is an $m\times n$ matrix $N$ will be a $\nu\times m$ matrix in dense
representation. The rows of $N$ are in lower-left reduced echelon form.
"""
function left_kernel(M::SMat{T}) where T <: FieldElement
  n, N = nullspace(transpose(M))
  return n, transpose(N)
end

@doc Markdown.doc"""
    right_kernel(M::SMat{T}) where {T <: FieldElement}

Return a tuple $\nu, N$ where $N$ is a matrix whose columns generate the
right kernel of $M$, i.e. $MN = 0$ and $\nu$ is the rank of the kernel.
If $M$ is an $m\times n$ matrix $N$ will be a $n \times \nu$ matrix in dense
representation. The columns of $N$ are in upper-right reduced echelon form.
"""
function right_kernel(M::SMat{T}) where T <: FieldElement
  return nullspace(M)
end

@doc Markdown.doc"""
    kernel(M::SMat{T}; side::Symbol = :right) where {T <: FieldElement}

Return a tuple $(n, N)$, where n is the rank of the kernel and $N$ is a
basis for it. If side is $:right$ or not specified, the right kernel is
computed, i.e. the matrix of columns whose span gives the right kernel
space. If side is $:left$, the left kernel is computed, i.e. the matrix
of rows whose span is the left kernel space.
"""
function kernel(M::SMat{T}; side::Symbol = :right) where T <: FieldElement
  if side == :right
    return right_kernel(M)
  elseif side == :left
    return left_kernel(M)
  else
    error("Unsupported argument: :$side for side: must be :left or :right")
  end
end