text:
__precompile__(true)
module IUP

export
    Ihandle,
    IupAlarm,
    IupBackgroundBox,
    IupButton,
    IupCanvas, IupCbox, IupClose, IupControlsOpen,
    IupDestroy, IupDialog,
    IupFill, IupFrame,
    IupGetAttribute, IupGetClassName, IupGetDialog, IupGetFloat, IupGetFile, IupGetGlobal, IupGetHandle, IupGetInt,
    IupGetDouble, IupGetDoubleId, IupGetInt2, IupGridBox, IupGetParam, IupGetParamv, IupMessagef,
    IupHbox,
    IupHide,
    IupImage, IupImageRGBA, IupItem,
    IupLabel, IupList,
    IupLoad,
    IupMainLoop, IupMap, IupMatrix,
    IupMenu,
    IupMessage,
    IupMultiLine,
    IupOpen,
    IupPopup, IupProgressBar,
    IupRadio,
    IupRefresh,
    IupSeparator,
    IupSetAttribute, IupSetAttributes, IupSetAttributeHandle, IupSetCallback, IupSetCallbacks, IupSetFunction,
    IupSetFloat, IupSetDouble, IupSetDoubleId, IupSetGlobal, IupSetHandle, IupSetLanguage,
    IupShow,
    IupShowXY,
    IupStoreAttribute,
    IupSubmenu,
    IupTabs, IupTabsv, IupText, IupToggle, IupTree,
    IupVal, IupVbox,
    IupDial, IupGauge, IupColorBrowser, IupColorbar,
    #
    # From iup_plot
    #
    IupPlotOpen,
    IupPlot,
    IupPlotBegin,
    IupPlotAdd,
    IupPlotAddStr,
    IupPlotEnd,
    IupPlotLoadData,
    IupPlotInsertStr,
    IupPlotInsert,
    IupPlotInsertStrPoints,
    IupPlotInsertPoints,
    IupPlotAddPoints,
    IupPlotAddStrPoints,
    IupPlotGetSample,
    IupPlotGetSampleStr,
    IupPlotGetSampleSelection,
    IupPlotSetSample,
    IupPlotSetSampleStr,
    IupPlotSetSampleSelection,
    IupPlotTransform,
    IupPlotTransformTo,
    IupPlotPaintTo,
    #
    # Need to add many more
    #
    OBJC_NEW_PROPERTIES, IUP_CLOSE, IUP_CONTINUE,
    IUP_ERROR, IUP_NOERROR, IUP_OPENED, IUP_INVALID, IUP_INVALID_ID, IUP_IGNORE,
    IUP_DEFAULT, IUP_CENTER, IUP_LEFT, IUP_TOP, IUP_BOTTOM, IUP_RIGHT, IUP_MOUSEPOS, IUP_CURRENT, IUP_CENTERPARENT,
    IUP_BUTTON1, IUP_BUTTON2, IUP_BUTTON3, IUP_BUTTON4, IUP_BUTTON5,
    IUP_IMAGE, IUP_TITLE, IUP_VALUE, IUP_ACTIVE, IUP_YES, IUP_NO, IUP_X, IUP_Y, IUP_VISIBLE

include("libiup.jl")
include("iupcontrols.jl")
include("iup_plot.jl")
include("iup_mglplot.jl")

include("IUP_IM.jl")
include("IUP_CD.jl")

end  # module


text:
@testset "Multivariate" begin

    @testset "Simple test" begin

        data = [1.0 1.0   2.0
                0.0 0.0   4.0
                7.0 1.5  -4.5
                2.0 2.0 -17.0 # outlier
                0.0 8.6  -4.6]

        # This is the most complete way of defining the arguments for model and gmodel!
        model = (x::Union{Vector{Float64}, SubArray}, θ::Vector{Float64}) -> θ[1] *
            x[1] + θ[2] * x[2] + θ[3]

        gmodel! = (g::SubArray, x::Union{Vector{Float64}, SubArray},
                   θ::Vector{Float64}) -> begin

            g[1] = x[1]
            g[2] = x[2]
            g[3] = 1.0

        end

        n = 3

        p = 4

        θsol = [- 1.0, - 1.0, 4.0]

        rout = lmlovo(model, gmodel!, data, n, p; ε=1.0e-8)

        @test rout.solution ≈ θsol atol=1.0e-4
        @test rout.outliers == [4]

        rout = raff(model, gmodel!, data, n)

        @test rout.solution ≈ θsol atol=1.0e-3
        @test rout.p == 4
        @test rout.outliers == [4]

        # Using automatic differentiation

        # This is the most complete way of defining the arguments for
        # model when automatic differentiation is being used
        model = (x::Union{Vector{Float64}, SubArray}, θ) -> θ[1] *
            x[1] + θ[2] * x[2] + θ[3]

        rout = lmlovo(model, data, n, p; ε=1.0e-8)

        @test rout.solution ≈ θsol atol=1.0e-4
        @test rout.outliers == [4]

        rout = raff(model, data, n)

        @test rout.solution ≈ θsol atol=1.0e-3
        @test rout.p == 4
        @test rout.outliers == [4]

    end

    @testset "Circle" begin

        data = [2.0              	1.0             	0.0
                1.766044443118978	1.6427876096865393	-3.3306690738754696e-16
                1.1736481776669305	1.9848077530122081	2.220446049250313e-16
                0.5000000000000002	1.8660254037844388	0.0
                0.06030737921409168	1.3420201433256689	-1.1102230246251565e-16
                0.06030737921409157	0.6579798566743313	0.0
                0.49999999999999956	0.13397459621556151	0.0
                3.17364817766693	0.015192246987791869	0.0 # noise
                1.766044443118978	0.3572123903134604	0.0]

        θsol = [1.0, 1.0, 1.0]

        model = (x::Union{Vector{Float64}, SubArray}, θ) -> (x[1] - θ[1])^2 + (x[2] - θ[2])^2 - θ[3]

        n = 3

        p = 8

        rout = lmlovo(model, [0.0, 0.0, 2.0], data, n, p; ε=1.0e-8)

        @test rout.solution ≈ θsol atol=1.0e-4
        @test rout.outliers == [8]

        rout = raff(model, data, n; ε=1.0e-8)

        @test rout.p == 8
        @test rout.outliers == [8]
        @test rout.solution ≈ θsol atol=1.0e-4

    end

end

ext: jl


text:
function _test_num_constraints(bridged_mock, F, S, n)
    @test MOI.get(bridged_mock, MOI.NumberOfConstraints{F,S}()) == n
    @test length(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{F,S}())) == n
    ret = (F, S) in MOI.get(bridged_mock, MOI.ListOfConstraintTypesPresent())
    @test ret == !iszero(n)
    return
end

function _warn_incomplete_list_num_constraints(BT, list_num_constraints)
    for (S,) in MOI.Bridges.added_constrained_variable_types(BT)
        F = MOI.Utilities.variable_function_type(S)
        if !any(c -> c[1] == F && c[2] == S, list_num_constraints)
            error(
                "Bridges of type $BT add constrained variable in $S but " *
                "their number is not tested.",
            )
        end
    end
    for (F, S) in MOI.Bridges.added_constraint_types(BT)
        if !any(c -> c[1] == F && c[2] == S, list_num_constraints)
            error(
                "Bridges of type $BT add $F-in-$S constraints but their " *
                "number is not tested.",
            )
        end
    end
    return
end

"""
    _test_delete_bridge(
        m::MOI.Bridges.AbstractBridgeOptimizer,
        ci::MOI.ConstraintIndex{F, S},
        nvars::Int,
        list_num_constraints::Tuple;
        used_bridges = 1,
        num_bridged = 1,
    ) where {F, S}

Test deletion of the constraint `ci` in model `m`.

* The number of variables added in `m` is `nvars` (ignoring the variables added
  by the bridges).
* For each `(Fi, Si, n)` in `list_num_constraints`, the number of `Fi`-in-`Si`
  constraints added to `m` is `n` (ignoring the constraints added by the
  bridges). All constraint types added by the bridge bridging `ci` should appear
  in `list_num_constraints`. Otherwise a error is thrown to ensure good test
  coverage.
* The number of bridges that will be deleted when deleted `ci` is
  `used_bridges`.
* The number of `F`-in-`S` constraints is `num_bridged`.
"""
function _test_delete_bridge(
    m::MOI.Bridges.AbstractBridgeOptimizer,
    ci::MOI.ConstraintIndex{F,S},
    nvars::Int,
    list_num_constraints::Tuple;
    used_bridges = 1,
    num_bridged = 1,
) where {F,S}
    _warn_incomplete_list_num_constraints(
        typeof(MOI.Bridges.bridge(m, ci)),
        list_num_constraints,
    )
    function num_bridges()
        return count(bridge -> true, values(MOI.Bridges.Constraint.bridges(m)))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    _test_num_constraints(m, F, S, num_bridged)
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    @test MOI.is_valid(m, ci)
    MOI.delete(m, ci)
    @test_throws MOI.InvalidIndex(ci) MOI.delete(m, ci)
    @test !MOI.is_valid(m, ci)
    @test num_bridges() == start_num_bridges - used_bridges
    _test_num_constraints(m, F, S, num_bridged - 1)
    # As the bridge has been removed, if the constraints it has created were not
    # removed, it wouldn't be there to decrease this counter anymore.
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    return
end

# Test deletion of variable bridge used for variable `vi`
function _test_delete_bridged_variable(
    m::MOI.Bridges.AbstractBridgeOptimizer,
    vi::MOI.VariableIndex,
    S::Type,
    nvars::Int,
    list_num_constraints::Tuple;
    used_bridges = 1,
    num_bridged = 1,
    used_constraints = 1,
)
    _warn_incomplete_list_num_constraints(
        typeof(MOI.Bridges.bridge(m, vi)),
        list_num_constraints,
    )
    function num_bridges()
        return count(bridge -> true, values(MOI.Bridges.Variable.bridges(m)))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    if S != MOI.Reals
        F =
            S <: MOI.AbstractScalarSet ? MOI.VariableIndex :
            MOI.VectorOfVariables
        _test_num_constraints(m, F, S, num_bridged)
    end
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    @test MOI.is_valid(m, vi)
    MOI.delete(m, vi)
    @test_throws MOI.InvalidIndex(vi) MOI.delete(m, vi)
    @test !MOI.is_valid(m, vi)
    @test num_bridges() == start_num_bridges - used_bridges
    if S != MOI.Reals
        _test_num_constraints(m, F, S, num_bridged - used_constraints)
    end
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars - 1
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars - 1
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    return
end

# Test deletion of variable bridge used for vector of variables `vis`
function _test_delete_bridged_variables(
    m::MOI.Bridges.AbstractBridgeOptimizer,
    vis::Vector{MOI.VariableIndex},
    S::Type,
    nvars::Int,
    list_num_constraints::Tuple;
    used_bridges = 1,
    num_bridged = 1,
)
    _warn_incomplete_list_num_constraints(
        typeof(MOI.Bridges.bridge(m, vis[1])),
        list_num_constraints,
    )
    function num_bridges()
        return count(bridge -> true, values(MOI.Bridges.Variable.bridges(m)))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    if S != MOI.Reals
        F =
            S <: MOI.AbstractScalarSet ? MOI.VariableIndex :
            MOI.VectorOfVariables
        _test_num_constraints(m, F, S, num_bridged)
    end
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    @test all(vi -> MOI.is_valid(m, vi), vis)
    MOI.delete(m, vis)
    @test_throws MOI.InvalidIndex(vis[1]) MOI.delete(m, vis)
    @test all(vi -> !MOI.is_valid(m, vi), vis)
    @test num_bridges() == start_num_bridges - used_bridges
    if S != MOI.Reals
        _test_num_constraints(m, F, S, num_bridged - 1)
    end
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars - length(vis)
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars - length(vis)
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    return
end

function _test_delete_objective(
    m::MOI.Bridges.AbstractBridgeOptimizer,
    nvars::Int,
    list_num_constraints::Tuple;
    used_bridges = 1,
)
    _warn_incomplete_list_num_constraints(
        typeof(
            MOI.Bridges.Objective.root_bridge(MOI.Bridges.Objective.bridges(m)),
        ),
        list_num_constraints,
    )
    function num_bridges()
        return length(MOI.Bridges.Objective.bridges(m))
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    MOI.set(m, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(m, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    @test num_bridges() == start_num_bridges - used_bridges
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    @test length(MOI.get(m, MOI.ListOfVariableIndices())) == nvars
    for num_constraints in list_num_constraints
        _test_num_constraints(m, num_constraints...)
    end
    return
end


text:
using ProgressBars, Plots, StatsBase, LaTeXStrings, JLD, Distributions, StatsPlots

Plots.default(titlefontsize=12, tickfontsize=10, labelfontsize=16,
    fontfamily="Computer Modern")

SampleNumber = 100000
NList = [5, 10, 100, 1000]
Data = [[sum(rand(N)) for i in ProgressBar(1:SampleNumber)] for N ∈ NList]

cd(dirname(@__FILE__))
save("../../Data/Q1/Q1-samples.jld","Samples", Data)

plots = []
colors = [:steelblue, :purple, :green, :brown]
for i ∈ 1:4
    Plt = begin
    histogram(Data[i],nbins = 60, normalize = true, c = colors[i], label = L"Histogram\ of\ Samples")
    plot!(Normal(mean(Data[i]),std(Data[i])), c = :black, lineWidth = 10,
        label = L"\mathcal{Normal}(\mu = %$(round(mean(Data[i]),digits = 2)), \sigma = %$(round(std(Data[i]),digits = 2)))")
    plot!(frame = :box, title = "N = $(NList[i]), $SampleNumber Samples",
        xlabel = L"X_N", ylabel = L"Counting\ ratio")
    end
    push!(plots,Plt)
end

plot(plots..., size = (1000,1000), plot_title = "Normalized Distribution of Samples")
savefig("../../Figs/Q1/Q1-Hist.pdf")


text:
using Markdown
using Test
using LinearAlgebra
using TestOptinum
using Optinum

include("../src/Algorithme_De_Newton.jl")
include("../src/Gradient_Conjugue_Tronque.jl")
include("../src/Lagrangien_Augmente.jl")
include("../src/Pas_De_Cauchy.jl")
include("../src/Regions_De_Confiance.jl")

TestOptinum.cacher_stacktrace()

affiche = false
println("affiche = ",affiche)
# Tester l'ensemble des algorithmes
@testset "Test SujetOptinum" begin
	# Tester l'algorithme de Newton
	tester_algo_newton(affiche,Algorithme_De_Newton)

	# Tester l'algorithme du pas de Cauchy
	tester_pas_de_cauchy(affiche,Pas_De_Cauchy)

	# Tester l'algorithme du gradient conjugué tronqué
	tester_gct(affiche,Gradient_Conjugue_Tronque)

	# Tester l'algorithme des Régions de confiance avec PasdeCauchy | GCT
	tester_regions_de_confiance(affiche,Regions_De_Confiance)

	# Tester l'algorithme du Lagrangien Augmenté
	tester_lagrangien_augmente(affiche,Lagrangien_Augmente)
end

max_stars_repo_path: 2A/S7/UE Optimisation et RO/Optimisation/Optinum/test/runtests.jl

max_issues_repo_name: Sajid-Badr/ENSEEIHT

max_forks_repo_name: Sajid-Badr/ENSEEIHT


text:
#!/home/jeremy/julia/julia

# initialise the engine
include("Astellarn.jl")
using .Astellarn
# using Astellarn

# start the uci loop
# Base.@ccallable function julia_main()::Cint
    uci_main()
#     0
# end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
using DrWatson
@quickactivate "ICN benchmarks"

println(
"""
Currently active project is: $(projectname())

Path of active project: $(projectdir())

Have fun with your new project!

You can help us improve DrWatson by opening
issues on GitHub, submitting feature requests,
or even opening your own Pull Requests!
"""
)


text:
# To be executed for local deployment.
# Make sure you have LiveServer in your environment,
# as well as all the dependencies of the generator project.
using LiveServer

servedocs(literate = "", skip_dir = joinpath("docs", "src"))

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
# this file setup is less than ideal. We will refactor this with packing macros
# to automatically generate the converters and packed data types.
# The eval functions are also currently spread out, they should be concentrated here as well.
# The Pose2D and Pose3D types will most likely be packaged with the RoME package in the future.

# Active constraint types listed below
# -------------

# define the simple 1D odo
# TODO -- rework to use Distributions rather than Z and Cov
type Odo <: FunctorPairwise
    Zij::Array{Float64,2} # 0rotations, 1translation in each column
    Cov::Array{Float64,2}
    W::Array{Float64,1}
    Odo() = new()
    Odo(x...) = new(x[1], x[2], x[3])
end
# TODO -- only computing first node
function (odo::Odo)(res::Vector{Float64},
    idx::Int,
    meas::Tuple{Array{Float64,2}},
    p1::Array{Float64},
    p2::Array{Float64}  )

  res[1] = meas[1][1,idx] - (p2[1,idx] - p1[1,idx])
  nothing
end
function getSample(odo::Odo, N::Int=1)
  (reshape(rand(Distributions.Normal(odo.Zij[1,1], odo.Cov[1,1]), N ),1,N),)
end
# function getSample(odo::Odo, N::Int=1)
#   ret = zeros(1,N)
#   if size(odo.Zij,2) > 1
#     error("getSample(::Odo,::Int) can't handle multi-column at present")
#   end
#   for i in 1:N
#     ret[1,i] = odo.Cov[1]*randn()+odo.Zij[1]
#   end
#   # rand(Distributions.Normal(odo.Zij[1],odo.Cov[1]), N)'
#   return ret
# end
type PackedOdo <: PackedInferenceType
    vecZij::Array{Float64,1} # 0rotations, 1translation in each column
    dimz::Int
    vecCov::Array{Float64,1}
    dimc::Int
    W::Array{Float64,1}
    PackedOdo() = new()
    PackedOdo(x...) = new(x[1], x[2], x[3], x[4], x[5])
end
function convert(::Type{Odo}, d::PackedOdo)
  Zij = reshapeVec2Mat(d.vecZij,d.dimz)
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Odo(Zij, Cov, d.W)
end
function convert(::Type{PackedOdo}, d::Odo)
  v1 = d.Zij[:];
  v2 = d.Cov[:];
  return PackedOdo(v1,length(v1),
                    v2,length(v2),
                    d.W)
end

type OdoMM <: Pairwise
    Zij::Array{Float64,2} # 0rotations, 1translation in each column
    Cov::Array{Float64,2}
    W::Array{Float64,1}
end
function getSample(odo::OdoMM, N::Int=1)
  ret = zeros(1,N)
  if size(odo.Zij,2) > 1
    error("getSample(::OdoMM,::Int) can't handle multi-column at present")
  end
  for i in 1:N
    ret[1,i] = odo.Cov[1]*randn()+odo.Zij[1]
  end
  # rand(Distributions.Normal(odo.Zij[1],odo.Cov[1]), N)'
  return (ret,)
end

type Ranged <: FunctorPairwise
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    Ranged() = new()
    Ranged(x...) = new(x[1], x[2], x[3])
end
type PackedRanged <: PackedInferenceType
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    PackedRanged() = new()
    PackedRanged(x...) = new(x[1], x[2], x[3])
end
function convert(::Type{Ranged}, r::PackedRanged)
  return Ranged(r.Zij, r.Cov, r.W)
end
function convert(::Type{PackedRanged}, r::Ranged)
  return PackedRanged(r.Zij, r.Cov, r.W)
end
function (ra::Ranged)(res::Vector{Float64},
    idx::Int,
    meas::Tuple{Array{Float64,2}},
    p1::Array{Float64},
    l1::Array{Float64})

  res[1] = meas[1][1,idx] - abs(l1[1,idx] - p1[1,idx])
  nothing
end
function getSample(ra::Ranged, N::Int=1)
  ret = zeros(1,N)
  for i in 1:N
    ret[1,i] = ra.Cov[1]*randn()+ra.Zij[1]
  end
  # rand(Distributions.Normal(odo.Zij[1],odo.Cov[1]), N)'
  return (ret,)
end

type GenericMarginal <: FunctorPairwise
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    GenericMarginal() = new()
    GenericMarginal(a,b,c) = new(a,b,c)
end
type PackedGenericMarginal <: PackedInferenceType
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    PackedGenericMarginal() = new()
    PackedGenericMarginal(a,b,c) = new(a,b,c)
end
function convert(::Type{PackedGenericMarginal}, d::GenericMarginal)
  return PackedGenericMarginal(d.Zij, d.Cov, d.W)
end
function convert(::Type{GenericMarginal}, d::PackedGenericMarginal)
  return GenericMarginal(d.Zij, d.Cov, d.W)
end

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

type Obsv2 <: FunctorSingleton
    pts::Array{Float64,2}
    bws::Array{Float64,2}
    W::Array{Float64,1}
    Obsv2() = new()
    Obsv2(x...) = new(x[1], x[2], x[3])
end
type PackedObsv2 <: PackedInferenceType
    vecZij::Array{Float64,1} # 0rotations, 1translation in each column
    dimz::Int
    vecCov::Array{Float64,1}
    dimc::Int
    W::Array{Float64,1}
    PackedObsv2() = new()
    PackedObsv2(x...) = new(x[1],x[2],x[3],x[4],x[5])
end
function convert(::Type{Obsv2}, d::PackedObsv2)
  Zij = reshapeVec2Mat(d.vecZij,d.dimz)
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Obsv2(Zij, Cov, d.W)
end
function convert(::Type{PackedObsv2}, d::Obsv2)
  v1 = d.pts[:];
  v2 = d.bws[:];
  return PackedObsv2(v1,size(d.pts,1),
                    v2,size(d.bws,1),
                    d.W)
end
function getSample(obs::Obsv2, N::Int=1)
  pd = kde!(obs.pts, obs.bws[:,1])
  return (KernelDensityEstimate.sample(pd,N)[1],)
end


text:
"""
Custom DB collection for storing SimString data using base Dictionary `Dict`
"""
struct DictDB{
    T1<:FeatureExtractor,
    T2<:AbstractString,
    T3<:AbstractDict,
    T4<:AbstractDict,
    T5<:AbstractDict,
    } <: AbstractSimStringDB

    feature_extractor::T1                       # NGram feature extractor
    string_collection::Vector{T2}               # Collection of strings in the DB
    string_size_map::T3                         # Index map of feature sizes
    string_feature_map::T4                      # Index map of all features with associated strings and sizes
    lookup_cache::T5                            # Cache for lookup results
end

"""
    DictDB(x::CharacterNGrams)

Initialize a dict DB with additional containers and Metadata for CharacterNGrams

# Arguments
* `x`: CharacterNGrams object

# Example
```julia
db = DictDB(CharacterNGrams(2, " "))
```

# Returns
* `DictDB`: A DictDB object with additional containers and Metadata for CharacterNGrams
"""
function DictDB(x::CharacterNGrams)
    DictDB(
        x,
        String[],
        DefaultDict{Int, Set{String}}( () -> Set{String}() ),
        DefaultDict{ Int, DefaultOrderedDict{Tuple{String, Int64}, Set{String}}  }( () -> DefaultOrderedDict{Tuple{String, Int64}, Set{String} }(Set{String})),
        DefaultDict{ Int, DefaultDict{Tuple{String, Int64}, Set{String}}  }( () -> DefaultDict{Tuple{String, Int64}, Set{String}}(Set{String}))
    )
end

"""
    DictDB(x::WordNGrams)

Initialize a dict DB with additional containers and Metadata for WordNGrams

# Arguments
* `x`: WordNGrams object

# Example
```julia
db = DictDB(WordNGrams(2, " ", " "))
```

# Returns
* `DictDB`: A DictDB object with additional containers and Metadata for WordNGrams
"""
function DictDB(x::WordNGrams)
    DictDB(
        x,
        String[],
        DefaultDict{Int, Set{String}}( () -> Set{String}() ),
        DefaultDict{ Int, DefaultOrderedDict{Tuple{SubArray{SubString{String}}, Int}, Set{String}}  }( () -> DefaultOrderedDict{Tuple{SubArray{SubString{String}}, Int}, Set{String} }(Set{String})),
        DefaultDict{ Int, DefaultDict{Tuple{SubArray{SubString{String}}, Int}, Set{String}} }( () -> DefaultDict{Tuple{SubArray{SubString{String}}, Int}, Set{String}}(Set{String}))
    )
end

################################## DictDB UTIL Functions  ############################
"""
    describe_collection(db::DictDB)

Basic summary stats for the DB

# Arguments
* `db`: DictDB object

# Example
```julia
db = DictDB(CharacterNGrams(2, " "));
append!(db, ["foo", "bar", "fooo"]);
describe_collection(db)
(total_collection = 3, avg_size_ngrams = 4.5, total_ngrams = 13)

# Returns
* NamedTuples: Summary stats for the DB
```

"""
function describe_collection(db::DictDB)

# Total number of strings in collection
∑ = length(db.string_collection)

# Average size of ngram features
n = [x for x in keys(db.string_size_map)]
μ = sum(n) / length(n)

# Total number of ngram features
total_ngrams = 0
for i in values(db.string_feature_map)
    total_ngrams += length(i)
end

return (total_collection = ∑, avg_size_ngrams = μ, total_ngrams = total_ngrams)
end

"""
Pretty print summary stats for the DB
"""
function Base.show(io::IO, x::DictDB)
    metrics = describe_collection(x)
    println(io, "DictDB($(x.feature_extractor))")
    println(io, "Total collection: ", metrics.total_collection)
    println(io, "Average number of ngram features: ", metrics.avg_size_ngrams)
    println(io, "Total number of ngram features: ", metrics.total_ngrams)
end

"""
Internal function to lookup feature sets by size and feature
"""
function lookup_feature_set_by_size_feature(db::DictDB, size, feature)
    if !haskey(db.lookup_cache[size], feature)
        db.lookup_cache[size][feature] = get(db.string_feature_map[size], feature, Set{String}())
    end
    return db.lookup_cache[size][feature]
end


text:
# wrappers of the fathom c library
# currently only Linux support for Astellarn

const TB_RESULT_FAILED = 0xFFFFFFFF

# constants defined in fathom
const TB_LOSS = 0
const TB_BLESSED_LOSS = 1
const TB_DRAW = 2
const TB_CURSED_WIN = 3
const TB_WIN = 4

const TB_PROMOTES_NONE = 0
const TB_PROMOTES_QUEEN = 1
const TB_PROMOTES_ROOK = 2
const TB_PROMOTES_BISHOP = 3
const TB_PROMOTES_KNIGHT = 4

const TB_RESULT_WDL_MASK = 0x0000000F
const TB_RESULT_TO_MASK = 0x000003F0
const TB_RESULT_FROM_MASK = 0x0000FC00
const TB_RESULT_PROMOTES_MASK = 0x00070000
const TB_RESULT_EP_MASK = 0x00080000
const TB_RESULT_DTZ_MASK = 0xFFF00000
const TB_RESULT_WDL_SHIFT = 0
const TB_RESULT_TO_SHIFT = 4
const TB_RESULT_FROM_SHIFT = 10
const TB_RESULT_PROMOTES_SHIFT = 16
const TB_RESULT_EP_SHIFT = 19
const TB_RESULT_DTZ_SHIFT = 20

function TB_GET_WDL(res::UInt32)
    Int((res & TB_RESULT_WDL_MASK) >> TB_RESULT_WDL_SHIFT)
end

function TB_GET_TO(res::UInt32)
    Int(((res & TB_RESULT_TO_MASK) >> TB_RESULT_TO_SHIFT) + 1)
end

function TB_GET_FROM(res::UInt32)
    Int(((res & TB_RESULT_FROM_MASK) >> TB_RESULT_FROM_SHIFT) + 1)
end

function TB_GET_PROMOTES(res::UInt32)
    Int((res & TB_RESULT_PROMOTES_MASK) >> TB_RESULT_PROMOTES_SHIFT)
end

# check output of this function
function TB_GET_EP(res::UInt32)
    Int(((res & TB_RESULT_EP_MASK) >> TB_RESULT_EP_SHIFT) + 1)
end

function TB_GET_DTZ(res::UInt32)
    Int((res & TB_RESULT_DTZ_MASK) >> TB_RESULT_DTZ_SHIFT)
end

function tb_init(syzygypath::S)::Bool where S <: Union{String, SubString}
    return ccall((:tb_init, FATHOM_PATH), UInt8, (Ptr{UInt8},), syzygypath)
end

function tb_free()
    return ccall((:tb_free, FATHOM_PATH), Cvoid, ())
end

# We can't call static inline function :tb_probe_wdl, so we use :tb_probe_wdl_impl as a workaround.
# Because of the above, care is taken to ensure we conform to syzygy rules.
function tb_probe_wdl(board::Board)::UInt32
    if iszero(board.enpass) === false
        return TB_RESULT_FAILED
    elseif iszero(board.castling) === false
        return TB_RESULT_FAILED
    elseif iszero(board.halfmovecount) === false
        return TB_RESULT_FAILED
    else
        return ccall((:tb_probe_wdl_impl, FATHOM_PATH), Cuint, (UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt32, UInt8,),
            board[WHITE].val, board[BLACK].val, board[KING].val, board[QUEEN].val, board[ROOK].val, board[BISHOP].val, board[KNIGHT].val, board[PAWN].val,
            0, board.turn == WHITE ? 1 : 0)
    end
end

# We can't call static inline function :tb_probe_root, so we use :tb_probe_root_impl as a workaround.
# Because of the above, care is taken to ensure we conform to syzygy rules.
function tb_probe_root(board::Board)::UInt32
    if iszero(board.castling) === false
        return TB_RESULT_FAILED
    else
        return ccall((:tb_probe_root_impl, FATHOM_PATH), Cuint, (UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt32, UInt32, UInt8, UInt32, ),
            board[WHITE].val, board[BLACK].val, board[KING].val, board[QUEEN].val, board[ROOK].val, board[BISHOP].val, board[KNIGHT].val, board[PAWN].val,
            board.halfmovecount, board.enpass, board.turn == WHITE ? 1 : 0, 0)
    end
end

# Called at the root node if a TB result is returned
function interpret_syzygy(thread::Thread, res::UInt32)
    thread.ss.tbhits += 1
    _eval = TB_GET_WDL(res)
    if iszero(_eval)
        eval = -MATE
    elseif 1 <= _eval <= 3 # blessed / cursed loss and wins are draws
        eval = 0
    else
        eval = MATE
    end
    move_from = TB_GET_FROM(res)
    move_to = TB_GET_TO(res)
    promotion = TB_GET_PROMOTES(res)
    if promotion !== TB_PROMOTES_NONE
        if promotion == TB_PROMOTES_QUEEN
            push!(thread.pv[1], Move(move_from, move_to, __QUEEN_PROMO))
            return eval
        elseif promotion == TB_PROMOTES_ROOK
            push!(thread.pv[1], Move(move_from, move_to, __ROOK_PROMO))
            return eval
        elseif promotion == TB_PROMOTES_BISHOP
            push!(thread.pv[1], Move(move_from, move_to, __BISHOP_PROMO))
            return eval
        elseif promotion == TB_PROMOTES_KNIGHT
            push!(thread.pv[1], Move(move_from, move_to, __KNIGHT_PROMO))
            return eval
        end
    else
        clear!(thread.pv[1])
        push!(thread.pv[1], Move(move_from, move_to, __NORMAL_MOVE))
        return eval
    end
end


text:
using Base: Bool
using LinearAlgebra
using SparseArrays

# functions

function unity(i,n)
    E = spzeros(Bool,n,n)
    E[i] = true
    return E
end

# xs = [rand(3,4) for _ in 1:5]
# reduce((x,y) -> cat(x,y,dims=3), A)
# idea: https://github.com/JuliaLang/julia/pull/37196
# reduce(cat(dims=4), [ones(2,3) for _ in 1:5]) |> size
#`reduce(vcat, [rand(2,3) for _ = 1:4]): A[i][a,b] -> A[(i-1)*len_a + a, b]`
#which is the same as restructuring to` B[i,a][b] via B = [A[i][a,:] for i = 1:dim_i, a = 1:dim_a] `
#and then taking `reduce(hcat,B'[:]')'`
#`reduce(hcat, [rand(2,3) for _ = 1:4]): A[i][a,b] -> A[a, (i-1)*len_b + b]`

function partial_trace(dim_1, dim_2, traced_system)
# traces out the system indicated by traced_system
# for C = kron(A,B), tracing over the first system (traced_system = 1);
# yields: reshape(p1*C[:], 4,4)/tr(A)  == B
#
# in colex notation:
# traced_system == 1 represents a trace over index 2 and 4 and
# traced_system == 2 represents a trace over index 1 and 3 and
    if traced_system == 1
        T = [kron(sparse(I, dim_1,dim_1), unity(i,dim_2))[:] for i = 1:dim_2^2]
        col = reduce(vcat,[findnz(T[i])[1] for i = 1:dim_2^2])
        row = kron(1:dim_2^2, ones(Int64,dim_1))
        return sparse(row, col, trues(dim_2^2*dim_1))
    else
        T = [kron(unity(i, dim_1), sparse(I, dim_2, dim_2))[:] for i = 1:dim_1^2]
        col = reduce(vcat, [findnz(T[i])[1] for i = 1:dim_1^2])
        row = kron(1:dim_1^2, ones(Int64,dim_2))
        return sparse(row, col, trues(dim_1^2*dim_2))
        #return    Tr_B = [reshape(kron(reshape(unity(i,dim_1), dim_1, dim_1), I(dim_2)), 1, dim_2^2 * dim_1^2)[1,j] for i = 1:dim_1^2, j =1:dim_1^2*dim_2^2]
    end
end

function choi(A, dim_1)
mm, nn = size(A)
n = mm ÷ dim_1
return reshape(permutedims(reshape(A,n,dim_1,n,dim_1),(1,3,2,4)), mm,mm)
end

function choi_rect(A,dim_1, dim_2, dim_3, dim_4)
return reshape(permutedims(reshape(A, dim_1, dim_2, dim_3, dim_4), (1,3,2,4)), dim_1*dim_3, dim_2*dim_4)
end

function ptr(ρ, dim_1, traced_system = 1)
    # get partial trace superoperator
    N,M = size(ρ)
    dim_2 = N÷dim_1
    PT = partial_trace(dim_1, dim_2, traced_system)
    new_dim = traced_system == 1 ? dim_2 : dim_1
    # apply superoperator to vectorized density matrix and reshape to new dimensions
    σ = reshape(PT * ρ[:], new_dim, new_dim)
    return σ
end

### apply split
function factorize(A, dim_1, dim_2)
    σ = prt(A, dim_1, dim_2, 2)
    ρ = prt(A, dim_1, dim_2, 1)
    A_tensor_prod = kron(σ, ρ)
    rho_c = A - A_tensor_prod
    return A_tensor_prod, rho_c
end

### density matrix

function density_matrix(n, neg = 0, real_valued = false)
    #trace 1 remains preserved!
    dia = rand(n)
    if neg > 0
        dia[rand(1:n, neg)] .*= -1
    end
    dia = dia/sum(dia)
    if real_valued
        temp = rand(n,n)
    else
        temp = rand(n,n) + rand(n,n) * im
    end
    ew, ev = eigen(Hermitian(temp'*temp))
    d_mat = Hermitian(ev* Diagonal(dia) * ev')
    return d_mat
end

function choi2kraus(C)
    dim_1,dim_2 = size(C)
    N= Int64(sqrt(dim_1))
    K = [zeros(ComplexF64, N,N) for _=1:N^2]
    ew,ev = eigen(C)
    for k = 1:N^2
        K[k] = reshape(conj.(ev[:,k]),N,N) .* sqrt(abs(ew[k]))
        if ew[k] < 0
            print("Warning, negative eigenvalue, k = " * string(k))
        end
    end
    return K
end

function kraus2choi(K)
    dim1,dim2 = size(K[1])
    C = zeros(ComplexF64, dim1^2,dim1^2)
    for k = 1:length(K)
        C += conj.(K[k])[:] * transpose(K[k][:])
    end
    return C
end
### ======================== functions end =======================

dim_1 = 2
dim_2 = 3
m = 7
n = 3

A = rand(m,m)
B = rand(n,n)

C = kron(A,B)

A_sup = kron(A[:], I(n^2))
C_sup = choi(C,m)

sum(abs.(A_sup * B[:] - C_sup[:]))

A_sup * B[:] - C[:]
reshape(A_sup * B[:],n*m,n*m) - C_sup
sum(abs.(ans[:]))

# create superoperator without applying Choi on the right hand side

G = reshape(permutedims(reshape(kron(A[:], I(n^2)), (n,n,m,m,n^2)), (1,3,2,4,5)), n^2*m^2,n^2)

reshape(G*B[:], m*n, m*n)
C
sum(abs.(G * B[:] - C[:]))

kron_fac_B = reshape(permutedims(reshape(kron(I(m^2), B[:]), (n,n,m,m,m^2)), (1,3,2,4,5)), n^2*m^2, m^2)
norm(reshape(kron_fac_B * A[:], n*m, n*m) - C)

function I_rect(n1,n2)
    if n1 >= n2
        temp = I(n1)
        temp = temp[:,1:n2]
    else
        temp = I(n2)
        temp = temp[1:n1,:]
    end
    return temp
end

m1 = 3
m2 = 4

n1 = 5
n2 = 6

A = rand(m1,m2)
B = rand(n1,n2)

C = kron(A,B)

#G = reshape(permutedims(reshape(kron(A[:], I(n1*n2)), (n1,n2,m1,m2,n1*n2)), (1,3,2,4,5)), n1*n2*m1*m2,n1*n2)

S_left = reshape(permutedims(reshape(kron(I(n1*n2), A[:]), (m1,m2,n1,n2,n1*n2)), (3,1,4,2,5)), prod([m1,m2,n1,n2]), n1*n2)
norm(S_left * B[:] - C[:])

S_right = reshape(permutedims(reshape(kron(B[:], I(m1*m2)), (m1,m2,n1,n2,m1*m2)), (3,1,4,2,5)), prod([m1,m2,n1,n2]), m1*m2)
norm(S_right * A[:] - C[:])

norm(kron(A[:], I(n1*n2)) * B[:] - choi_rect(C, n1,m1, n2, m2)[:])

norm(kron(I(m1*m2), B[:]) * A[:] - choi_rect(C,n1,m1,n2,m2)[:])

norm(kron(A'[:], I(n1*n2)) * B'[:] - C'[:])

# elegant version to show right hand side as a Choi matrix for first and second factor:
norm(kron(I(n1*n2) ,A[:]) * B[:] - choi_rect(C, n1, m1, n2, m2)'[:])
norm(kron(B[:], I(m1*m2)) * A[:] - choi_rect(C, n1, m1, n2, m2)'[:])

# check embeddings

A = rand(2)
B = rand(3)
C = rand(4)
D = rand(5)

Q = kron(A, B')
kron(kron(A, B'), ones(4,5))

kron(kron(A,B'), kron(C,D')) - kron(kron(A, B'), I(4)) * kron(I(4), kron(C,D'))

Q = kron(A, C')
kron(kron(A, ones(3)'), kron(C, ones(5)'))- kron(Q'[:], ones(3*5)')

Q = kron(A, D')
kron(kron(A, ones(3)'), kron(ones(4), D'))-kron(kron(ones(3)', Q), ones(4))

Q = kron(B, C')
kron(kron(ones(2), B'), kron(C, ones(5)')) - kron(kron(ones(2),Q'), ones(5)')

Q = kron(B, D')
kron(kron(ones(2), B'), kron(ones(4), D'))-kron(ones(2*4), Q[:]')

# multiplication table for nxn unity base matrices

n = 4
A = zeros(Int64,n^2,n^2)
function unity(i,n)
    E = spzeros(Bool,n,n)
    E[i] = true
    return E
end

for i = 1:n^2, j=1:n^2
    L = findall((unity(i,n)*unity(j,n))[:] .> 0)
    if isempty(L)
        A[i,j] = 0
    else
        A[i,j] = L[1]
    end
end
A

# try to check \sum_{i != j} a_ii a_jj - abs(a_ij)^2 == sum_{i \neq j} \lambda_i lambda_j

A = density_matrix(10,0);
#Diagonal(real(A))
ew,ev = eigen(A);
L1 = [ A[i,i]*A[j,j] - abs(A[i,j])^2 for i = 1:5, j = 1:5];
L2 = [ew[i] * ew[j] * (i != j) for i = 1:5, j=1:5];

sum(real(L1) .< 0)
real(ew)'

# can every density matrix be represented as the sum of tensor products of valid density matrices?

A = density_matrix(4)
B = density_matrix(3)

ew, ev = eigen(kron(A,B))
sum(ew .< 0)

a = rand(4)

# the Kraus decomposition for initially correlated dynamical maps:
A6 = density_matrix(6)
A2 = density_matrix(2)
A3 = density_matrix(3)
A4 = density_matrix(4)

A = density_matrix(9)

n_runs = 80

tr_weight = zeros(n_runs)

A =  kron(A3, A4)
for k = 1:n_runs
    T1, C1 = factorize(A./tr(A), 2,6)

    A, M1 = split(C1)
    tr_weight[k] = tr(A)
end
tr_weight
sum(cumprod(tr_weight))

# entangled state

A = [1 0 0 1; 0 0 0 0; 0 0 0 0; 1 0 0 1]/2
t,c = factorize(A/tr(A),2,2)
p,m = split(c)

ω = tr(p)
p = p/tr(p)
m = m/tr(m)
round)
function split(rho)
    ew, ev = eigen(rho)
    rho_p = Hermitian(ev[:,ew .> 0] * Diagonal(ew[ew.> 0]) * ev[:,ew .> 0]')
    rho_m = Hermitian(ev[:,ew .< 0] * Diagonal(-ew[ew.< 0]) * ev[:,ew .< 0]')
    return rho_p, rho_m
end

P1 = partial_trace(2,3,1)
P2 = partial_trace(2,3,2)

A1 = reshape(P2*A[:],2,2)
A2 = reshape(P1*A[:],3,3)

rho_c = A - kron(A1,A2)
rho_c = Hermitian(rho_c)

ew,ev = eigen(rho_c)

rho_p = Hermitian(ev[:,ew .> 0] * Diagonal(ew[ew.> 0]) * ev[:,ew .> 0]')
rho_m = Hermitian(ev[:,ew .< 0] * Diagonal(-ew[ew.< 0]) * ev[:,ew .< 0]')

# now decompose both again, maybe it does converge

A = rand(4,4)
B = rand(4,4)

kron(A,B)' - kron(A', B')

# Choi to Kraus representation

A = density_matrix(16) * 4

ew, ev = eigen(A)
K = [zeros(ComplexF64,4,4) for _ in 1:16]
for k = 1:16
    K[k] = reshape(ev[:,k], 4,4) .* sqrt(ew[k])
end

Test = zeros(ComplexF64,4,4)
for k=1:16
    Test = Test + K[k] * K[k]'
end
Test
tr(Test)
B = density_matrix(4)
B2 = zeros(ComplexF64,4,4)
for k = 1:16
    B2 += K[k] * B * K[k]'
end
tr(B)
tr(B2)

A = rand(4,4)
A = A

exp(-im*A) * exp(-im*A)'

# general way to construct a valid quantum operation:
# find A \in C^{d^2 x d^2}
# with A^\dag = A and A>0
# then A2 = A + (I - P_13 A) \kron I/d is a valid Choi matrix

# create choi matrix directly fulfilling the trace condition

A = density_matrix(16)*4
P13 = partial_trace(4,4,2) #tracing out the second system - flip to switch to colexicographic
P24 = partial_trace(4,4,1)

A2 = A + kron(I(4) - reshape(P13*A[:], 4,4), I(4)/4)

ew, ev = eigen(A2); ew
round(sum(abs.(ew)) - tr(A2), digits=14)

reshape(P13*A2[:], 4,4)
# unital map: Tr_24 yields I(N)
A3 = A2 + kron(I(4)/4, I(4) - reshape(P24*A[:],4,4))

round.(reshape(P13*A3[:], 4,4), digits = 14)

round.(reshape(P24*A3[:], 4,4), digits = 14)

# get Kraus operators from A2 from eigenvectors
#(complex conjugation can be ignored since C is hermitian):
K = [zeros(ComplexF64, 4,4) for _=1:16]
for k = 1:16
    K[k] = reshape(ev[:,k],4,4) * sqrt(ew[k])

end
K = choi2kraus(A2)
id_test = zeros(ComplexF64,4,4)
for k = 1:16
    id_test += K[k]' * K[k]
end

id_test

C = kraus2choi(K)
sum(abs.(C - (A2)))
C
A2

#Ginibre matrix NxN
#G = randn(N,N)
# 1/(2 pi )^(N^2/2) * exp( -1/2 Tr(G * transpose(G)))

N= 1000
G = randn(N,N)
1/(2*pi)^(N^2/2) * exp(- 1/2 * tr(G * transpose(G)))

# create random Kraus operators from large unitary matrix
N = 4
A = rand(N^2*N, N^2*N)
ew, ev = eigen(Hermitian(A))
K = [zeros(ComplexF64, 4,4) for _=1:16]
for k=1:16
    K[k] = ev[1:4, (1:4) .+ (k-1)*4]'
end

# get the choi matrix from the Kraus operators
C = zeros(ComplexF64, 16,16)

for k = 1:16
    C = C + K[k]'[:] * K[k]'[:]'
end
ew, ev = eigen(C);
tr(C) - sum(ew)
sum(abs.(C - C'))

# get linear superoperator
D = choi(C,4)

p1 = partial_trace(4,4,1)

Test  = zeros(4,4)
reshape(sum(p1, dims = 1),4,4)

A = randn(12,16) + randn(12,16) * 0.1*im
C_ = A'*A/tr(A'*A)*4
w,v = eigen(C_); w

#C_ = density_matrix(16) * 4
C = (C_ + kron(I(4) - ptr(C_, 4, 2), I(4)/4))
tr(C)
ew, ev = eigen(C);ew

ptr(C,4,2)
reshape(P13 * C[:],4,4)

# Kraus operators of partial transpose

#transpose operator
T = [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]

C = choi(T,2)
ew, ev = eigen(C)

K = choi2kraus(C)
A = density_matrix(2)
temp = zeros(ComplexF64,2,2)
for k = 1:4
    temp += K[k] * A * K[k]'
    if k == 1
        temp *=-1
    end
end
temp - A
A
Id_test  = zeros(ComplexF64,2,2)
for k = 1:4
 Id_test += K[k]' * K[k]
end
Id_test

K[1]' * K[1]

A = rand(4,4)

B = rand(3,3)

# tensor product as superoperator

# L(X) = A ⊗ X

A = rand(2,2)
X = rand(3,3)
# aim:
L * X[:] == kron(A,X)[:]

CL = kron(A[:], I(9))
L = reshape(permutedims(reshape(CL, 3,3,2,2,9),[1,3,2,4,5]), 36,9)
L * X[:]
kron(A,X)[:]

L2 = kron(kron(A,I(3))[:], I(3)[:]')
L

i = 16

n = [2,3,2,3]

a = [(i-1) % n[1]+1,
((i-1) % prod(n[1:2]) ÷ n[1]) + 1,
(i-1) % prod(n[1:3]) ÷ prod(n[1:2]) + 1,
(i-1) ÷ prod(n[1:3]) + 1]

a[1] + n[1] * (a[2]-1 + n[2] * (a[3] - 1 + n[3] * (a[4]-1)))

function choi_index(i, n)
    return  (i.-1) .% n[1] .+ 1 .+
    n[1] .* ((i.-1) .% prod(n[1:3]) .÷ prod(n[1:2]) .+
    n[2] .* ((i.-1) .% prod(n[1:2]) .÷ n[1] .+
    n[3] .* ((i.-1) .÷ prod(n[1:3]))))
end
T = reshape(1:36,2,3,2,3)


text:
module Neet

# provides the @unimpl macro
include("macros.jl")

# provides the Abstract interfaces
include("core.jl")

# provides the state space types
include("statespace.jl")

end


text:
using VisualRegressionTests
using Base.Test

# load an image for testing
import Images, TestImages
img = TestImages.testimage("cameraman")

# save a temporary copy
reffn = tempname() * ".png"
Images.save(reffn, img)

# this is the test function, which saves an image to the given location
func = fn -> Images.save(fn, img)

# do the test
result = test_images(VisualTest(func, reffn))

@test isa(result, VisualTestResult)
@test result.status == EXACT_MATCH
@test result.err == nothing
@test result.diff == 0


text:
##########   I/O   ###########

import .LibMagick
import .LibOSXNative

abstract ImageFileType

# Database of image file extension, magic bytes, and file type stubs
# Add new items to this database using add_image_file_format (below)
#
# Filename extensions do not uniquely specify the image type, so we
# have extensions pointing to a list of candidate formats. On reading,
# this list is really just a hint to improve efficiency: in the end,
# it's the set of magic bytes in the file that determine the file
# format.
fileext = Dict{ByteString, Vector{Int}}()
filemagic = Array(Vector{Uint8}, 0)
filetype = Array(Any, 0)
filesrcloaded = Array(Bool, 0)
filesrc = Array(String, 0)

function _loadformat(index::Int)
    filename = joinpath("ioformats", filesrc[index])
    if !isfile(filename)
        filename = joinpath(Pkg.dir(), "Images", "src", "ioformats", filesrc[index])
    end
    include(filename)
    filesrcloaded[index] = true
end
function loadformat{FileType<:ImageFileType}(::Type{FileType})
    indx = find(filetype .== FileType)
    if length(indx) == 1
        _loadformat(indx[1])
    elseif isempty(indx)
        error("File format ", FileType, " not found")
    else
        error("File format ", FileType, " is entered multiple times")
    end
    nothing
end

function add_image_file_format{ImageType<:ImageFileType}(ext::ByteString, magic::Vector{Uint8}, ::Type{ImageType}, filecode::ASCIIString)
    # Check to see whether these magic bytes are already in the database
    for i in 1:length(filemagic)
        if magic == filemagic[i]
            fileext[ext] = [i]  # update/add extension lookup
            filetype[i] = ImageType
            return
        end
    end
    # Add a new entry on the list
    push!(filemagic, magic)
    len = length(filemagic)
    push!(filetype, ImageType)
    if haskey(fileext, ext)
        push!(fileext[ext], len)
    else
        fileext[ext] = [len]
    end
    push!(filesrcloaded, isempty(filecode))
    push!(filesrc, filecode)
end
add_image_file_format{ImageType<:ImageFileType}(ext::ByteString, magic::Vector{Uint8}, ::Type{ImageType}) = add_image_file_format(ext, magic, ImageType, "")
add_image_file_format{ImageType<:ImageFileType}(ext::ByteString, ::Type{ImageType}, filecode::ASCIIString) = add_image_file_format(ext, b"", ImageType, filecode)

# Define our fallback file formats now, because we need them in generic imread.
# This has no extension (and is not added to the database), because it is always used as a stream.
type ImageMagick <: ImageFileType end
type OSXNative <: ImageFileType end

function imread(filename::String)
    _, ext = splitext(filename)
    ext = lowercase(ext)

    img = open(filename, "r") do stream
        # Tries to read the image using a set of potential type candidates.
        # Returns the image if successful, `nothing` else.
        function tryread(candidates)
            if (index = image_decode_magic(stream, candidates)) > 0
                # Position to end of this type's magic bytes
                seek(stream, length(filemagic[index]))
                if !filesrcloaded[index]
                    _loadformat(index)
                end
                imread(stream, filetype[index])
            end
        end

        # Use the extension as a hint to determine file type
        if haskey(fileext, ext) && (img = tryread(fileext[ext])) != nothing
            return img
        end

        # Extension wasn't helpful, look at all known magic bytes
        if (img = tryread(1:length(filemagic))) != nothing
            return img
        end
    end

    if img != nothing
        return img
    end

    @osx_only begin
        img = imread(filename, OSXNative)
        if img != nothing
            return img
        end
    end

    # There are no registered readers for this type. Try using ImageMagick if available.
    if have_imagemagick
        return imread(filename, ImageMagick)
    else
        error("Do not know how to read file ", filename, ". Is ImageMagick installed properly? See README.")
    end
end

imread{C<:ColorValue}(filename::String, ::Type{C}) = imread(filename, ImageMagick, C)

# Identify via magic bytes
function image_decode_magic{S<:IO}(stream::S, candidates::AbstractVector{Int})
    maxlen = 0
    for i in candidates
        len = length(filemagic[i])
        maxlen = (len > maxlen) ? len : maxlen
    end
    # If there are no magic bytes, simply use the file extension
    if maxlen == 0 && length(candidates) == 1
        return candidates[1]
    end

    magicbuf = zeros(Uint8, maxlen)
    for i=1:maxlen
        if eof(stream) break end
        magicbuf[i] = read(stream, Uint8)
    end
    for i in candidates
        if length(filemagic[i]) == 0
            continue
        end
        ret = ccall(:memcmp, Int32, (Ptr{Uint8}, Ptr{Uint8}, Int), magicbuf, filemagic[i], min(length(filemagic[i]), length(magicbuf)))
        if ret == 0
            return i
        end
    end
    return -1
end

function imwrite(img, filename::String; kwargs...)
    _, ext = splitext(filename)
    ext = lowercase(ext)
    if haskey(fileext, ext)
        # Write using specific format
        candidates = fileext[ext]
        index = candidates[1]  # TODO?: use options, don't default to first
        if !filesrcloaded[index]
            _loadformat(index)
        end
        imwrite(img, filename, filetype[index]; kwargs...)
    elseif have_imagemagick
        # Fall back on ImageMagick
        imwrite(img, filename, ImageMagick; kwargs...)
    else
        error("Do not know how to write file ", filename)
    end
end

function imwrite{T<:ImageFileType}(img, filename::String, ::Type{T}; kwargs...)
    open(filename, "w") do s
        imwrite(img, s, T; kwargs...)
    end
end

# only mime writeable to PNG if 2D (used by IJulia for example)
import Base.mimewritable
Base.mimewritable(::MIME"image/png", img::AbstractImage) = sdims(img) == 2 && timedim(img) == 0

function _writemime(stream::IO, ::MIME"image/png", img::AbstractImage; scalei=scaleinfo_uint(img))
    assert2d(img)
    if isa(img, AbstractImageIndexed)
        # For now, convert to direct
        img = convert(Image, img)
    end
    A = data(img)
    nc = ncolorelem(img)
    npix = length(A)/nc
    while npix > 1e6
        A = restrict(A, coords_spatial(img))
        npix = length(A)/nc
    end
    if eltype(A) != eltype(img)
        A = truncround(eltype(img), A)
    end
    wand = image2wand(share(img, A), scalei)
#     LibMagick.setimagecompression(wand, LibMagick.NoCompression)
    blob = LibMagick.getblob(wand, "png")
    write(stream, blob)
end
writemime(stream::IO, mime::MIME"image/png", img::AbstractImage{RGB}; scalei = scaleinfo_uint(img)) = _writemime(stream, mime, img, scalei=scalei)
writemime{C<:ColorValue}(stream::IO, mime::MIME"image/png", img::AbstractImage{C}; kwargs...) = writemime(stream, mime, convert(Image{RGB}, img); kwargs...)
writemime(stream::IO, mime::MIME"image/png", img::AbstractImage; scalei = scaleinfo_uint(img)) = _writemime(stream, mime, img, scalei=scalei)

#### Implementation of specific formats ####

#### OSX Native readers from CoreGraphics

imread(filename::String, ::Type{OSXNative}) = LibOSXNative.imread(filename)

#### ImageMagick library

typedict = [1 => Uint8, 8 => Uint8, 16 => Uint16, 32 => Uint32]

function imread(filename::String, ::Type{ImageMagick})
    wand = LibMagick.MagickWand()
    LibMagick.readimage(wand, filename)
    LibMagick.resetiterator(wand)
    imtype = LibMagick.getimagetype(wand)
    # Determine what we need to know about the image format
    sz = size(wand)
    n = LibMagick.getnumberimages(wand)
    if n > 1
        sz = tuple(sz..., n)
    end
    havealpha = LibMagick.getimagealphachannel(wand)
    cs = LibMagick.getimagecolorspace(wand)
    nc, cs = LibMagick.nchannels(imtype, cs, havealpha)
    depth = LibMagick.getimagedepth(wand)
    T = typedict[iceil(depth/8)*8]
    # Allocate the buffer and get the pixel data
    buf = (nc > 1) ? Array(T, nc, sz...) : Array(T, sz...)
    LibMagick.exportimagepixels!(buf, wand, cs)
    # Set up the properties
#     spatialorder = (n > 1) ? ["x", "y", "z"] : ["x", "y"]
    prop = ["colorspace" => cs, "spatialorder" => ["x", "y"], "limits" => (zero(T), typemax(T))]
    if nc > 1
        prop["colordim"] = 1
    end
    if n > 1
        prop["timedim"] = ndims(buf)
    end
    Image(buf, prop)
end

imread{C<:ColorValue}(filename::String, ::Type{ImageMagick}, ::Type{C}) = convert(Image{C}, imread(filename, ImageMagick))

function imwrite(img, filename::String, ::Type{ImageMagick}; scalei = scaleinfo_uint(img))
    wand = image2wand(img, scalei)
    LibMagick.writeimage(wand, filename)
end

function image2wand(img, scalei)
    if isa(img, AbstractImageIndexed)
        # For now, convert to direct
        img = convert(Image, img)
    end
    imgw = scale(scalei, img)
    imgw = permutedims_canonical(imgw)
    have_color = colordim(imgw)!=0
    if ndims(imgw) > 3+have_color
        error("At most 3 dimensions are supported")
    end
    wand = LibMagick.MagickWand()
    LibMagick.constituteimage(to_explicit(to_contiguous(data(imgw))), wand, colorspace(img))
    LibMagick.resetiterator(wand)
    wand
end

# Make the data contiguous in memory, this is necessary for imagemagick since it doesn't handle stride.
to_contiguous(A::AbstractArray) = A
to_contiguous(A::SubArray) = copy(A)

to_explicit(A::AbstractArray) = A
to_explicit(A::AbstractArray{RGB8}) = reinterpret(Uint8, A, tuple(3, size(A)...))

# Write grayscale values in horizontal-major order
function writegray(stream, img, scalei::ScaleInfo)
    assert2d(img)
    xfirst = isxfirst(img)
    firstindex, spsz, spstride, csz, cstride = iterate_spatial(img)
    isz, jsz = spsz
    istride, jstride = spstride
    A = parent(img)
    if xfirst
        for j = 1:jsz
            k = firstindex + (j-1)*jstride
            for i = 0:istride:(isz-1)*istride
                write(stream, scale(scalei, A[k+i]))
            end
        end
    else
        for i = 1:isz
            for j = 1:jsz
                k = firstindex+(i-1)*istride+(j-1)*jstride
                write(stream, scale(scalei,A[k]))
            end
        end
    end
end

# Write RGB or RGBA values in horizontal-major order
# The A value is written if present, otherwise skipped
function writecolor(stream, img, scalei::ScaleInfo)
    assert2d(img)
    xfirst = isxfirst(img)
    firstindex, spsz, spstride, csz, cstride = iterate_spatial(img)
    isz, jsz = spsz
    istride, jstride = spstride
    A = parent(img)
    cs = colorspace(img)
    if cs == "Gray"
        if xfirst
            for j = 1:jsz
                k = firstindex + (j-1)*jstride
                for i = 0:istride:(isz-1)*istride
                    v = scale(scalei, A[k+i])
                    write(stream, v)
                    write(stream, v)
                    write(stream, v)
                end
            end
        else
            for i = 1:isz
                for j = 1:jsz
                    k = firstindex+(i-1)*istride+(j-1)*jstride
                    v = scale(scalei, A[k])
                    write(stream, v)
                    write(stream, v)
                    write(stream, v)
                end
            end
        end
    elseif cs == "RGB"
        if xfirst
            for j = 1:jsz
                k = firstindex + (j-1)*jstride
                for i = 0:istride:(isz-1)*istride
                    ki = k+i
                    write(stream, scale(scalei, A[ki]))
                    write(stream, scale(scalei, A[ki+cstride]))
                    write(stream, scale(scalei, A[ki+2cstride]))
                end
            end
        else
            for i = 1:isz
                for j = 1:jsz
                    k = firstindex+(i-1)*istride+(j-1)*jstride
                    write(stream, scale(scalei, A[k]))
                    write(stream, scale(scalei, A[k+cstride]))
                    write(stream, scale(scalei, A[k+2cstride]))
                end
            end
        end
    elseif cs == "ARGB"
        if xfirst
            for j = 1:jsz
                k = firstindex + (j-1)*jstride
                for i = 0:istride:(isz-1)*istride
                    ki = k+i
                    write(stream,scale(scalei, A[ki+cstride]))
                    write(stream,scale(scalei, A[ki+2cstride]))
                    write(stream,scale(scalei, A[ki+3cstride]))
                    write(stream,scale(scalei, A[ki]))
                end
            end
        else
            for i = 1:isz
                for j = 1:jsz
                    k = firstindex+(i-1)*istride+(j-1)*jstride
                    write(stream,scale(scalei, A[k+cstride]))
                    write(stream,scale(scalei, A[k+2cstride]))
                    write(stream,scale(scalei, A[k+3cstride]))
                    write(stream,scale(scalei, A[k]))
                end
            end
        end
    elseif cs == "RGBA"
        if xfirst
            for j = 1:jsz
                k = firstindex + (j-1)*jstride
                for i = 0:istride:(isz-1)*istride
                    ki = k+i
                    write(stream,scale(scalei, A[ki]))
                    write(stream,scale(scalei, A[ki+cstride]))
                    write(stream,scale(scalei, A[ki+2cstride]))
                    write(stream,scale(scalei, A[ki+3cstride]))
                end
            end
        else
            for i = 1:isz
                for j = 1:jsz
                    k = firstindex+(i-1)*istride+(j-1)*jstride
                    write(stream,scale(scalei, A[k]))
                    write(stream,scale(scalei, A[k+cstride]))
                    write(stream,scale(scalei, A[k+2cstride]))
                    write(stream,scale(scalei, A[k+3cstride]))
                end
            end
        end
    elseif cs == "RGB24"
        if xfirst
            for j = 1:jsz
                k = firstindex + (j-1)*jstride
                for i = 0:istride:(isz-1)*istride
                    v = A[k+i]
                    write(stream, redval(v))
                    write(stream, greenval(v))
                    write(stream, blueval(v))
                end
            end
        else
            for i = 1:isz
                for j = 1:jsz
                    k = firstindex+(i-1)*istride+(j-1)*jstride
                    v = A[k]
                    write(stream, redval(v))
                    write(stream, greenval(v))
                    write(stream, blueval(v))
                end
            end
        end
    elseif cs == "ARGB32"
        if xfirst
            for j = 1:jsz
                k = firstindex + (j-1)*jstride
                for i = 0:istride:(isz-1)*istride
                    v = A[k+i]
                    write(stream, redval(v))
                    write(stream, greenval(v))
                    write(stream, blueval(v))
                    write(stream, alphaval(v))
                end
            end
        else
            for i = 1:isz
                for j = 1:jsz
                    k = firstindex+(i-1)*istride+(j-1)*jstride
                    v = A[k]
                    write(stream, redval(v))
                    write(stream, greenval(v))
                    write(stream, blueval(v))
                    write(stream, alphaval(v))
                end
            end
        end
    else
        error("colorspace ", cs, " not yet supported")
    end
end

## PPM, PGM, and PBM ##
type PPMBinary <: ImageFileType end
type PGMBinary <: ImageFileType end
type PBMBinary <: ImageFileType end

add_image_file_format(".ppm", b"P6", PPMBinary)
#add_image_file_format(".ppm", b"P3", PPMASCII)
add_image_file_format(".pgm", b"P5", PGMBinary)
#add_image_file_format(".pgm", b"P2", PGMASCII)
add_image_file_format(".pbm", b"P4", PBMBinary)
#add_image_file_format(".pbm", b"P1", PBMASCII)

function parse_netpbm_size(stream::IO)
    szline = strip(readline(stream))
    while isempty(szline) || szline[1] == '#'
        szline = strip(readline(stream))
    end
    parseints(szline, 2)
end

function parse_netpbm_maxval(stream::IO)
    skipchars(stream, isspace, linecomment='#')
    maxvalline = strip(readline(stream))
    parseint(maxvalline)
end

function imread{S<:IO}(stream::S, ::Type{PPMBinary})
    w, h = parse_netpbm_size(stream)
    maxval = parse_netpbm_maxval(stream)
    local dat
    if maxval <= 255
        dat = read(stream, Uint8, 3, w, h)
    elseif maxval <= typemax(Uint16)
        dat = Array(Uint16, 3, w, h)
        # there is no endian standard, but netpbm is big-endian
        if ENDIAN_BOM == 0x01020304
            for indx = 1:3*w*h
                dat[indx] = read(stream, Uint16)
            end
        else
            for indx = 1:3*w*h
                dat[indx] = bswap(read(stream, Uint16))
            end
        end
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
    T = eltype(dat)
    Image(dat, ["colorspace" => "RGB", "colordim" => 1, "spatialorder" => ["x", "y"], "limits" => (zero(T),convert(T,maxval))])
end

function imread{S<:IO}(stream::S, ::Type{PGMBinary})
    w, h = parse_netpbm_size(stream)
    maxval = parse_netpbm_maxval(stream)
    local dat
    if maxval <= 255
        dat = read(stream, Uint8, w, h)
    elseif maxval <= typemax(Uint16)
        dat = Array(Uint16, w, h)
        if ENDIAN_BOM == 0x01020304
            for indx = 1:w*h
                dat[indx] = read(stream, Uint16)
            end
        else
            for indx = 1:w*h
                dat[indx] = bswap(read(stream, Uint16))
            end
        end
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
    T = eltype(dat)
    Image(dat, ["colorspace" => "Gray", "spatialorder" => ["x", "y"], "limits" => (zero(T),convert(T,maxval))])
end

function imread{S<:IO}(stream::S, ::Type{PBMBinary})
    w, h = parse_netpbm_size(stream)
    dat = BitArray(w, h)
    nbytes_per_row = iceil(w/8)
    for irow = 1:h, j = 1:nbytes_per_row
        tmp = read(stream, Uint8)
        offset = (j-1)*8
        for k = 1:min(8, w-offset)
            dat[offset+k, irow] = (tmp>>>(8-k))&0x01
        end
    end
    Image(dat, ["spatialorder" => ["x", "y"]])
end

function imwrite(img, filename::String, ::Type{PPMBinary})
    open(filename, "w") do stream
        imwrite(img, stream, PPMBinary)
    end
end

function imwrite(img, s::IO, ::Type{PPMBinary})
    write(s, "P6\n")
    write(s, "# ppm file written by Julia\n")
    w, h = widthheight(img)
    bitdepth = 8*sizeof(eltype(img))
    if eltype(img) <: FloatingPoint
        bitdepth = 8
    end
    cs = colorspace(img)
    if cs == "RGB24"
        bitdepth = 8
    elseif cs == "ARGB32"
        bitdepth = 8
    end
    T = typedict[bitdepth]
    mx = int(typemax(T))
    scalei = scaleinfo(T, img)
    write(s, "$w $h\n$mx\n")
    writecolor(s, img, scalei)
end

# ## PNG ##
# type PNGFile <: ImageFileType end
#
# add_image_file_format(".png", [0x89,0x50,0x4E,0x47,0x0D,0x0A,0x1A,0x0A], PNGFile)
#
# libpng = dlopen("libpng")
# libimage_jl = dlopen(strcat(JULIA_HOME, "/../../extras/libjlimage"))
# png_major   = ccall(dlsym(libimage_jl, :jl_png_libpng_ver_major), Int32, ())
# png_minor   = ccall(dlsym(libimage_jl, :jl_png_libpng_ver_minor), Int32, ())
# png_release = ccall(dlsym(libimage_jl, :jl_png_libpng_ver_release), Int32, ())
# jl_png_read_init           = dlsym(libimage_jl, :jl_png_read_init)
# jl_png_write_init          = dlsym(libimage_jl, :jl_png_write_init)
# png_set_sig_bytes          = dlsym(libpng, :png_set_sig_bytes)
# png_read_info              = dlsym(libpng, :png_read_info)
# png_write_info             = dlsym(libpng, :png_write_info)
# png_read_update_info       = dlsym(libpng, :png_read_update_info)
# png_get_image_width        = dlsym(libpng, :png_get_image_width)
# png_get_image_height       = dlsym(libpng, :png_get_image_height)
# png_get_bit_depth          = dlsym(libpng, :png_get_bit_depth)
# png_get_color_type         = dlsym(libpng, :png_get_color_type)
# #png_get_filter_type        = dlsym(libpng, :png_get_filter_type)
# #png_get_compression_type   = dlsym(libpng, :png_get_compression_type)
# #png_get_interlace_type     = dlsym(libpng, :png_get_interlace_type)
# #png_set_interlace_handling = dlsym(libpng, :png_set_interlace_handling)
# png_set_expand_gray_1_2_4_to_8 = convert(Ptr{Void}, 0)
# try
#     png_set_expand_gray_1_2_4_to_8 = dlsym(libpng, :png_set_expand_gray_1_2_4_to_8)
# catch
#     png_set_expand_gray_1_2_4_to_8 = dlsym(libpng, :png_set_gray_1_2_4_to_8)
# end
# png_set_swap               = dlsym(libpng, :png_set_swap)
# png_set_IHDR               = dlsym(libpng, :png_set_IHDR)
# png_get_valid              = dlsym(libpng, :png_get_valid)
# png_get_rowbytes           = dlsym(libpng, :png_get_rowbytes)
# jl_png_read_image          = dlsym(libimage_jl, :jl_png_read_image)
# jl_png_write_image          = dlsym(libimage_jl, :jl_png_write_image)
# jl_png_read_close          = dlsym(libimage_jl, :jl_png_read_close)
# jl_png_write_close          = dlsym(libimage_jl, :jl_png_write_close)
# const PNG_COLOR_MASK_PALETTE = 1
# const PNG_COLOR_MASK_COLOR   = 2
# const PNG_COLOR_MASK_ALPHA   = 4
# const PNG_COLOR_TYPE_GRAY = 0
# const PNG_COLOR_TYPE_PALETTE = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_PALETTE
# const PNG_COLOR_TYPE_RGB = PNG_COLOR_MASK_COLOR
# const PNG_COLOR_TYPE_RGB_ALPHA = PNG_COLOR_MASK_COLOR | PNG_COLOR_MASK_ALPHA
# const PNG_COLOR_TYPE_GRAY_ALPHA = PNG_COLOR_MASK_ALPHA
# const PNG_INFO_tRNS = 0x0010
# const PNG_INTERLACE_NONE = 0
# const PNG_INTERLACE_ADAM7 = 1
# const PNG_COMPRESSION_TYPE_BASE = 0
# const PNG_FILTER_TYPE_BASE = 0
#
# png_color_type(::Type{CSGray}) = PNG_COLOR_TYPE_GRAY
# png_color_type(::Type{RGB}) = PNG_COLOR_TYPE_RGB
# png_color_type(::Type{CSRGBA}) = PNG_COLOR_TYPE_RGB_ALPHA
# png_color_type(::Type{CSGrayAlpha}) = PNG_COLOR_TYPE_GRAPH_ALPHA
#
# # Used by finalizer to close down resources
# png_read_cleanup(png_voidp) = ccall(jl_png_read_close, Void, (Ptr{Ptr{Void}},), png_voidp)
# png_write_cleanup(png_voidp) = ccall(jl_png_write_close, Void, (Ptr{Ptr{Void}},), png_voidp)
#
# function imread{S<:IO}(stream::S, ::Type{PNGFile})
#     png_voidp = ccall(jl_png_read_init, Ptr{Ptr{Void}}, (Ptr{Void},), fd(stream))
#     if png_voidp == C_NULL
#         error("Error opening PNG file ", stream.name, " for reading")
#     end
#     finalizer(png_voidp, png_read_cleanup)  # gracefully handle errors
#     png_p = pointer_to_array(png_voidp, (3,))
#     # png_p[1] is the png_structp, png_p[2] is the png_infop, and
#     # png_p[3] is a FILE* created from stream
#     # Tell the library how many header bytes we've already read
#     ccall(png_set_sig_bytes, Void, (Ptr{Void}, Int32), png_p[1], position(stream))
#     # Read the rest of the header
#     ccall(png_read_info, Void, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#     # Determine image parameters
#     width = ccall(png_get_image_width, Uint32, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#     height = ccall(png_get_image_height, Uint32, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#     bit_depth = ccall(png_get_bit_depth, Int32, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#     if bit_depth <= 8
#         T = Uint8
#     elseif bit_depth == 16
#         T = Uint16
#     else
#         error("Unknown pixel type")
#     end
#     color_type = ccall(png_get_color_type, Int32, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#     if color_type == PNG_COLOR_TYPE_GRAY
#         n_color_channels = 1
#         cs = CSGray
#     elseif color_type == PNG_COLOR_TYPE_RGB
#         n_color_channels = 3
#         cs = RGB
#     elseif color_type == PNG_COLOR_TYPE_GRAY_ALPHA
#         n_color_channels = 2
#         cs = CSGrayAlpha
#     elseif color_type == PNG_COLOR_TYPE_RGB_ALPHA
#         n_color_channels = 4
#         cs = CSRGBA
#     else
#         error("Color type not recognized")
#     end
#     # There are certain data types we just don't want to handle in
#     # their raw format. Check for these, and if needed ask libpng to
#     # transform them.
#     if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
#         ccall(png_set_expand_gray_1_2_4_to_8, Void, (Ptr{Void},), png_p[1])
#         bit_depth = 8
#     end
#     if ccall(png_get_valid, Uint32, (Ptr{Void}, Ptr{Void}, Uint32), png_p[1], png_p[2], PNG_INFO_tRNS) > 0
#         call(png_set_tRNS_to_alpha, Void, (Ptr{Void},), png_p[1])
#     end
#     # TODO: paletted images
#     if bit_depth > 8 && ENDIAN_BOM == 0x04030201
#         call(png_set_swap, Void, (Ptr{Void},), png_p[1])
#     end
#     # TODO? set a background
#     ccall(png_read_update_info, Void, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#     # Allocate the buffer
#     pixbytes = sizeof(T) * n_color_channels
#     rowbytes = pixbytes*width
#     if rowbytes != ccall(png_get_rowbytes, Uint32, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#         error("Number of bytes per row is wrong")
#     end
#     if n_color_channels == 1
#         data = Array(T, width, height)
#         order = "xy"
#     else
#         data = Array(T, n_color_channels, width, height)
#         order = "cxy"
#     end
#     # Read the image
#     status = ccall(jl_png_read_image, Int32, (Ptr{Void}, Ptr{Void}, Ptr{Uint8}), png_p[1], png_p[2], data)
#     if status < 0
#         error("PNG error: read failed")
#     end
#     # No need to clean up, the finalizer will do it
#     # Build the image
#     return ImageCS{T,2,cs}(data,order)
# end
#
# function imwrite{S<:IO}(img, stream::S, ::Type{PNGFile}, opts::Options)
#     dat = uint(img)
#     cs = colorspace(img)
#     # Ensure PNG storage order
#     if cs != CSGray
#         dat = permute_c_first(dat, storageorder(img))
#     end
#     imwrite(dat, cs, stream, PNGFile, opts)
# end
# imwrite{S<:IO}(img, stream::S, ::Type{PNGFile}) = imwrite(img, stream, PNGFile, Options())
#
# function imwrite(img, filename::ByteString, ::Type{PNGFile}, opts::Options)
#     s = open(filename, "w")
#     imwrite(img, s, PNGFile, opts)
#     gc()  # force the finalizer to run
# end
# imwrite(img, filename::ByteString, ::Type{PNGFile}) = imwrite(img, filename, PNGFile, Options())
#
# function imwrite{T<:Unsigned, CS<:ColorSpace, S<:IO}(data::Array{T}, ::Type{CS}, stream::S, ::Type{PNGFile}, opts::Options)
#     @defaults opts interlace=PNG_INTERLACE_NONE compression=PNG_COMPRESSION_TYPE_BASE filter=PNG_FILTER_TYPE_BASE
#     png_voidp = ccall(jl_png_write_init, Ptr{Ptr{Void}}, (Ptr{Void},), fd(stream))
#     if png_voidp == C_NULL
#         error("Error opening PNG file ", stream.name, " for writing")
#     end
#     finalizer(png_voidp, png_write_cleanup)
#     png_p = pointer_to_array(png_voidp, (3,))
#     # Specify the parameters and write the header
#     if CS == CSGray
#         height = size(data, 2)
#         width = size(data, 1)
#     else
#         height = size(data, 3)
#         width = size(data, 2)
#     end
#     bit_depth = 8*sizeof(T)
#     ccall(png_set_IHDR, Void, (Ptr{Void}, Ptr{Void}, Uint32, Uint32, Int32, Int32, Int32, Int32, Int32), png_p[1], png_p[2], width, height, bit_depth, png_color_type(CS), interlace, compression, filter)
#     # TODO: support comments via a comments=Array(PNGComments, 0) option
#     ccall(png_write_info, Void, (Ptr{Void}, Ptr{Void}), png_p[1], png_p[2])
#     if bit_depth > 8 && ENDIAN_BOM == 0x04030201
#         call(png_set_swap, Void, (Ptr{Void},), png_p[1])
#     end
#     # Write the image
#     status = ccall(jl_png_write_image, Int32, (Ptr{Void}, Ptr{Void}, Ptr{Uint8}), png_p[1], png_p[2], data)
#     if status < 0
#         error("PNG error: read failed")
#     end
#     # Clean up now, to prevent problems that otherwise occur when we
#     # try to immediately use the file from an external program before
#     # png_voidp gets garbage-collected
# #    finalizer_del(png_voidp)  # delete the finalizer
# #    png_write_cleanup(png_voidp)
# end
# imwrite(data, cs, stream, ::Type{PNGFile}) = imwrite(data, cs, stream, PNGFile, Options())

function parseints(line, n)
    ret = Array(Int, n)
    pos = 1
    for i = 1:n
        pos2 = search(line, ' ', pos)
        if pos2 == 0
            pos2 = length(line)+1
        end
        ret[i] = parseint(line[pos:pos2-1])
        pos = pos2+1
        if pos > length(line) && i < n
            error("Line terminated without finding all ", n, " integers")
        end
    end
    tuple(ret...)

end

# Permute to a color, horizontal, vertical, ... storage order (with time always last)
function permutedims_canonical(img)
    cd = colordim(img)
    td = timedim(img)
    p = spatialpermutation(["x", "y"], img)
    if cd != 0
        insert!(p, 1, cd)
    end
    if td != 0
        push!(p, td)
    end
    permutedims(img, p)
end

### Register formats for later loading here
type Dummy <: ImageFileType; end
add_image_file_format(".dummy", b"Dummy Image", Dummy, "dummy.jl")

# NRRD image format
type NRRDFile <: ImageFileType end
add_image_file_format(".nrrd", b"NRRD", NRRDFile, "nrrd.jl")
add_image_file_format(".nhdr", b"NRRD", NRRDFile, "nrrd.jl")

# Andor Technologies SIF file format
type AndorSIF <: Images.ImageFileType end
add_image_file_format(".sif", b"Andor Technology Multi-Channel File", AndorSIF, "SIF.jl")

# Imagine file format (http://holylab.wustl.edu, "Software" tab)
type ImagineFile <: ImageFileType end
add_image_file_format(".imagine", b"IMAGINE", ImagineFile, "Imagine.jl")

# PCO b16 image format
type B16File <: ImageFileType end
add_image_file_format(".b16", b"PCO-", B16File, "b16.jl")


text:
mutable struct piecewise_polynomial_t  <: LCMType
    timestamp::Int64
    num_breaks::Int32
    breaks::Vector{Float64}
    num_segments::Int32
    polynomial_matrices::Vector{polynomial_matrix_t}
end

@lcmtypesetup(piecewise_polynomial_t,
    breaks => (num_breaks,),
    polynomial_matrices => (num_segments,)
)


text:
using LightGraphs
include("../src/utils.jl")

"""
    lattice_nearest_neighbor(Lx::Int,Ly::Int)

lattice with nearest neighbor interactions, where nearest neighbor here
also includes diagonal neighbors.
"""
function lattice_nearest_neighbor(Lx::Int,Ly::Int)
    N=Lx*Ly
    lattice=SimpleGraph(N)
    for y in 1:Ly
        for x in 1:Lx
            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x-1,y+1,Lx,Ly))
            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x  ,y+1,Lx,Ly))
            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x+1,y+1,Lx,Ly))
            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x-1,y  ,Lx,Ly))

            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x+1,y  ,Lx,Ly))
            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x-1,y-1,Lx,Ly))
            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x  ,y-1,Lx,Ly))
            add_edge!(lattice, pbc_index(x,y,Lx,Ly), pbc_index(x+1,y-1,Lx,Ly))
        end
    end
    return lattice
end

"""
    step!(rng, state, state_new, lattice, p_rec[, state_rev, val_ref])

recurrent update step with stochastic activations of `state_new` from `state`
based on interactions defined in `lattice` with probability `p_rec` using
random number generator `rng`.

External activations can be added to `state_new` before calling `step!`.

Optionally, one can pass the refractory state `state_rev` that is reinitalized
upon activation to `val_ref`, decreases one value per time step, and makes
neurons only excitable if `state_ref[i]>0`.
"""
function step!(
        rng::AbstractRNG,
        state::Array{Bool},
        state_new::Array{Bool},
        lattice::SimpleGraph{Int64},
        p_rec
    )
    prepare!(rng, state, state_new, lattice, p_rec)

    # swap
    state .= state_new
    state_new .= false
    return nothing
end
function step!(
        rng::AbstractRNG,
        state::Array{Bool},
        state_new::Array{Bool},
        lattice::SimpleGraph{Int64},
        p_rec,
        state_ref::Array{Int},
        val_ref::Int,
    )
    prepare!(rng, state, state_new, lattice, p_rec)

    # if still in refractory period then no update
    for i in 1:length(state)
        if state_ref[i] > 0
            state_new[i] = false
        end
    end

    # swap
    state .= state_new
    state_new .= false

    # update refractory period
    for i in 1:length(state)
        if state[i]
            state_ref[i] = val_ref
        else
            state_ref[i] -= 1
        end
    end
    return nothing
end
function prepare!(rng, state,state_new, lattice, p_rec)
    for i in 1:length(state)
        if state[i]==true
            for j in outneighbors(lattice, i)
                if rand(rng) < p_rec
                    state_new[j] = true
                end
            end
        end
    end
end


text:
using Unitful: ustrip, dimension
using Crystals.Utilities: eldimension, CellParameters

@testset "> Hart-Forcade" begin
    lattice = Crystal([0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm", position=[0, 0, 0]u"nm")
    crystal = Crystal([1 0 0; 0 1 0; 0 0 1]u"nm",
                      position=[0, 0, 0]u"nm",
                      position=[0.5, 0.5, 0]u"nm",
                      position=[0.5, 0, 0.5]u"nm",
                      position=[0, 0.5, 0.5]u"nm")
    hf = @inferred(hart_forcade(lattice.cell, crystal.cell))
    @test hf.quotient == [1, 2, 2]
    for i in eachindex(crystal)
        @test all(isinteger, hf.transform * crystal[i, :position])
    end
    for i=1:hf.quotient[1], j=1:hf.quotient[2], k=1:hf.quotient[3]
        @test all(isinteger, inv(lattice.cell) * (inv(hf.transform) * [i, j, k]))
    end
end

@testset "> eldimension" begin
    @test @inferred(eldimension(1u"m")) === typeof(dimension(u"m"))
    @test @inferred(eldimension(1)) === typeof(dimension(1))
    @test @inferred(eldimension([1, 1]u"m")) === typeof(dimension(u"m"))
    @test @inferred(eldimension([1, 1])) === typeof(dimension(1))
end

# @testset "> is_position_array" begin
#     @test !is_position_array([1, 1])
#     @test !is_position_array(Position(1, 1))
#     @test is_position_array([Position(1, 1)])
#     @test is_position_array([1 1; 1 1])
# end
#
@testset "> to fractional" begin
  cell = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm"
  pos = [0.1, 0.1, 0.2]u"nm"
  @test @inferred(to_fractional(pos, cell)) ≈ inv(cell) * pos
  uposition = ustrip.(pos)
  @test @inferred(to_fractional(uposition, cell)) === uposition
end

@testset "> to cartesian" begin
    cell = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm"
    pos = [0.1, 0.1, 0.2]u"nm"
    @test @inferred(to_cartesian(pos, cell)) === pos
    uposition = ustrip.(pos)
    @test @inferred(to_cartesian(uposition, cell)) ≈ cell * uposition
end

@testset "> Periodic images" begin
  cell = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm"

  for i in 1:1 #10
    pos = cell * rand(3)
    image = pos + cell * rand(Int8, (3, ))
    @test @inferred(is_periodic(pos, image, cell))
    @test !@inferred(is_periodic(pos, image + (rand(3) * 1e-4) * u"nm", cell))

    positions = cell * rand((3, 10))
    images = positions + cell * rand(Int8, (3, 10))
    positions[:, [1, 3, 4]] += rand((3, 3))u"nm"
    @test @inferred(is_periodic(positions, images, cell))[2]
    @test all(is_periodic(positions, images, cell)[5:end])
    @test !any(is_periodic(positions, images, cell)[[1, 3, 4]])
  end
end

@testset "> Fold back into cell" begin
  cell = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm"

  for i in 1:10
    pos = cell * (rand(3) + rand(Int8, (3, )))
    folded = into_cell(pos, cell)
    @test @inferred(is_periodic(pos, folded, cell))
    @test all(0 .≤ inv(cell) * folded .< 1)

    pos = cell * (rand((3, 10)) + rand(Int8, (3, 10)))
    folded = into_cell(pos, cell)
    @test all(@inferred(is_periodic(pos, folded, cell)))
    @test all(0 .≤ inv(cell) * folded .< 1)
  end
end

@testset "> Fold back around origin" begin
  cell = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm"

  for i in 1:10
    pos = cell * (rand(3) + rand(Int8, (3, )))
    folded = @inferred(origin_centered(pos, cell))
    @test is_periodic(pos, folded, cell)
    @test all(-0.5 .< inv(cell) * folded .≤ 0.5)

    pos = cell * (rand((3, 10)) + rand(Int8, (3, 10)))
    folded = @inferred(origin_centered(pos, cell))
    @test all(is_periodic(pos, folded, cell))
    @test all(-0.5 .< inv(cell) * folded .≤ 0.5)
  end
end

function none_smaller(pos::Vector, cell::Matrix)
  const d = norm(pos)
  for i = -2:2, j = -2:2, k = -2:2
    i == j == k == 0 && continue
    norm(pos + cell * [i, j, k]) < d * (1 - 1e-12)  && return false
  end
  true
end

@testset "> Fold back into voronoi" begin
  cell = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm"

  for i in 1:10
    pos = cell * (rand(3) + rand(Int8, (3, )))
    folded = into_voronoi(pos, cell)
    @test is_periodic(pos, folded, cell)
    @test none_smaller(folded, cell)

    pos = cell * (rand((3, 10)) + rand(Int8, (3, 10)))
    folded = into_voronoi(pos, cell)
    @test all(is_periodic(pos, folded, cell))
    for i = 1:size(pos, 2)
      @test none_smaller(folded[:, i], cell)
    end
  end
end

@testset "> Supercell" begin
  lattice = Crystal(
    [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm",
    position=[0 0.25; 0 0.25; 0 0.25]u"nm", species=["In", "Ga"]
  )

  result = supercell(lattice, lattice.cell * [-1 1 1; 1 -1 1; 1 1 -1])
  @test result.cell ≈ eye(3)u"nm"
  @test nrow(result) % nrow(lattice) == 0
  const ncells =  nrow(result) // nrow(lattice)
  @test ncells == 4
  @test names(result) == ∪(names(result), (:site_id, ))
  @test all(1 .≤ result[:site_id] .≤ nrow(lattice))
  actual = convert(Array, result[:position])
  expected = convert(Array, lattice[result[:site_id], :position])
  @test all(is_periodic(actual, expected, lattice.cell))
  @test length(unique(result[:site_id])) == nrow(lattice)
  for i in unique(result[:site_id])
    @test countnz(result[:site_id] .== i) == ncells
  end
  for atom in eachatom(result)
      @test countnz(is_periodic(atom[:cartesian], lattice[:cartesian], lattice.cell)) == 1
      @test countnz(is_periodic(atom[:cartesian], result[:cartesian], result.cell)) == 1
  end
end

@testset "> cell parameters" begin
    cells = Matrix[
        [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]u"nm",
    ]
    parameters = CellParameters[
      CellParameters(sqrt(0.5)u"nm", sqrt(0.5)u"nm", sqrt(0.5)u"nm", 60u"°", 60u"°", 60u"°")
    ]
    for (cell, params) in zip(cells, parameters)
        actual_params = cell_parameters(cell)
        @test actual_params.a ≈ params.a
        @test actual_params.b ≈ params.b
        @test actual_params.c ≈ params.c
        @test actual_params.α ≈ params.α
        @test actual_params.β ≈ params.β
        @test actual_params.γ ≈ params.γ

        round_trip = cell_parameters(cell_parameters(actual_params))
        @test round_trip.a ≈ params.a
        @test round_trip.b ≈ params.b
        @test round_trip.c ≈ params.c
        @test round_trip.α ≈ params.α
        @test round_trip.β ≈ params.β
        @test round_trip.γ ≈ params.γ
        @test det(inv(cell) * cell_parameters(actual_params)) ≈ 1
    end
end


text:
using CompressionStreams

function dotest(fn, N)

    buf = Vector{UInt8}(2^20)

    for i in 1:N
        s = cstream(fn)
        red = readbytes!(s, buf)
    end

    assert(all(buf .== 0))

end

#@time dotest("1M.bin", 10)
#@time dotest("1M.bin.gz", 10)

@time dotest("1G.bin", 2)
@time dotest("1G.bin.gz", 2)


text:
using Random

function generate_erdosrenyi(n)

	dismat 			=	zeros(n,n)
	permlength 		=	(n^2 - n)/2
	permlength 		=	round(Int64, permlength)
	perm 			=	randperm(permlength)

	counter 		=	0
	for i = 1:n
		for j = (i+1):n
			counter =	counter+1
			dismat[i,j] 	=	perm[counter]
			dismat[j,i] 	=	perm[counter]
		end
	end

	return dismat

end


text:
struct VolumeCoarsening <: AbstractAlgebraicCoarsening
    Q::Float64
    η::Float64
    order::Int
end

# TODO maybe reuse order vec from choose_seeds to save space
function split(vc::VolumeCoarsening, W, volume)
    isseed, vol_order = _choose_seeds(vc, W, volume)
    nc = count(isseed)
    nf = length(isseed) - nc

    C = zeros(Int, nc)
    F = zeros(Int, nf)
    invseed = vol_order

    idxc = 0
    idxf = 0
    for (i, t) in enumerate(isseed)
        if t
            idxc += 1
            C[idxc] = i
            invseed[i] = idxc
        else
            idxf += 1
            F[idxf] = i
            invseed[i] = idxf + nc
        end
    end

    # invseed equivalent to invperm([C; F])
    return C, F, invseed
end

# TODO verify this is working correctly
function _choose_seeds(vc::VolumeCoarsening, W, volume)
    n = size(W, 1)
    isseed = falses(n)
    coarse_connection = zeros(n)
    total_strength = sumrows(W)
    fvolume = futurevolume(W, volume, total_strength)
    mean_volume = sum(fvolume) / length(fvolume)
    vol_order = sortperm(fvolume; rev=true)

    for col in vol_order
        isabovemean = fvolume[col] > mean_volume * vc.η
        wellconnected = coarse_connection[col] / total_strength[col] > vc.Q
        if isabovemean || !wellconnected
            isseed[col] = true
            addcol!(coarse_connection, W, col)
        end
    end
    return isseed, vol_order
end

""
function futurevolume(W, volume, total_strength)
    future = copy(volume)
    rows = rowvals(W)
    vals = nonzeros(W)
    for j in axes(W, 2)
        for i in nzrange(W, j)
            row = rows[i]
            val = vals[i]
            future[j] += volume[row] * val / total_strength[row]
        end
    end
    return future
end

function coarseop(vc::VolumeCoarsening, A; volume, strength, _...)
    C, F, invseed = split(vc, strength, volume)
    if length(F) > 0
        N = coarse_neighborhoods(strength, C, invseed, order(vc))
        P = formcoarseop(vc, C, F, N)
    else
        n = length(C)
        P = sparse(I, n, n)
    end
    return P, C, F
end

function coarsen(vc::VolumeCoarsening, A; volume, strength, _...)
    P, C, F = coarseop(vc, A; volume=volume, strength=strength)
    Ac = P' * (A * P)
    return Ac, P, C, F
end


text:
module BouncingSphere3D

using Modia

BouncingSphere = Model(
    boxHeigth = 0.1,
    world      = Object3D(feature=Scene(enableContactDetection = true, # default value
    animationFile="BouncingSphere.json")),
    ground     = Object3D(parent=:world, translation=:[0.0,-boxHeigth/2,0.0],
                    feature=Solid(shape=Box(lengthX=4.0, lengthY=:boxHeigth, lengthZ=0.7),
                        visualMaterial=VisualMaterial(color="DarkGreen"),
                        solidMaterial="Steel", # for mass and force computation
                        collision=true)),      # enable collision flag
    sphere     = Object3D(feature=Solid(shape=Sphere(diameter=0.2),
                        visualMaterial=VisualMaterial(color="Blue"),
                        solidMaterial="Steel", # for mass and force computation
                        massProperties=MassPropertiesFromShapeAndMass(mass=0.001),
                        collision=true)),      # enable collision flag
    free       = FreeMotion(obj1=:world, obj2=:sphere, r=Var(init=[0.0, 1.0, 0.0]))
)

bouncingSphere = @instantiateModel(buildModia3D(BouncingSphere), unitless=true)
simulate!(bouncingSphere, stopTime=2.2, dtmax=0.1)

@usingModiaPlot
plot(bouncingSphere, "free.r", figure=1)

end


text:
"""
    amr = loadfortq(filename::AbstractString, ncol::Integer; vartype::Symbol=:surface,
                    params::VisClaw.GeoParam=VisClaw.GeoParam(), runup::Bool=true,
                    xlims=(-Inf,Inf), ylims=(-Inf,Inf), region="", AMRlevel=[])

Function: fort.qxxxx reader.
"""
function loadfortq(filename::AbstractString, ncol::Integer; vartype::Symbol=:surface,
                   params::VisClaw.GeoParam=VisClaw.GeoParam(), runup::Bool=true,
                   xlims=(-Inf,Inf), ylims=(-Inf,Inf), region="", AMRlevel=[])
    ## check
    !any(map(sym -> vartype == sym, [:surface, :current, :storm])) && error("kwarg 'vartype' is invalid")

    ## set range
    if isa(region, VisClaw.AbstractTopo); xlims=extrema(region.x); ylims=extrema(region.y); end
    if isa(region, VisClaw.Region); xlims=region.xlims; ylims=region.ylims; end

    ## file open
    f = open(filename,"r")
    txtorg = readlines(f)
    close(f) #close

    ## count the number of lines and grids
    nlineall = length(txtorg)
    idx = occursin.("grid_number",txtorg)
    ngrid = length(txtorg[idx])

    if vartype==:surface
        amr = Array{VisClaw.SurfaceHeight}(undef,ngrid) ## preallocate
    elseif vartype==:current
        amr = Array{VisClaw.Velocity}(undef,ngrid)
    elseif vartype==:storm
        amr = Array{VisClaw.Storm}(undef,ngrid)
    end

    l = 1
    i = 1
    while l < nlineall
        ## read header
        #header = txtorg[1:8]
        header = txtorg[l:l+7]
        header = map(strip,header)
        gridnumber = parse(Int64, split(header[1],r"\s+")[1])
        AMRlevel_load = parse(Int64, split(header[2],r"\s+")[1])
        mx = parse(Int64, split(header[3],r"\s+")[1])
        my = parse(Int64, split(header[4],r"\s+")[1])
        xlow = parse(Float64, split(header[5],r"\s+")[1])
        ylow = parse(Float64, split(header[6],r"\s+")[1])
        dx = parse(Float64, split(header[7],r"\s+")[1])
        dy = parse(Float64, split(header[8],r"\s+")[1])
        ## read variables
        body = txtorg[l+9:l+9+(mx+1)*my-1]

        # the next tile
        l = l+9+(mx+1)*my
        ## check AMRlevel
        if !isempty(AMRlevel); if isempty(findall(AMRlevel .== AMRlevel_load)); i += 1; continue; end; end

        ## check whether the tile is on the domain
        if (xlow+dx*mx < xlims[1]) | (xlims[2] < xlow); i += 1; continue; end
        if (ylow+dy*my < ylims[1]) | (ylims[2] < ylow); i += 1; continue; end

        if vartype==:surface
            elev = [parse(Float64, body[(i-1)*(mx+1)+j][26*(ncol-1)+1:26*ncol]) for i=1:my, j=1:mx]
            depth = [parse(Float64, body[(i-1)*(mx+1)+j][1:26]) for i=1:my, j=1:mx]
            # wet condition
            land = (elev-depth) .>= params.dmin

            # sea surface anomaly
            (params.eta0 != 0.0) && (elev[.!land] = elev[.!land].-params.eta0)

            # inundation depth if wet
            runup && (elev[land] = depth[land])

            # NaN if dry
            elev[depth.<=0.0] .= NaN

            ## array
            amr[i] = VisClaw.SurfaceHeight(gridnumber,AMRlevel_load,mx,my,xlow,ylow,dx,dy,elev)

        elseif vartype==:current
            ucol = ncol
            vcol = ncol+1
            # read
            depth = [parse(Float64, body[(i-1)*(mx+1)+j][1:26]) for i=1:my, j=1:mx]
            u = [parse(Float64, body[(i-1)*(mx+1)+j][26*(ucol-1)+1:26*ucol]) for i=1:my, j=1:mx]
            v = [parse(Float64, body[(i-1)*(mx+1)+j][26*(vcol-1)+1:26*vcol]) for i=1:my, j=1:mx]
            # replace to NaN
            mask = depth.<=0.0
            depth[mask] .= NaN
            u[mask] .= NaN
            v[mask] .= NaN
            # calc
            u = u./depth
            v = v./depth
            vel = sqrt.(u.^2 .+ v.^2)
            ## array
            amr[i] = VisClaw.Velocity(gridnumber,AMRlevel_load,mx,my,xlow,ylow,dx,dy,u,v,vel)

        elseif vartype==:storm
            ucol = ncol
            vcol = ncol+1
            pcol = ncol+2
            u = [parse(Float64, body[(i-1)*(mx+1)+j][26*(ucol-1)+1:26*ucol]) for i=1:my, j=1:mx]
            v = [parse(Float64, body[(i-1)*(mx+1)+j][26*(vcol-1)+1:26*vcol]) for i=1:my, j=1:mx]
            p = [parse(Float64, body[(i-1)*(mx+1)+j][26*(pcol-1)+1:26*pcol]) for i=1:my, j=1:mx]
            p = p./1e+2

            # u[(abs.(u).<=1e-2) .& (abs.(v).<=1e-2)] .= NaN
            # v[(abs.(u).<=1e-2) .& (abs.(v).<=1e-2)] .= NaN

            ## array
            amr[i] = VisClaw.Storm(gridnumber,AMRlevel_load,mx,my,xlow,ylow,dx,dy,u,v,p)
        end
        ## print
        #@printf("%d, ",gridnumber)

        ## counter; go to the next grid
        i += 1
    end
    amr = amr[filter(i -> isassigned(amr, i), 1:length(amr))]
    ## return
    return amr
end
#################################
"""
    amr = loadforta(filename::AbstractString, ncol::Integer; kwargs...)

Function: fort.axxxx reader. See also [`loadfortq`](@ref).
"""
loadforta(filename::AbstractString, ncol::Integer; kwargs...) = loadfortq(filename, ncol; vartype=:storm, kwargs...)
#################################

#################################
"""
    timelaps = loadfortt(filename::AbstractString)

Function: fort.txxxx reader.
"""
function loadfortt(filename::AbstractString)
    ## file open
    f = open(filename,"r")
    txtorg = readlines(f)
    close(f) #close
    ## parse timelaps from the 1st line
    timelaps = parse(Float64, txtorg[1][1:18])
    ## return
    return timelaps
end
#################################

#######################################
"""
    amrs = loadsurface(loaddir::AbstractString, filesequence::AbstractVector; kwargs...)
    amrs = loadsurface(loaddir::AbstractString, fileid::Integer; kwargs...)

Function: load time-series of water surface.
          The keyword arguments follow [`loadfortq`](@ref).
          See also: [`loadfortt`](@ref).
"""
function loadsurface(loaddir::AbstractString, filesequence::AbstractVector=0:0; vartype::Symbol=:surface, kwargs...)

    # check
    !any(map(sym -> vartype == sym, [:surface, :current, :storm])) && error("kwarg 'vartype' is invalid")

    ## define the filepath & filename
    if vartype==:surface
        fnamekw = r"^fort\.q\d+$"
        col=4
    elseif vartype==:current
        fnamekw = r"^fort\.q\d+$"
        col=2
    elseif vartype==:storm
        fnamekw = r"^fort\.a\d+$"
        col=5
    end

    ## make a list
    isdir(loaddir) || error("Directory $loaddir doesn't exist")
    flist = readdir(loaddir)
    filter!(x->occursin(fnamekw, x), flist)
    isempty(flist) && error("File named $fnamekw was not found")

    # load geoclaw.data
    params = VisClaw.geodata(loaddir)

    ## the number of files
    nfile = length(flist)

    # file sequence to be loaded
    if filesequence==0:0; filesequence = 1:nfile; end
    (any(filesequence .< 1) || any(filesequence .> nfile)) && error("Incorrect file sequence was specified. (This must be from 1 to $nfile)")

    ## the number of files (to be loaded)
    nfile = length(filesequence)

    ## preallocate
    if vartype==:surface
        amr = Vector{AbstractVector{VisClaw.SurfaceHeight}}(undef,nfile)
    elseif vartype==:current
        amr = Vector{AbstractVector{VisClaw.Velocity}}(undef,nfile)
    elseif vartype==:storm
        amr = Vector{AbstractVector{VisClaw.Storm}}(undef,nfile)
    end

    ## load all files
    tlap = vec(zeros(nfile,1))
    cnt = 0
    for it = filesequence
        cnt += 1
        if vartype==:surface
            amr[cnt] = VisClaw.loadfortq(joinpath(loaddir,flist[it]), col; vartype=vartype, kwargs...)
        elseif vartype==:current
            amr[cnt] = VisClaw.loadfortq(joinpath(loaddir,flist[it]), col; vartype=vartype, kwargs...)
        elseif vartype==:storm
            amr[cnt] = VisClaw.loadforta(joinpath(loaddir,flist[it]), col; kwargs...)
        end
        tlap[cnt] = VisClaw.loadfortt(joinpath(loaddir,replace(flist[it],r"\.." => ".t")))
    end

    ## AMR Array
    amrs = VisClaw.AMR(nfile,tlap,amr)

    ## return value
    return amrs
end
#######################################
loadsurface(loaddir::AbstractString, fileid::Integer; kwargs...) =
loadsurface(loaddir, fileid:fileid; kwargs...)
#######################################

######################################
"""
    amrs = loadstorm(loaddir::AbstractString, filesequence::AbstractVector=0:0; kwargs...)
    amrs = loadstorm(loaddir::AbstractString, fileid::Integer; kwargs...)

Function: load time-series of storm data.
          The keyword arguments follow [`loadfortq`](@ref).
          See also: [`loadfortt`](@ref).
"""
loadstorm(loaddir::AbstractString, filesequence::AbstractVector=0:0; kwargs...) =
loadsurface(loaddir, filesequence; vartype=:storm, kwargs...)
#######################################
loadstorm(loaddir::AbstractString, fileid::Integer; kwargs...) =
loadsurface(loaddir, fileid:fileid; vartype=:storm, kwargs...)
#######################################

#######################################
"""
    amrs = loadcurrent(loaddir::AbstractString, filesequence::AbstractVector=0:0; kwargs...)
    amrs = loadcurrent(loaddir::AbstractString, fileid::Integer; kwargs...)

Function: load time-series of ocean current data.
          The keyword arguments follow [`loadfortq`](@ref).
          See also: [`loadfortt`](@ref).
"""
loadcurrent(loaddir::AbstractString, filesequence::AbstractVector=0:0; kwargs...) =
loadsurface(loaddir, filesequence, vartype=:current, kwargs...)
#######################################
loadcurrent(loaddir::AbstractString, fileid::Integer; kwargs...) =
loadsurface(loaddir, fileid:fileid; vartype=:current, kwargs...)
#######################################