text:
"""
    ComtradeData

`ComtradeData` holds data and configuration information from a set of COMTRADE
files. Fields include:

* `filename::UTF8String` -- Base file name (without extension).
* `cfg::ComtradeCfg` -- Configuration information.
* `dat::DataFrame` -- Table of results with each channel as a column.
* `hdr::UTF8String` -- Contents of the "hdr" file.
* `inf::UTF8String` -- Contents of the "inf" file.
"""
type ComtradeData
    filename::UTF8String
    cfg::ComtradeCfg
    dat::DataFrame
    hdr::UTF8String
    inf::UTF8String
end

function addcol!(d::DataFrame, x)
    push!(d.columns, x)
    push!(DataFrames.index(d), symbol(string("x", ncol(d) + 1)))
    d
end

"""
    read_comtrade(basename)

Read contents of a set of COMTRADE files with base file name `basename` (no
extension).

Returns a `COMTRADE.ComtradeData` with the following fields:

* `filename::UTF8String` -- Base file name (without extension).
* `cfg::ComtradeCfg` -- Configuration information.
* `dat::DataFrame` -- Table of results with each channel as a column.
* `hdr::UTF8String` -- Contents of the "hdr" file.
* `inf::UTF8String` -- Contents of the "inf" file.

See also `COMTRADE.ComtradeCfg` for a description of the `cfg` field.
"""
function read_comtrade(basename)
    cfg = read_cfg("$basename.cfg")
    if cfg.ft == "ASCII"
        d = readtable("$basename.dat", separator = ',', header = false)
        for i in 1:cfg.nA+2   # convert DataArray columns to Arrays
            d.columns[i] = d[i].data
        end
        for i in (1:cfg.nD) + 2 + cfg.nA
            d.columns[i] = map(Bool, d[i].data)
        end
    elseif cfg.ft == "BINARY"
        x = reinterpret(Int16, readbytes(open("$basename.dat")))
        n = cfg.nA + div(cfg.nD, 16) + 2 + 2
        x = reshape(x[1:(n * cfg.npts)], (n, cfg.npts))
        global d = DataFrame()
        addcol!(d, reinterpret(UInt32, x[1:2, :], (cfg.npts,)))
        addcol!(d, reinterpret(UInt32, x[3:4, :], (cfg.npts,)))
        for i = 1:cfg.nA
            addcol!(d, reshape(x[i+4, :], (cfg.npts,)))
        end
        for i = 1:cfg.nD
            addcol!(d, Bool[(x >>> (mod1(i, 16) - 1)) & 0x01 for x in x[cfg.nA + 4 + cld(i, 16), :]])
            # addcol!(d, Bool[(x >>> (16 - mod1(i, 16))) & 0x01 for x in x[cfg.nA + 4 + cld(i, 16), :]])
        end
    else
        error("Unsupported data type $(x.ft)")
    end
    newnames = DataFrames.make_unique(vcat(:n, :time, Symbol[symbol(x) for x in vcat(cfg.A[:ch_id], cfg.D[:ch_id])]))
    names!(d, newnames)
    # analog multipliers
    for i in 1:cfg.nA
        d.columns[i+2] = cfg.A[i,:a] * d[i+2] + cfg.A[i,:b]
    end
    # analog multipliers
    d.columns[2] = 1e-6 * cfg.timemult * d[2]
    hdr = isreadable("$basename.hdr") ? readall("$basename.hdr") : ""
    inf = isreadable("$basename.inf") ? readall("$basename.inf") : ""
    return ComtradeData(basename, cfg, d, hdr, inf)
end

"""
    S"some name"

Returns a symbol. Equivalent to `symbol("some name")`. This helps when COMTRADE
data columns use nonstandard names.
"""
macro S_str(arg)
    Meta.quot(symbol(arg))
end


text:
using BasicPOMCP
using Test
using POMDPs
using POMDPModels
using NBInclude
using D3Trees
using Random
using POMDPSimulators
using POMDPModelTools
using POMDPTesting
using POMDPLinter: @requirements_info, @show_requirements, requirements_info
using ParticleFilters: n_particles, particles, particle, weights, weighted_particles, weight_sum, weight

import POMDPs:
	transition,
	observation,
    reward,
    discount,
	initialstate,
	updater,
	states,
	actions,
	observations

struct ConstObsPOMDP <: POMDP{Bool, Symbol, Bool} end
updater(problem::ConstObsPOMDP) = DiscreteUpdater(problem)
initialstate(::ConstObsPOMDP) = BoolDistribution(0.0)
transition(p::ConstObsPOMDP, s::Bool, a::Symbol) = BoolDistribution(0.0)
observation(p::ConstObsPOMDP, a::Symbol, sp::Bool) = BoolDistribution(1.0)
reward(p::ConstObsPOMDP, s::Bool, a::Symbol, sp::Bool) = 1.
discount(p::ConstObsPOMDP) = 0.9
states(p::ConstObsPOMDP) = (true, false)
actions(p::ConstObsPOMDP) = (:the_only_action,)
observations(p::ConstObsPOMDP) = (true, false)

@testset "POMDPTesting" begin
	pomdp = BabyPOMDP()
	test_solver(POMCPSolver(), BabyPOMDP())
end;

@testset "type stability" begin
	pomdp = BabyPOMDP()
	solver = POMCPSolver(rng = MersenneTwister(1))
	planner = solve(solver, pomdp)
	b = initialstate(pomdp)
    tree = BasicPOMCP.POMCPTree(pomdp, b, solver.tree_queries)
    node = BasicPOMCP.POMCPObsNode(tree, 1)

    r = @inferred BasicPOMCP.simulate(planner, rand(MersenneTwister(1), initialstate(pomdp)), node, 20)
end;

@testset "belief dependent actions" begin
	pomdp = ConstObsPOMDP()
	function POMDPs.actions(m::ConstObsPOMDP, b::LeafNodeBelief)
		@test currentobs(b) == true
        @test history(b)[end].o == true
        @test history(b)[end].a == :the_only_action
        @test mean(b) == 0.0
        @test mode(b) == 0.0
        @test only(support(b)) == false
        @test pdf(b, false) == 1.0
        @test pdf(b, true) == 0.0
        @test rand(b) == false
        @test n_particles(b) == 1
        @test only(particles(b)) == false
        @test only(weights(b)) == 1.0
        @test only(weighted_particles(b)) == (false => 1.0)
        @test weight_sum(b) == 1.0
        @test weight(b, 1) == 1.0
        @test particle(b, 1) == false

        # old type name - this can be removed when upgrading versions
        @test b isa AOHistoryBelief
        return actions(m)
	end

	solver = POMCPSolver(rng = MersenneTwister(1))
	planner = solve(solver, pomdp)
	b = initialstate(pomdp)
    tree = BasicPOMCP.POMCPTree(pomdp, b, solver.tree_queries)
    node = BasicPOMCP.POMCPObsNode(tree, 1)

    @inferred BasicPOMCP.simulate(planner, rand(MersenneTwister(1), initialstate(pomdp)), node, 20)
end;

@testset "simulation" begin
	pomdp = BabyPOMDP()
	solver = POMCPSolver(rng = MersenneTwister(1))
	planner = solve(solver, pomdp)
    solver = POMCPSolver(max_time=0.1, tree_queries=typemax(Int), rng = MersenneTwister(1))
    planner = solve(solver, pomdp)
	b = initialstate(pomdp)

    a, info = action_info(planner, b)
    println("time below should be about 0.1 seconds")
    etime = @elapsed a, info = action_info(planner, b)
    @show etime
    @test etime < 0.2
    @show info[:search_time_us]

    sim = HistoryRecorder(max_steps=10)
    simulate(sim, pomdp, planner, updater(pomdp))
end;

@testset "d3t" begin
	pomdp = BabyPOMDP()
    solver = POMCPSolver(max_time=0.1, tree_queries=typemax(Int), rng = MersenneTwister(1))
    planner = solve(solver, pomdp)
	b = initialstate(pomdp)
    a, info = action_info(planner, b, tree_in_info=true)

    d3t = D3Tree(info[:tree], title="test tree")
    # inchrome(d3t)
    show(stdout, MIME("text/plain"), d3t)

    solver = POMCPSolver(max_time=0.1, tree_queries=typemax(Int), rng=MersenneTwister(1), tree_in_info=true)
    planner = solve(solver, pomdp)
    a, info = action_info(planner, b)

    d3t = D3Tree(info[:tree], title="test tree (tree_in_info solver option)")
end;

@testset "Minimal_Example" begin
    @nbinclude(joinpath(dirname(@__FILE__), "..", "notebooks", "Minimal_Example.ipynb"))
end;

@testset "consistency" begin
    # test consistency when rng is specified
    pomdp = BabyPOMDP()
    solver = POMCPSolver(rng = MersenneTwister(1))
    planner = solve(solver, pomdp)
    hist1 = simulate(HistoryRecorder(max_steps=1000, rng=MersenneTwister(3)), pomdp, planner)

    solver = POMCPSolver(rng = MersenneTwister(1))
    planner = solve(solver, pomdp)
    hist2 = simulate(HistoryRecorder(max_steps=1000, rng=MersenneTwister(3)), pomdp, planner)

    @test discounted_reward(hist1) == discounted_reward(hist2)
end;

@testset "requires" begin
    # REQUIREMENTS
    solver = POMCPSolver()
    pomdp = TigerPOMDP()

    println("============== @requirements_info with only solver:")
    requirements_info(solver)
    println("============== @requirements_info with solver and pomdp:")
    requirements_info(solver, pomdp)
    @show_requirements POMDPs.solve(solver, pomdp)
end;

@testset "errors" begin
    struct TerminalPOMDP <: POMDP{Int, Int, Int} end
    POMDPs.isterminal(::TerminalPOMDP, s) = true
    POMDPs.actions(::TerminalPOMDP) = [1,2,3]

    solver = POMCPSolver()
    planner = solve(solver, TerminalPOMDP())
    @test_throws AllSamplesTerminal action(planner, Deterministic(1))
    let ex = nothing
        try
            action(planner, Deterministic(1))
        catch ex
        end

        @test sprint(showerror, ex) == """
        Planner failed to choose an action because all states sampled from the belief were terminal.

        To see the belief, catch this exception as ex and see ex.belief.

        To specify an action for this case, use the default_action solver parameter.
        """
    end
end

ext: jl


text:
# ----- PyTorch Model ----- #
using PyCall
struct PyTorchModel <: AbstractDifferentiableModel
    nn::Any
end

function logits(M::PyTorchModel, X::AbstractArray)
    py"""
    import torch
    from torch import nn
    """
    nn = M.nn
    if !isa(X, Matrix)
      X = reshape(X, length(X), 1)
    end
    ŷ = py"$nn(torch.Tensor($X).T).detach().numpy()"
    ŷ = isa(ŷ, AbstractArray) ? ŷ : [ŷ]
    return ŷ'
end

probs(M::PyTorchModel, X::AbstractArray)= σ.(logits(M, X))


text:
@testset "getindex" begin
    @testset "getindex(::Matrix{<:Number},...)" begin
        x = [1.0 2.0 3.0; 10.0 20.0 30.0]

        @testset "single element" begin
            test_rrule(getindex, x, 2)
            test_rrule(getindex, x, 2, 1)
            test_rrule(getindex, x, 2, 2)

            test_rrule(getindex, x, CartesianIndex(2, 3))
        end

        @testset "slice/index postions" begin
            test_rrule(getindex, x, 2:3)
            test_rrule(getindex, x, 3:-1:2)
            test_rrule(getindex, x, [3, 2])
            test_rrule(getindex, x, [2, 3])

            test_rrule(getindex, x, 1:2, 2:3)
            test_rrule(getindex, x, (:), 2:3)

            test_rrule(getindex, x, 1:2, 1)
            test_rrule(getindex, x, 1, 1:2)

            test_rrule(getindex, x, 1:2, 2:3)
            test_rrule(getindex, x, (:), 2:3)

            test_rrule(getindex, x, (:), (:))
            test_rrule(getindex, x, (:))
        end

        @testset "masking" begin
            test_rrule(getindex, x, trues(size(x)))
            test_rrule(getindex, x, trues(length(x)))

            mask = falses(size(x))
            mask[2, 3] = true
            mask[1, 2] = true
            test_rrule(getindex, x, mask)

            test_rrule(getindex, x, [true, false], (:))
        end

        @testset "By position with repeated elements" begin
            test_rrule(getindex, x, [2, 2])
            test_rrule(getindex, x, [2, 2, 2])
            test_rrule(getindex, x, [2, 2], [3, 3])
        end
    end
end


text:
@testset "Test Parser" begin
    @testset "Empty program" begin
        @test isnothing(m.parse(""))
    end

    @testset "Test Parsing Let Statements" begin
        for (code, expected_ident, expected_value) in [
            ("let x = 5;", "x", 5),
            ("let y = true;", "y", true),
            ("let foobar = y;", "foobar", "y"),
            ("let a = null;", "a", nothing),
        ]
            test_let_statement(ls::m.LetStatement, name::String) =
                ls.name.value == name && m.token_literal(ls.name) == name

            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1

            statement = program.statements[1]
            test_let_statement(statement, expected_ident)

            val = statement.value
            test_literal_expression(val, expected_value)
        end
    end

    @testset "Test Parsing Return Statements" begin
        for (code, expected_value) in
            [("return 5;", 5), ("return false;", false), ("return y;", "y")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1

            statement = program.statements[1]
            test_literal_expression(statement.return_value, expected_value)
        end
    end

    @testset "Test Parsing Identifier Expression" begin
        for (code, value) in [("foobar;", "foobar")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            ident = statement.expression
            test_literal_expression(ident, value)
        end
    end

    @testset "Test Parsing BooleanLiteral Expression" begin
        for (code, value) in [("true;", true), ("false;", false)]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            bool = statement.expression
            test_literal_expression(bool, value)
        end
    end

    @testset "Test Parsing Integer Literal Expression" begin
        for (code, value) in [("5;", 5)]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            il = statement.expression
            test_literal_expression(il, value)
        end
    end

    @testset "Test Parsing Prefix Expressions" begin
        for (code, operator, right_value) in [
            ("!5", "!", 5),
            ("-15", "-", 15),
            ("-a", "-", "a"),
            ("!true", "!", true),
            ("!false", "!", false),
        ]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            expr = statement.expression
            @test isa(expr, m.PrefixExpression)
            @test expr.operator == operator

            test_literal_expression(expr.right, right_value)
        end
    end

    @testset "Test Parsing Infix Expressions" begin
        for (code, left_value, operator, right_value) in [
            ("5 + 5", 5, "+", 5),
            ("5 - 5", 5, "-", 5),
            ("5 * 5", 5, "*", 5),
            ("5 / 5", 5, "/", 5),
            ("5 > 5", 5, ">", 5),
            ("5 < 5", 5, "<", 5),
            ("5 == 5", 5, "==", 5),
            ("5 != 5", 5, "!=", 5),
            ("true == true", true, "==", true),
            ("true != false", true, "!=", false),
            ("false == false", false, "==", false),
        ]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            expr = statement.expression
            test_infix_expression(expr, left_value, operator, right_value)
        end
    end

    @testset "Test Parsing If Expression" begin
        for (code) in [("if (x < y) { x }")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            expr = statement.expression

            @test isa(expr, m.IfExpression)

            test_infix_expression(expr.condition, "x", "<", "y")

            @test length(expr.consequence.statements) == 1

            consequence = expr.consequence.statements[1]
            @test isa(consequence, m.ExpressionStatement)

            test_identifier(consequence.expression, "x")

            @test isnothing(expr.alternative)
        end
    end

    @testset "Test Parsing If Else Expression" begin
        for (code) in [("if (x < y) { x } else { y }")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            expr = statement.expression

            @test isa(expr, m.IfExpression)

            test_infix_expression(expr.condition, "x", "<", "y")

            @test length(expr.consequence.statements) == 1

            consequence = expr.consequence.statements[1]
            @test isa(consequence, m.ExpressionStatement)

            test_identifier(consequence.expression, "x")

            alternative = expr.alternative.statements[1]
            @test isa(alternative, m.ExpressionStatement)

            test_identifier(alternative.expression, "y")
        end
    end

    @testset "Test Parsing While Statement" begin
        for (code) in [("while (x < y) { x = x + 1 }")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.WhileStatement)

            statement = program.statements[1]

            @test isa(statement, m.WhileStatement)

            test_infix_expression(statement.condition, "x", "<", "y")

            @test length(statement.body.statements) == 1

            ls = statement.body.statements[1]
            @test isa(ls, m.LetStatement)
            @test ls.reassign

            test_identifier(ls.name, "x")
            test_infix_expression(ls.value, "x", "+", 1)
        end
    end

    @testset "Test Parsing Break Statement" begin
        for code in ["break", "break;"]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.BreakStatement)
        end
    end

    @testset "Test Parsing Continue Statement" begin
        for code in ["continue", "continue;"]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ContinueStatement)
        end
    end

    @testset "Test Parsing Functional Literal" begin
        for (code) in [("fn(x, y) { x + y; }")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)
            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            fn = statement.expression

            @test isa(fn, m.FunctionLiteral)
            @test length(fn.parameters) == 2

            test_literal_expression(fn.parameters[1], "x")
            test_literal_expression(fn.parameters[2], "y")

            @test length(fn.body.statements) == 1

            body_statement = fn.body.statements[1]
            @test isa(body_statement, m.ExpressionStatement)
            test_infix_expression(body_statement.expression, "x", "+", "y")
        end
    end

    @testset "Test Parsing Functional Literal With Name" begin
        for (code) in [("let myFunction = fn() { };")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.LetStatement)

            fl = program.statements[1].value
            @test fl.name == "myFunction"
        end
    end

    @testset "Test Parsing Function Parameters" begin
        for (code, expected) in
            [("fn() {};", []), ("fn(x) {};", ["x"]), ("fn(x, y, z) {};", ["x", "y", "z"])]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)
            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            fn = statement.expression

            @test length(fn.parameters) == length(expected)

            for (parameter, expected_parameter) in zip(fn.parameters, expected)
                test_literal_expression(parameter, expected_parameter)
            end
        end
    end

    @testset "Test Parsing Call Expression" begin
        for (code) in [("add(1, 2 * 3, 4 + 5)")]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)
            check_parser_errors(p)

            @test length(program.statements) == 1
            @test isa(program.statements[1], m.ExpressionStatement)

            statement = program.statements[1]
            expr = statement.expression
            @test isa(expr, m.CallExpression)
            test_identifier(expr.fn, "add")

            @test length(expr.arguments) == 3
            test_literal_expression(expr.arguments[1], 1)
            test_infix_expression(expr.arguments[2], 2, "*", 3)
            test_infix_expression(expr.arguments[3], 4, "+", 5)
        end
    end

    @testset "Test Parsing String Literal Expression" begin
        code = "\"hello world\""
        l = m.Lexer(code)
        p = m.Parser(l)
        program = m.parse!(p)

        check_parser_errors(p)

        expr = program.statements[1].expression
        @test isa(expr, m.StringLiteral)
        @test expr.value == "hello world"
    end

    @testset "Test Parsing Array Literal" begin
        code = "[1, 2 * 2, 3 + 3]"
        l = m.Lexer(code)
        p = m.Parser(l)
        program = m.parse!(p)

        check_parser_errors(p)

        arr = program.statements[1].expression
        @test isa(arr, m.ArrayLiteral)
        @test length(arr.elements) == 3

        test_integer_literal(arr.elements[1], 1)
        test_infix_expression(arr.elements[2], 2, "*", 2)
        test_infix_expression(arr.elements[3], 3, "+", 3)
    end

    @testset "Test Parsing Index Expression" begin
        code = "myArray[1 + 1]"
        l = m.Lexer(code)
        p = m.Parser(l)
        program = m.parse!(p)
        check_parser_errors(p)

        expr = program.statements[1].expression
        @test isa(expr, m.IndexExpression)

        test_identifier(expr.left, "myArray")
        test_infix_expression(expr.index, 1, "+", 1)
    end

    @testset "Test Parsing Hash Literal" begin
        @testset "Test Parsing Hash Literal with String Keys" begin
            code = """{"one": 1, "two": 2, "three": 3}"""
            expected = Dict("one" => 1, "two" => 2, "three" => 3)
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)
            check_parser_errors(p)

            hash = program.statements[1].expression
            @test isa(hash, m.HashLiteral)
            @test length(hash.pairs) == 3

            for (key, value) in hash.pairs
                @test isa(key, m.StringLiteral)
                @test key.value ∈ keys(expected)
                test_integer_literal(value, expected[key.value])
            end
        end

        @testset "Test Parsing Hash Literal with String Keys and Expression Values" begin
            code = """{"one": 0 + 1, "two": 10 - 8, "three": 15 / 5}"""
            tests = Dict(
                "one" => x -> test_infix_expression(x, 0, "+", 1),
                "two" => x -> test_infix_expression(x, 10, "-", 8),
                "three" => x -> test_infix_expression(x, 15, "/", 5),
            )

            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            hash = program.statements[1].expression
            @test isa(hash, m.HashLiteral)

            @test length(hash.pairs) == 3

            for (key, value) in hash.pairs
                @test isa(key, m.StringLiteral)

                @test key.value ∈ keys(tests)

                tests[key.value](value)
            end
        end

        @testset "Test Parsing Hash Literal with Integer Keys" begin
            code = "{1: 1, 2: 2, 3: 3}"
            expected = Dict(1 => 1, 2 => 2, 3 => 3)
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            hash = program.statements[1].expression
            @test isa(hash, m.HashLiteral)

            @test length(hash.pairs) == 3

            for (key, value) in hash.pairs
                @test isa(key, m.IntegerLiteral)

                @test key.value ∈ keys(expected)

                test_integer_literal(value, expected[key.value])
            end
        end

        @testset "Test Parsing Hash Literal with Boolean Keys" begin
            code = "{false: 0, true: 1}"
            expected = Dict(false => 0, true => 1)
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            hash = program.statements[1].expression
            @test isa(hash, m.HashLiteral)

            @test length(hash.pairs) == 2

            for (key, value) in hash.pairs
                @test isa(key, m.BooleanLiteral)

                @test key.value ∈ keys(expected)

                test_integer_literal(value, expected[key.value])
            end
        end

        @testset "Test Parsing Empty Hash Literal" begin
            code = "{}"
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            hash = program.statements[1].expression
            @test isa(hash, m.HashLiteral)

            @test length(hash.pairs) == 0
        end
    end

    @testset "Test Parsing Macro Literal" begin
        code = "macro(x, y) { x + y; }"
        l = m.Lexer(code)
        p = m.Parser(l)
        program = m.parse!(p)

        check_parser_errors(p)

        ml = program.statements[1].expression
        @test isa(ml, m.MacroLiteral)

        @test length(ml.parameters) == 2

        test_literal_expression(ml.parameters[1], "x")
        test_literal_expression(ml.parameters[2], "y")

        @test length(ml.body.statements) == 1

        test_infix_expression(ml.body.statements[1].expression, "x", "+", "y")
    end

    @testset "Test Operator Order" begin
        for (code, expected) in [
            ("-a * b", "((-a) * b)"),
            ("!-a", "(!(-a))"),
            ("a + b + c", "((a + b) + c)"),
            ("a + b - c", "((a + b) - c)"),
            ("a * b * c", "((a * b) * c)"),
            ("a * b / c", "((a * b) / c)"),
            ("a + b / c", "(a + (b / c))"),
            ("a + b * c + d / e - f", "(((a + (b * c)) + (d / e)) - f)"),
            ("3 + 4; -5 * 5", "(3 + 4)((-5) * 5)"), #FIXME
            ("5 > 4 == 3 < 4", "((5 > 4) == (3 < 4))"),
            ("5 < 4 != 3 > 4", "((5 < 4) != (3 > 4))"),
            ("3 + 4 * 5 == 3 * 1 + 4 * 5", "((3 + (4 * 5)) == ((3 * 1) + (4 * 5)))"),
            ("true", "true"),
            ("false", "false"),
            ("3 > 5 == false", "((3 > 5) == false)"),
            ("3 < 5 == true", "((3 < 5) == true)"),
            ("1 + (2 + 3) + 4", "((1 + (2 + 3)) + 4)"),
            ("(5 + 5) * 2", "((5 + 5) * 2)"),
            ("2 / (5 + 5)", "(2 / (5 + 5))"),
            ("-(5 + 5)", "(-(5 + 5))"),
            ("!(true == true)", "(!(true == true))"),
            ("a + add(b * c) + d", "((a + add((b * c))) + d)"),
            (
                "add(a, b, 1, 2 * 3, 4 + 5, add(6, 7 * 8))",
                "add(a, b, 1, (2 * 3), (4 + 5), add(6, (7 * 8)))",
            ),
            ("add(a + b + c * d / f + g)", "add((((a + b) + ((c * d) / f)) + g))"),
            ("a * [1, 2, 3, 4][b * c] * d", "((a * ([1, 2, 3, 4][(b * c)])) * d)"),
            (
                "add(a * b[2], b[1], 2 * [1, 2][1])",
                "add((a * (b[2])), (b[1]), (2 * ([1, 2][1])))",
            ),
        ]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)

            check_parser_errors(p)

            @test string(program) == expected
        end
    end

    @testset "Test Parsing Errors" begin
        for (code, expected) in [
            (
                "100000000000000000000",
                ["ERROR: parse error: could not parse 100000000000000000000 as integer"],
            ),
            ("2!", ["ERROR: parse error: no prefix parse function for EOF found"]),
            (
                "{1:2",
                ["ERROR: parse error: expected next token to be COMMA, got EOF instead"],
            ),
            (
                "{1}",
                [
                    "ERROR: parse error: expected next token to be COLON, got RBRACE instead",
                    "ERROR: parse error: no prefix parse function for RBRACE found",
                ],
            ),
            (
                "[1",
                ["ERROR: parse error: expected next token to be RBRACKET, got EOF instead"],
            ),
            (
                "{1:2}[1",
                ["ERROR: parse error: expected next token to be RBRACKET, got EOF instead"],
            ),
            (
                "fn x",
                ["ERROR: parse error: expected next token to be LPAREN, got IDENT instead"],
            ),
            (
                "fn (x {x}",
                [
                    "ERROR: parse error: expected next token to be RPAREN, got LBRACE instead",
                ],
            ),
            (
                "fn (x) x",
                ["ERROR: parse error: expected next token to be LBRACE, got IDENT instead"],
            ),
            ("fn (x) {x", ["ERROR: parse error: braces must be closed"]),
            (
                "if x",
                ["ERROR: parse error: expected next token to be LPAREN, got IDENT instead"],
            ),
            (
                "if (x",
                ["ERROR: parse error: expected next token to be RPAREN, got EOF instead"],
            ),
            (
                "if (x) c",
                ["ERROR: parse error: expected next token to be LBRACE, got IDENT instead"],
            ),
            (
                "if (x) { 1 } else 2",
                ["ERROR: parse error: expected next token to be LBRACE, got INT instead"],
            ),
            ("if (true) {", ["ERROR: parse error: braces must be closed"]),
            (
                "while",
                ["ERROR: parse error: expected next token to be LPAREN, got EOF instead"],
            ),
            (
                "while (true",
                ["ERROR: parse error: expected next token to be RPAREN, got EOF instead"],
            ),
            (
                "while (true)",
                ["ERROR: parse error: expected next token to be LBRACE, got EOF instead"],
            ),
            ("while (true) {", ["ERROR: parse error: braces must be closed"]),
            (
                "let 5",
                ["ERROR: parse error: expected next token to be IDENT, got INT instead"],
            ),
            (
                "let x 3",
                ["ERROR: parse error: expected next token to be ASSIGN, got INT instead"],
            ),
            (
                "macro x",
                ["ERROR: parse error: expected next token to be LPAREN, got IDENT instead"],
            ),
            (
                "macro (x",
                [
                    "ERROR: parse error: expected next token to be RPAREN, got EOF instead",
                    "ERROR: parse error: expected next token to be LBRACE, got EOF instead",
                ],
            ),
            (
                "macro (x)",
                ["ERROR: parse error: expected next token to be LBRACE, got EOF instead"],
            ),
            ("macro (x) {x", ["ERROR: parse error: braces must be closed"]),
        ]
            l = m.Lexer(code)
            p = m.Parser(l)
            program = m.parse!(p)
            @test map(string, p.errors) == expected
        end
    end

    @testset "Test Directly Parsing" begin
        for (code, expected_output) in [
            (
                "100000000000000000000",
                "ERROR: parser has 1 error\nERROR: parse error: could not parse 100000000000000000000 as integer\n",
            ),
        ]
            c = IOCapture.capture() do
                program = m.parse(code)
                @test isnothing(program)
            end

            @test c.output == expected_output
        end
    end
end


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity

using OceanTurbulenceParameterEstimation
using OceanTurbulenceParameterEstimation.Observations: FieldTimeSeriesCollector, initialize_simulation!, observation_times
using OceanTurbulenceParameterEstimation.InverseProblems: transpose_model_output

using Distributions

@testset "Forward map tests" begin

    Nz = 8
    experiment_name = "forward_map_test"

    # Generate synthetic observations

    default_closure() = ConvectiveAdjustmentVerticalDiffusivity(; convective_κz=1.0,
                                                                  background_κz=1e-5,
                                                                  convective_νz=0.9,
                                                                  background_νz=1e-4)

    function build_simulation(size=Nz)
        N² = 1e-5
        grid = RegularRectilinearGrid(size=size, z=(-128, 0), topology=(Flat, Flat, Bounded))
        u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(-1e-4))
        b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(1e-7), bottom = GradientBoundaryCondition(N²))
        model = HydrostaticFreeSurfaceModel(grid = grid,
                                            tracers = :b,
                                            buoyancy = BuoyancyTracer(),
                                            boundary_conditions = (; u=u_bcs, b=b_bcs),
                                            closure = default_closure())

        set!(model, b = (x, y, z) -> N² * z)

        simulation = Simulation(model; Δt=20.0, stop_iteration=3*60*10)

        return simulation
    end

    # Make observations

    truth_simulation = build_simulation()
    model = truth_simulation.model
    stop_iteration = truth_simulation.stop_iteration
    truth_simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                                schedule = IterationInterval(round(Int, stop_iteration / 10)),
                                                                prefix = experiment_name,
                                                                array_type = Array{Float64},
                                                                field_slicer = nothing,
                                                                force = true)

    run!(truth_simulation)

    data_path = experiment_name * ".jld2"
    observations = OneDimensionalTimeSeries(data_path, field_names=(:u, :b), normalize=ZScore)

    #####
    ##### Make model data
    #####

    @testset "Single-member forward map" begin
        # First test transpose_model_output
        test_simulation = build_simulation()
        collected_fields = (u = test_simulation.model.velocities.u, b = test_simulation.model.tracers.b)

        time_series_collector = FieldTimeSeriesCollector(collected_fields, observation_times(observations))

        initialize_simulation!(test_simulation, observations, time_series_collector)
        run!(test_simulation)

        output = transpose_model_output(time_series_collector, observations)[1]

        test_b = output.field_time_serieses.b
        truth_b = observations.field_time_serieses.b
        test_u = output.field_time_serieses.u
        truth_u = observations.field_time_serieses.u

        @test interior(test_b) == interior(truth_b)
        @test interior(test_u) == interior(truth_u)

        # Now test forward_map and observation_map
        closure = default_closure()

        priors = (
            convective_κz = Normal(closure.convective_κz, 0.05),
            background_κz = Normal(closure.background_κz, 1e-5),
            convective_νz = Normal(closure.convective_νz, 0.01),
            background_νz = Normal(closure.background_νz, 1e-5),
        )

        free_parameters = FreeParameters(priors)
        calibration = InverseProblem(observations, ensemble_simulation, free_parameters)
        optimal_parameters = [getproperty(closure, p) for p in keys(priors)]

        x = forward_map(calibration, optimal_parameters)[1:1, :]
        y = observation_map(calibration)

        @test x == y
    end

    @testset "Two-member transposition of model output" begin
        ensemble_size = ColumnEnsembleSize(Nz=Nz, ensemble=(2, 1), Hz=1)
        test_simulation = build_simulation(ensemble_size)
        collected_fields = (u = test_simulation.model.velocities.u, b = test_simulation.model.tracers.b)
        time_series_collector = FieldTimeSeriesCollector(collected_fields, observation_times(observations))
        initialize_simulation!(test_simulation, observations, time_series_collector)
        run!(test_simulation)
        output = transpose_model_output(time_series_collector, observations)

        truth_b = observations.field_time_serieses.b
        truth_u = observations.field_time_serieses.u

        test_b1 = output[1].field_time_serieses.b
        test_u1 = output[1].field_time_serieses.u
        test_b2 = output[2].field_time_serieses.b
        test_u2 = output[2].field_time_serieses.u

        @test interior(test_b1) == interior(truth_b)
        @test interior(test_u1) == interior(truth_u)
        @test interior(test_b2) == interior(truth_b)
        @test interior(test_u2) == interior(truth_u)
    end

    @testset "Four-member transposition of model output" begin
        ensemble_size = ColumnEnsembleSize(Nz=Nz, ensemble=(2, 2), Hz=1)
        test_simulation = build_simulation(ensemble_size)
        collected_fields = (u = test_simulation.model.velocities.u, b = test_simulation.model.tracers.b)
        time_series_collector = FieldTimeSeriesCollector(collected_fields, observation_times(observations))
        initialize_simulation!(test_simulation, observations, time_series_collector)
        run!(test_simulation)
        output = transpose_model_output(time_series_collector, observations)

        truth_b = observations.field_time_serieses.b
        truth_u = observations.field_time_serieses.u

        test_b1 = output[1].field_time_serieses.b
        test_u1 = output[1].field_time_serieses.u
        test_b2 = output[2].field_time_serieses.b
        test_u2 = output[2].field_time_serieses.u
        test_b3 = output[3].field_time_serieses.b
        test_u3 = output[3].field_time_serieses.u
        test_b4 = output[4].field_time_serieses.b
        test_u4 = output[4].field_time_serieses.u

        @test interior(test_b1) == interior(truth_b)
        @test interior(test_u1) == interior(truth_u)
        @test interior(test_b2) == interior(truth_b)
        @test interior(test_u2) == interior(truth_u)
        @test interior(test_b3) == interior(truth_b)
        @test interior(test_u3) == interior(truth_u)
        @test interior(test_b4) == interior(truth_b)
        @test interior(test_u4) == interior(truth_u)
    end
end


text:
module PIMDtools
    using LinearAlgebra
    #using HDF5
    export Trajectory,a0,get_relative_positions,calc_MSD,toBinary_r,
            calc_MSD_fromfile
    const a0 = 0.5291772 #Bohr length [Angstrom]

    const nbytes_Int64= sizeof(Int64)
    const nbytes_Float64= sizeof(Float64)

    struct Atom_vFE
        vx::Float64
        vy::Float64
        vz::Float64
        Fx::Float64
        Fy::Float64
        Fz::Float64
        E::Float64
    end

    struct Atom
        x::Float64 #[a.u.]
        y::Float64
        z::Float64
        vFE::Union{Nothing,Atom_vFE}

        Atom(A) = Atom(A,ronly=true)

        function Atom(A;ronly=true)
            if ronly
                return new(A[1],A[2],A[3],nothing)
            else
                return new(A[1],A[2],A[3],Atom_vFE(A[4],A[5],A[6],A[7],A[8],A[9],A[10]))
            end
        end
        #Atom(A) = new(A[1],A[2],A[3],A[4],A[5],A[6],A[7],A[8],A[9],A[10])
    end

    struct Snapshot
        atoms::Array{Atom,1}
        Snapshot(atoms) = new(atoms)
    end

    struct Trajectory
        data::Array{Snapshot,1}
        numtrj::Int64
        numatoms::Int64
        atomnames::Array{String,1}
        numeachatom::Array{Int64,1}
        periodic::Bool
        unitcell::Union{Nothing,Array{Array{Float64,1},1}}
        rotatematrix::Union{Nothing,Array{Float64,2}}
        rotatematrix_inv::Union{Nothing,Array{Float64,2}}
        itrj_start::Int64
    end

    Base.length(trj::Trajectory) = trj.numtrj

    Base.length(snapshot::Snapshot) = length(snapshot.atoms)

    function Base.getindex(x::Snapshot,i)
        return x.atoms[i]
    end

    function Base.getindex(x::Trajectory,i)
        return x.data[i]
    end

    function Base.getindex(trj::Trajectory,iatom,itrj)
        #println(trj[i][j].x)
        return trj[itrj][iatom]
    end

    #atominfo::Array{Tuple,1}
    #like: [("Pd",20),("Ru",16),("Al",92)] for Pd20Ru16Al92

    function Trajectory(filename,itrj_start,itrj_end,atominfo;
                            periodic=true,
                            unitcell=nothing,
                            ronly=true
                        )

        numkinds = length(atominfo)
        print("System: ")
        atomnames = String[]
        numeachatom = Int64[]
        df = zeros(Float64,10)
        numtrj = itrj_end-itrj_start+1
        for (kind,num) in atominfo
            print("$(kind)$(num)")
            push!(atomnames,kind)
            push!(numeachatom,num)
        end
        println("\t")

        if periodic
            @assert unitcell != nothing "periocic system should have unitcell!"
            rotatematrix = zeros(Float64,3,3)
            for i=1:3
                rotatematrix[:,i] = unitcell[i]/norm(unitcell[i])
            end
            rotatematrix_inv = inv(rotatematrix)

            println("unitcell is ")
            for i=1:3
                println(unitcell[i]," [a.u.]")
            end
            println("rotatematrix is ",rotatematrix)
        else
            unitcell = nothing
            rotatematrix = nothing
            rotatematrix_inv = nothing
        end

        numatoms = sum(numeachatom)
        println("number of total atoms: ",numatoms)
        println("number of total data: ",numtrj)

        data = Array{Snapshot,1}(undef,numtrj)
        atoms = Array{Atom,1}(undef,numatoms)

        fp = open(filename,"r")
        for itrj=1:itrj_start-1
            for iatom=1:numatoms
                readline(fp)
            end
        end

        count = 0
        for itrj=itrj_start:itrj_end
            count += 1
            for iatom=1:numatoms
                u = split(readline(fp))
                di = parse(Int64,u[1])
                for i=2:11
                    df[i-1] = parse(Float64,u[i])
                end
                atoms[iatom] = Atom(df[1:10],ronly=ronly)
            end
            data[count] = deepcopy(Snapshot(atoms))
        end
        close(fp)

        return Trajectory(
            data,
            numtrj,
            numatoms,
            atomnames,
            numeachatom,
            periodic,
            unitcell,
            rotatematrix,
            rotatematrix_inv,
            itrj_start
        )
        #=
            data::Array{Snapshot,1}
            numtrj::Int64
            numatoms::Int64
            atomnames::Array{String,1}
            numeachatom::Array{Int64,1}
            periodic::Bool
            unitcell::Union{Nothing,Array{Array{Float64,1},1}}
            rotatematrix::Union{Nothing,Array{Float64,2}}
            rotatematrix_inv::Union{Nothing,Array{Float64,2}}
            itrj_start::Int64
        =#
    end

    """
        Binary:
        num_of_eachdata::Int64
        for itrj=start_itrj:start_itrj+num_of_eachdata-1
            itrj::Int64
            numatoms::Int64
            for iatom=1:numatoms
                x::Float64
                y::Float64
                z::Float64
            end
        end
    """
    function toBinary_r(filename,dirname,itrj_start,itrj_end,num_of_structures,atominfo)
        numkinds = length(atominfo)
        print("System: ")
        atomnames = String[]
        numeachatom = Int64[]
        df = zeros(Float64,10)
        numtrj = itrj_end-itrj_start+1
        for (kind,num) in atominfo
            print("$(kind)$(num)")
            push!(atomnames,kind)
            push!(numeachatom,num)
        end
        println("\t")

        numatoms = sum(numeachatom)
        println("number of total atoms: ",numatoms)
        println("number of total data: ",numtrj)

        fp = open(filename,"r")
        for itrj=1:itrj_start-1
            for iatom=1:numatoms
                readline(fp)
            end
        end

        count = 0
        setcount = itrj_start
        countstruct = 0
        setstart = setcount
        setend  =setcount +num_of_structures - 1

        dataname = dirname*"/$(lpad(setstart,8,"0"))_$(lpad(setend,8,"0"))_trjset_r.dat"
        #ifp1 = h5open(dataname, "w")
        ifp1 = open(dataname, "w")
        num_of_eachdata = setend -setstart + 1
        write(ifp1,num_of_eachdata)
        #ifp1 = open(dirname*"/$(lpad(setstart,8,"0"))_$(lpad(setend,8,"0"))_trjset_r.dat","w")
        for itrj=itrj_start:itrj_end
            count += 1
            countstruct += 1
            write(ifp1,itrj)
            write(ifp1,numatoms)
            for iatom=1:numatoms
                u = split(readline(fp))
                di = parse(Int64,u[1])
                for i=2:11
                    df[i-1] = parse(Float64,u[i])
                end
                #atoms[iatom] = Atom(df[1:10],ronly=ronly)
                write(ifp1,df[1])
                write(ifp1,df[2])
                write(ifp1,df[3])
                #write(ifp1,"$(lpad(itrj,8,"0"))/$(lpad(iatom,4,"0"))/r",Float64[df[1],df[2],df[3]])
            end
            setcount += 1
            if countstruct == num_of_structures
                println("$dataname is done. $num_of_eachdata strctures")
                close(ifp1)
                setstart = setcount

                if setstart <= itrj_end
                    setend  =setcount +num_of_structures - 1
                    if setend > itrj_end
                        setend =itrj_end
                    end
                    num_of_eachdata = setend -setstart + 1

                    dataname = dirname*"/$(lpad(setstart,8,"0"))_$(lpad(setend,8,"0"))_trjset_r.dat"
                    #ifp1 = h5open(dataname, "w")
                    ifp1 = open(dataname, "w")
                    write(ifp1,num_of_eachdata)
                    countstruct = 0
                end
            end
        end
        println("$dataname is done. $num_of_eachdata strctures")
        close(ifp1)
        close(fp)

    end

    function toBinary_r(filename,dirname,itrj_start,itrj_end,num_of_structures,atominfo,
        unitcell
        )

        numkinds = length(atominfo)
        print("System: ")
        atomnames = String[]
        numeachatom = Int64[]
        df = zeros(Float64,10)
        numtrj = itrj_end-itrj_start+1
        for (kind,num) in atominfo
            print("$(kind)$(num)")
            push!(atomnames,kind)
            push!(numeachatom,num)
        end
        println("\t")

        rotatematrix = zeros(Float64,3,3)
        for i=1:3
            rotatematrix[:,i] = unitcell[i]/norm(unitcell[i])
        end
        rotatematrix_inv = inv(rotatematrix)

        println("unitcell is ")
        for i=1:3
            println(unitcell[i]," [a.u.]")
        end
        println("rotatematrix is ",rotatematrix)

        unorm = zeros(Float64,3)
        dr = zeros(Float64,3)
        tempvec = zeros(Float64,3)
        dur = zeros(Float64,3)
        for i=1:3
            unorm[i] = norm(unitcell[i])
        end

        numatoms = sum(numeachatom)
        println("number of total atoms: ",numatoms)
        println("number of total data: ",numtrj)

        fp = open(filename,"r")
        for itrj=1:itrj_start-1
            for iatom=1:numatoms
                readline(fp)
            end
        end

        count = 0
        setcount = itrj_start
        countstruct = 0
        setstart = setcount
        setend  =setcount +num_of_structures - 1

        df_struct = zeros(Float64,3,numatoms)
        df_struct_old = zeros(Float64,3,numatoms)
        dataname = dirname*"/$(lpad(setstart,8,"0"))_$(lpad(setend,8,"0"))_trjset_r.dat"
        #ifp1 = h5open(dataname, "w")
        ifp1 = open(dataname, "w")
        num_of_eachdata = setend -setstart + 1
        write(ifp1,num_of_eachdata)
        #ifp1 = open(dirname*"/$(lpad(setstart,8,"0"))_$(lpad(setend,8,"0"))_trjset_r.dat","w")
        for itrj=itrj_start:itrj_end
            count += 1
            countstruct += 1
            write(ifp1,itrj)
            write(ifp1,numatoms)
            for iatom=1:numatoms
                u = split(readline(fp))
                di = parse(Int64,u[1])
                for i=2:11
                    df[i-1] = parse(Float64,u[i])
                end
                for i=1:3
                    df_struct[i,iatom] = df[i]
                end

                if itrj > 1
                    dx = df_struct[1,iatom]-df_struct_old[1,iatom]
                    dy = df_struct[2,iatom]-df_struct_old[2,iatom]
                    dz = df_struct[3,iatom]-df_struct_old[3,iatom]

                    dr[1] = dx
                    dr[2] = dy
                    dr[3] = dz
                    mul!(dur,rotatematrix_inv,dr)

                    for i=1:3
                        tempvec .= 0
                        tempvec[i] = ifelse(dur[i] > unorm[i]/3,-unorm[i],0)
                        mul!(dr,rotatematrix,tempvec)
                        df_struct[:,iatom] .+= dr

                        tempvec .= 0
                        tempvec[i] = ifelse(dur[i] < -unorm[i]/3,unorm[i],0)
                        mul!(dr,rotatematrix,tempvec)
                        df_struct[:,iatom] .+= dr
                    end
                end
            end

            df_struct_old[:,:] = df_struct

            for iatom=1:numatoms
                #atoms[iatom] = Atom(df[1:10],ronly=ronly)
                for i=1:3
                    write(ifp1,df_struct[i,iatom])
                    #df_struct[i,iatom] = df[i]
                end
                #write(ifp1,"$(lpad(itrj,8,"0"))/$(lpad(iatom,4,"0"))/r",Float64[df[1],df[2],df[3]])
            end
            setcount += 1
            if countstruct == num_of_structures
                println("$dataname is done. $num_of_eachdata strctures")
                close(ifp1)
                setstart = setcount

                if setstart <= itrj_end
                    setend  =setcount +num_of_structures - 1
                    if setend > itrj_end
                        setend =itrj_end
                    end
                    num_of_eachdata = setend -setstart + 1

                    dataname = dirname*"/$(lpad(setstart,8,"0"))_$(lpad(setend,8,"0"))_trjset_r.dat"
                    #ifp1 = h5open(dataname, "w")
                    ifp1 = open(dataname, "w")
                    write(ifp1,num_of_eachdata)
                    countstruct = 0
                end
            end
        end
        println("$dataname is done. $num_of_eachdata strctures")
        close(ifp1)
        close(fp)

    end

    function get_relative_positions(trj::Trajectory)
        positions = get_positions(trj)
        centerofmass = calc_centerofmass(trj,positions)
        for itrj=1:length(trj)
            for iatom=1:trj.numatoms
                positions[:,iatom,itrj] -= centerofmass[:,itrj]
            end
        end
        return positions
    end

    function calc_Galilei_frame!(trj::Trajectory,positions)
        centerofmass = calc_centerofmass(trj,positions)
        for itrj=1:length(trj)
            for iatom=1:trj.numatoms
                positions[:,iatom,itrj] -= centerofmass[:,itrj]
            end
        end
    end

    function get_relative_positions(positions_old,snapshot::Snapshot,
        unitcell,rotatematrix,rotatematrix_inv)
        positions,numatoms = get_positions(positions_old,snapshot,
                        unitcell,rotatematrix,rotatematrix_inv)
        centerofmass = calc_centerofmass(positions,numatoms)
        for iatom=1:numatoms
            for i=1:3
                positions[i,iatom] -= centerofmass[i]
            end
        end
        return positions
    end

    function get_relative_positions(snapshot::Snapshot)
        numatoms = length(snapshot)
        positions = zeros(Float64,3,numatoms)

        for iatom=1:numatoms
            positions[1,iatom] = snapshot[iatom].x
            positions[2,iatom] = snapshot[iatom].y
            positions[3,iatom] = snapshot[iatom].z
        end

        centerofmass = calc_centerofmass(positions,numatoms)
        for iatom=1:numatoms
            for i=1:3
                positions[i,iatom] -= centerofmass[i]
            end
        end
        return positions
    end

    function calc_centerofmass(trj::Trajectory,positions )
        centerofmass = zeros(Float64,3,length(trj))
        for itrj=1:length(trj)
            for iatom=1:trj.numatoms
                for i=1:3
                    centerofmass[i,itrj] += positions[i,iatom,itrj]/trj.numatoms
                end
            end
        end
        return centerofmass
    end

    function calc_centerofmass(positions,numatoms)
        centerofmass = zeros(Float64,3)
        for iatom=1:numatoms
            for i=1:3
                centerofmass[i] += positions[i,iatom]/numatoms
            end
        end
        return centerofmass
    end

    #=
    function get_positions(positions_old,snapshot::Snapshot,
        unitcell,rotatematrix,rotatematrix_inv)

        unorm = zeros(Float64,3)
        dr = zeros(Float64,3)
        tempvec = zeros(Float64,3)
        dur = zeros(Float64,3)
        for i=1:3
            unorm[i] = norm(unitcell[i])
        end

        numatoms = length(snapshot)
        positions = zeros(Float64,3,numatoms)

        for iatom=1:numatoms
            positions[1,iatom] = snapshot[iatom].x
            positions[2,iatom] = snapshot[iatom].y
            positions[3,iatom] = snapshot[iatom].z
            dx = snapshot[iatom].x - positions_old[1,iatom]
            dy = snapshot[iatom].y - positions_old[2,iatom]
            dz = snapshot[iatom].z - positions_old[3,iatom]

            dr[1] = dx
            dr[2] = dy
            dr[3] = dz
            mul!(dur,rotatematrix_inv,dr)

            for i=1:3
                tempvec .= 0
                tempvec[i] = ifelse(dur[i] > unorm[i]/3,-unorm[i],0)
                mul!(dr,rotatematrix,tempvec)
                positions[:,iatom] .+= dr

                tempvec .= 0
                tempvec[i] = ifelse(dur[i] < -unorm[i]/3,unorm[i],0)
                mul!(dr,rotatematrix,tempvec)
                positions[:,iatom] .+= dr
            end
        end
        return positions,numatoms
    end
    =#

    function get_positions(trj::Trajectory)
        positions = zeros(Float64,3,trj.numatoms,length(trj))
        dr = zeros(Float64,3)
        tempvec = zeros(Float64,3)
        dur = zeros(Float64,3)
        unorm = zeros(Float64,3)
        for i=1:3
            unorm[i] = norm(trj.unitcell[i])
        end

        for itrj=1:length(trj)
            for iatom=1:trj.numatoms
                #println(itrj,"\t",iatom)
                positions[1,iatom,itrj] = trj[itrj][iatom].x
                positions[2,iatom,itrj] = trj[itrj][iatom].y
                positions[3,iatom,itrj] = trj[itrj][iatom].z
                if itrj > 1
                    dx = positions[1,iatom,itrj]-positions[1,iatom,itrj-1]
                    dy = positions[2,iatom,itrj]-positions[2,iatom,itrj-1]
                    dz = positions[3,iatom,itrj]-positions[3,iatom,itrj-1]

                    dr[1] = dx
                    dr[2] = dy
                    dr[3] = dz
                    mul!(dur,trj.rotatematrix_inv,dr)

                    for i=1:3
                        tempvec .= 0
                        tempvec[i] = ifelse(dur[i] > unorm[i]/3,-unorm[i],0)
                        mul!(dr,trj.rotatematrix,tempvec)
                        positions[:,iatom,itrj] .+= dr

                        tempvec .= 0
                        tempvec[i] = ifelse(dur[i] < -unorm[i]/3,unorm[i],0)
                        mul!(dr,trj.rotatematrix,tempvec)
                        positions[:,iatom,itrj] .+= dr
                    end

                end
            end
        end
        return positions
    end

    struct Trjfile
        inputdirname::String
        filelist::Array{String,1}
        numdata_list::Array{Int64,1}
        index_list::Array{Int64,1}
        numatoms_list::Array{Int64,1}

        function Trjfile(inputdirname)
            filelist = filter(x -> x[end-3:end] == ".dat",readdir(inputdirname))
            println("filelist: ",filelist)
            numdata_list = Int64[]
            index_list = Int64[]
            numatoms_list = Int64[]
            count = 0
            for file in filelist
                println(file)
                fp = open(inputdirname*"/"*file,"r")
                #println(u)
                push!(numdata_list,read(fp,Int64))
                push!(index_list,read(fp,Int64))
                push!(numatoms_list,read(fp,Int64))
                seek(fp,0)
                close(fp)
            end

            return new(inputdirname,filelist,numdata_list,index_list,numatoms_list)
        end
    end

    function get_lastindex(trjfile::Trjfile)
        return trjfile.index_list[end] + trjfile.numdata_list[end] -1
    end

    function get_structure(trjfile::Trjfile,itrj)
        #fileindex = searchsortedfirst(trjfile.index_list,itrj)
        #println(filter(x -> x <= itrj,trjfile.index_list)[end])
        istart = filter(x -> x <= itrj,trjfile.index_list)[end]
        fileindex = searchsortedfirst(trjfile.index_list,istart)
        fp = open(trjfile.inputdirname*"/"*trjfile.filelist[fileindex],"r")

        seekposition = 0
        count = 0
        seekposition += nbytes_Int64
        #println(istart)
        for i=istart:itrj-1
            count += 1
            seekposition += 2*nbytes_Int64
            for iatom=1:trjfile.numatoms_list[fileindex]
                seekposition += 3*nbytes_Float64
            end
        end
        seek(fp,seekposition)
        index = read(fp,Int64)
        numatoms = read(fp,Int64)
        @assert numatoms == trjfile.numatoms_list[fileindex] "data structure might be wrong"
        @assert index == itrj
        #println("\t $index $numatoms")

        structure = Array{Atom,1}(undef,numatoms)

        #structure = zeros(Float64,3,numatoms)
        for iatom=1:numatoms
            x = read(fp,Float64)
            y = read(fp,Float64)
            z = read(fp,Float64)
            structure[iatom] = Atom(Float64[x,y,z])
        end
        seek(fp,0)
        close(fp)
        #println(fileindex)

        positions = get_relative_positions(Snapshot(structure))

        for iatom=1:numatoms
            structure[iatom] = Atom(positions[:,iatom])
        end

        return Snapshot(structure),fileindex

    end

    function calc_MSD_fromfile(inputdirname,dirname,itrj_start,
                        period,start_itrj,maxsteps)

        trjfile = Trjfile(inputdirname)
        println(trjfile.numdata_list)
        println(trjfile.index_list)
        println(trjfile.numatoms_list)

        snapshot,fileindex = get_structure(trjfile,itrj_start)
        snapshot_old = deepcopy(snapshot)

        numatoms = trjfile.numatoms_list[fileindex]

        count = 0

        positions = zeros(Float64,3,numatoms,maxsteps)
        positions_old = zeros(Float64,3,numatoms)

        lastindex = get_lastindex(trjfile)

        msd = zeros(maxsteps)
        phi = zeros(numatoms,maxsteps)
        ri = zeros(3)
        drj = zeros(3)
        dphi = 0

        count = 0
        for itrj=itrj_start:period:lastindex

            if itrj+maxsteps-1 > lastindex
                break
            end
            count += 1
            snapshot,fileindex = get_structure(trjfile,itrj)

            for it=2:maxsteps
                jtrj = itrj+it-1
                snapshot_j,fileindex_j = get_structure(trjfile,jtrj)

                for iatom=1:numatoms
                    dphi_x = (snapshot_j[iatom].x*a0-snapshot[iatom].x*a0)^2
                    dphi_y = (snapshot_j[iatom].y*a0-snapshot[iatom].y*a0)^2
                    dphi_z = (snapshot_j[iatom].z*a0-snapshot[iatom].z*a0)^2
                    dphi = dphi_x+dphi_y+dphi_z
                    msd[it] += dphi
                    phi[iatom,it] += dphi
                end
            end

            #println(msd[2],"\t",count)

            fp = open(dirname*"/atomic_MSD_from$(itrj_start)-th_file.txt","w")
            fp2 = open(dirname*"/MSD_from$(itrj_start)-th_file.txt","w")

            for it = 1:maxsteps
                println(fp2,it,"\t",msd[it]/(count*numatoms))
                print(fp,it,"\t")
                for iatom=1:numatoms
                    print(fp,phi[iatom,it]/count,"\t")
                end
                println(fp,"\t")
            end
            close(fp)

            close(fp2)
            println("itrj = ",itrj)

        end

    end

    function calc_MSD(trj::Trajectory,dirname,period,start_itrj,maxsteps)
        positions = get_relative_positions(trj)
        calc_MSD(trj,dirname,positions,period,start_itrj,maxsteps)
    end

    function calc_MSD(trj::Trajectory,dirname,positions,period,start_itrj,maxsteps)
        msd = zeros(maxsteps)
        phi = zeros(trj.numatoms,maxsteps)
        ri = zeros(3)
        drj = zeros(3)
        dphi = 0
        count = 0
        istart = start_itrj - trj.itrj_start + 1
        @assert istart > 0 "no data available!"
        for itrj=istart:period:length(trj)
            if itrj+maxsteps-1 > length(trj)
                break
            end
            count += 1

            for iatom=1:trj.numatoms
                for i=1:3
                    ri[i] = positions[i,iatom,itrj]*a0
                end
                for it = 1:maxsteps
                    jtrj = itrj+it-1

                    dphi = 0.0
                    for i=1:3
                        dphi += (positions[i,iatom,jtrj]*a0 - ri[i])^2
                        #drj[i] = positions[i,iatom,jtrj]*a0 - ri[i]
                    end
                    msd[it] += dphi
                    phi[iatom,it] += dphi
                end

            end
            #println(msd[2],"\t",count)
            fp = open(dirname*"/atomic_MSD_from$(start_itrj)-th.txt","w")
            fp2 = open(dirname*"/MSD_from$(start_itrj)-th.txt","w")

            for it = 1:maxsteps

                #println(positions[:,:,istart+1])
                #return

                println(fp2,it,"\t",msd[it]/(count*trj.numatoms))
                print(fp,it,"\t")
                for iatom=1:trj.numatoms
                    print(fp,phi[iatom,it]/count,"\t")
                end
                println(fp,"\t")
            end
            close(fp)

            close(fp2)
            println("itrj = ",itrj)
        end
        return msd,phi
    end

# Write your package code here.

end


text:
using BinaryBuilder

# Collection of sources required to build Xsum
sources = [
    "https://gitlab.com/radfordneal/xsum.git" =>
    "63bd4c13e4a8da4932ded362e2921ed99079f432", # version 2018-12-02 "1.0"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/xsum
mkdir -p ${prefix}/lib
if [[ $target == i686-* ]]; then
    xsumfpmath="-mfpmath=sse"
fi
${CC} ${LDFLAGS} -shared -fPIC -O3 -std=c99 ${xsumfpmath} -fno-inline-functions -o ${prefix}/lib/libxsum.${dlext} xsum.c pbinary.c
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms() # build on all supported platforms

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libxsum", :libxsum),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "xsum", v"1.0", sources, script, platforms, products, dependencies)

println("Contents of ", pwd(), " = ", readdir("."))
if isdir("products")
    println("PRODUCTS = ", readdir("products"))
else
    println("products/ directory not found")
end


text:
__precompile__()

using Base.Filesystem
import FileIO
using Images
using Logging
using Random

export
    run_and_get_output,
    run_tesseract

# Tesseract settings
command = "tesseract"
psm_valid_range = 0:14
oem_valid_range = 0:4

function check_tesseract_installed()
    try
        read(`$command --version`, String);
        @info "Tesseract is properly installed!"
    catch
        @error "Tesseract is not properly installed!"
    end
end

check_tesseract_installed()

# TODO: function to set command for executing tesseract?

function get_tesseract_version()
    info = read(`$command --version`, String)
    version = split(info, "\n")[1]
    return version
end

@info "Using Tesseract version: $(get_tesseract_version())"

function image_to_boxes(image; lang="eng", config="", nice=0, boxes=false)
    nothing
end

function image_to_osd(image; lang="eng", config="", nice=0, boxes=false)
    nothing
end

"""
Wrapper function to run Tesseract for a image stored in disk, and write the results in a given path.

# Arguments
- `input_path::String`: Path to the image to be OCRed
- `output_path::String`: Path to the text result to be stored
- `psm::Int`: Page segmentation modes (PSM):
    0.    Orientation and script detection (OSD) only.
    1.    Automatic page segmentation with OSD.
    2.    Automatic page segmentation, but no OSD, or OCR.
    3.    Fully automatic page segmentation, but no OSD. (Default)
    4.    Assume a single column of text of variable sizes.
    5.    Assume a single uniform block of vertically aligned text.
    6.    Assume a single uniform block of text.
    7.    Treat the image as a single text line.
    8.    Treat the image as a single word.
    9.    Treat the image as a single word in a circle.
    10.    Treat the image as a single character.
    11.    Sparse text. Find as much text as possible in no particular order.
    12.    Sparse text with OSD.
    13.    Raw line. Treat the image as a single text line,
        bypassing hacks that are Tesseract-specific.
- `oem::Int`: OCR Engine modes (OEM):
    0.    Legacy engine only.
    1.    Neural nets LSTM engine only. (Default)
    2.    Legacy + LSTM engines.
    3.    Default, based on what is available.

- `lang::String` Language to be configured in Tesseract (can be 'nothing').

# Examples
```julia-repl
 julia> using OCReract;
 julia> img_path = "/path/to/img.png";
 julia> out_path = "/tmp/tesseract_result";
 julia> run_tesseract(img_path, out_path, psm=3, oem=1)
```

"""
function run_tesseract(input_path::String, output_path::String; lang=nothing, psm=3, oem=1, nice=0)
    # TODO: allow to handle user-words, user-patterns
    if !isfile(input_path)
        @error "Input path '$input_path' doesn't exist!"
        return false
    end

    output_file = output_path*".txt"

    if isfile(output_file)
        @warn "Output path '$output_file' already exists!"
    end

    cmd = "tesseract $input_path $output_path"

    if lang != nothing
        cmd = join([cmd, "-l $lang"], " ")
    end

    if ! (psm in psm_valid_range)
        psm = 3
        @warn "PSM parameter not in valid range: [$(string(psm_valid_range))]. Changing to default value: PSM=$psm"
    end

    if ! (oem in oem_valid_range)
        oem = 1
        @warn "OEM parameter not in valid range: [$(string(oem_valid_range))]. Changing to default value: OEM=$oem"
    end

    cmd = join([cmd, "--oem $oem --psm $psm"], " ")

    @debug "Running command '$cmd' ..."

    try
        run(`$(split(cmd))`)
    catch e
        @info "Error ocurred while running Tesseract! $e"
    end

    return true

end

function get_tmp_path(;extension=".png")

    basepath = "/tmp/"

    path = basepath * randstring(10) * extension

    while isfile(path)
        path = basepath * randstring(10) * extension
    end

    return path
end

"""
Function to run Tesseract over an image in memory, and get the results in a String.

# Arguments
- `ìmage::Image`: Image to be OCR-ed, in Image format.
- `psm::Int`: Page segmentation modes (PSM):
    0.    Orientation and script detection (OSD) only.
    1.    Automatic page segmentation with OSD.
    2.    Automatic page segmentation, but no OSD, or OCR.
    3.    Fully automatic page segmentation, but no OSD. (Default)
    4.    Assume a single column of text of variable sizes.
    5.    Assume a single uniform block of vertically aligned text.
    6.    Assume a single uniform block of text.
    7.    Treat the image as a single text line.
    8.    Treat the image as a single word.
    9.    Treat the image as a single word in a circle.
    10.    Treat the image as a single character.
    11.    Sparse text. Find as much text as possible in no particular order.
    12.    Sparse text with OSD.
    13.    Raw line. Treat the image as a single text line,
        bypassing hacks that are Tesseract-specific.
- `oem::Int`: OCR Engine modes (OEM):
    0.    Legacy engine only.
    1.    Neural nets LSTM engine only. (Default)
    2.    Legacy + LSTM engines.
    3.    Default, based on what is available.

- `lang::String` Language to be configured in Tesseract (can be 'nothing').

# Examples

```jldoctest

julia> using Images;
julia> using OCReract;
julia> img_path = "/path/to/img.png";
julia> img = Images.load(img_path);
julia> res_text = run_and_get_output(img, psm=3, oem=1);
julia> println(res_text);
```
"""
function run_and_get_output(image; psm=3, oem=1, lang=nothing)

    # TODO: handle image type
    input_path = get_tmp_path(extension=".png")  # PNG image resulting
    output_path = get_tmp_path(extension="")  # No extension since we need to set the basename

    # Save image into disk to be handled by Tesseract
    FileIO.save(input_path, image)

    # Run Tesseract!
    run_tesseract(input_path, output_path, lang=lang, psm=psm, oem=oem)

    # Read output into a string
    output_filename = output_path*".txt"
    txt = ""

    open(output_filename, "r") do f
        txt = read(f, String)
    end

    # Now remove these temporary files generated
    try
        @debug "Removing temporary files generated..."
        rm(input_path)
        rm(output_filename)
    catch e
        @warn "Error ocurred while removing temporary files! $e"
    end

    return txt
end


text:
using DataDeps
using Test

@testset "Manual included" begin
    register(ManualDataDep(
        "Example",
        """
        This manual datadep should have been installed with the package.
        If you are seeing this message something has gone wrong.
        Try removing and then adding the package back again.
        If that doesn't work raise an issue on the repo.
        """
    ))

    content = read(datadep"Example"*"/loremipsum.txt", String)
    @test startswith(content, "lorem ipsum dolor sit amet")
end

if DataDeps.env_bool("DATADEPS_ENABLE_MANUAL_TESTS")
    @testset "Manual nonincluded" begin
        register(ManualDataDep(
            "Kafta",
            """
            Please go to
                https://www.gutenberg.org/ebooks/5200
                and download the Plain Text version as `mm.txt`
                Note this must be done manually as Project Gutenberg does not allow directly linking to their ebooks.
                only to the canonical catalog entry (See https://www.gutenberg.org/wiki/Gutenberg:Information_About_Linking_to_our_Pages)
            """
        ))

        @test first(readlines(datadep"Kafta"*"/mm.txt")) == "\ufeffThe Project Gutenberg EBook of Metamorphosis, by Franz Kafka"
    end
end


text:
"""utilities of solvers"""

"""
    SearchNode{State<:AbsState}

used in space-timed A*
"""
@kwdef struct SearchNode{State<:AbsState}
    v::Node{State}  # corresponding vertex in roadmap
    t::Int64  # time
    id::String = @sprintf("%d-%d", v.id, t)  # id of search node
    parent_id::Union{Nothing,String} = nothing  # parent in the search tree
    g::Float64 = 0.0  # g-value
    h::Float64 = 0.0  # h-value
    f::Float64 = g + h  # f-value
    unique_id::Int64 = 1   # to avoid duplication
end

"""
    gen_g_func(;
        greedy::Bool = false,           # true -> returns constant function that always returns 0
        stay_penalty::Float64 = 0.0     # adds cost for "stay" actions
    )::Function

generate cost function based on distance (g-value)
"""
function gen_g_func(; greedy::Bool = false, stay_penalty::Float64 = 0.0)::Function
    f(Q_from::Vector, Q_to::Vector)::Float64 = begin
        greedy ? 0 : dist(Q_from, Q_to)
    end

    f(q_from, q_to, i::Int64 = 1)::Float64 = begin
        greedy ? 0 : (dist(q_from, q_to) + (q_from == q_to ? stay_penalty : 0))
    end

    return f
end

"""
    find_timed_path(
        roadmap::Vector{Node{State}},   # roadmap
        invalid::Function,              # collision check function
        check_goal::Function,           # goal check function
        h_func::Function,               # heuristic function
        g_func::Function;               # cost function
        TIME_LIMIT::Union{Nothing,Float64} = nothing,  # time limit (sec)
    )::Union{
        Nothing,                        # case not found
        Vector{Node{State}}             # path
    } where {State<:AbsState}

perform space-time A*
"""
function find_timed_path(
    roadmap::Vector{Node{State}},
    invalid::Function,
    check_goal::Function,
    h_func::Function,
    g_func::Function;
    TIME_LIMIT::Union{Nothing,Float64} = nothing,
)::Union{Nothing,Vector{Node{State}}} where {State<:AbsState}

    # timeout
    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    CLOSE = Dict{String,SearchNode{State}}()
    OPEN = PriorityQueue{SearchNode{State},Float64}()

    num_generated_nodes = 1

    # setup initial node
    v_init = roadmap[1]
    S_init =
        SearchNode(v = v_init, t = 1, h = h_func(v_init), unique_id = num_generated_nodes)
    # case: not connected
    S_init.h == typemax(Float64) && return nothing
    # insert initial node
    enqueue!(OPEN, S_init, S_init.f)

    # main loop
    while !isempty(OPEN) && !timeover()

        # pop
        S = dequeue!(OPEN)
        haskey(CLOSE, S.id) && continue
        CLOSE[S.id] = S

        # check goal condition
        if check_goal(S)
            path = []
            while S.parent_id != nothing
                pushfirst!(path, S.v)
                S = CLOSE[S.parent_id]
            end
            pushfirst!(path, v_init)
            return path
        end

        # expand
        for k in vcat(S.v.neighbors, S.v.id)
            u = roadmap[k]
            num_generated_nodes += 1
            S_new = SearchNode(
                v = u,
                t = S.t + 1,
                parent_id = S.id,
                g = g_func(S, u.q),
                h = h_func(u),
                unique_id = num_generated_nodes,
            )
            (haskey(CLOSE, S_new.id) || invalid(S, S_new)) && continue
            enqueue!(OPEN, S_new, S_new.f)
        end
    end

    return nothing
end

"""
    convert_paths_to_configurations(
        paths::Vector{Vector{Node{State}}},
    )::Vector{Vector{Node{State}}} where {State<:AbsState}

transform paths (vector of paths) to configurations (vector of configurations)
"""
function convert_paths_to_configurations(
    paths::Vector{Vector{Node{State}}},
)::Vector{Vector{Node{State}}} where {State<:AbsState}

    max_len = maximum(map(e -> length(e), paths))
    return [map(path -> path[min(t, length(path))], paths) for t = 1:max_len]
end

"""

    get_distance_table(
        roadmap::Vector{Node{State}};
        g_func::Function = dist,
        goal_node::Node{State} = roadmap[2],
    )::Vector{Float64} where {State<:AbsState}

given one roadamp, compute distance table by Dijkstra method
"""
function get_distance_table(
    roadmap::Vector{Node{State}};
    g_func::Function = dist,
    goal_node::Node{State} = roadmap[2],
)::Vector{Float64} where {State<:AbsState}

    table = fill(typemax(Float64), length(roadmap))
    OPEN = PriorityQueue{Int64,Float64}()

    # setup initial vertex
    table[goal_node.id] = 0
    enqueue!(OPEN, goal_node.id, 0)

    while !isempty(OPEN)
        # pop
        v_id = dequeue!(OPEN)
        v = roadmap[v_id]
        d = table[v_id]

        # expand
        for u_id in v.neighbors
            u = roadmap[u_id]
            g = g_func(v.q, u.q) + d
            # update distance
            if g < table[u_id]
                haskey(OPEN, u_id) && delete!(OPEN, u_id)
                table[u_id] = g
                enqueue!(OPEN, u_id, g)
            end
        end
    end

    return table
end

"""dubins case requires neighbor computation"""
function get_distance_table(
    roadmap::Vector{Node{MRMP.StateDubins}};
    g_func::Function = dist,
    goal_node::Node{MRMP.StateDubins} = roadmap[2],
)::Vector{Float64}

    table = fill(typemax(Float64), length(roadmap))
    OPEN = PriorityQueue{Int64,Float64}()

    # setup initial vertex
    table[goal_node.id] = 0
    enqueue!(OPEN, goal_node.id, 0)

    while !isempty(OPEN)
        # pop
        v_id = dequeue!(OPEN)
        v = roadmap[v_id]
        d = table[v_id]

        # expand
        for u in filter(w -> v_id in w.neighbors, roadmap)
            g = g_func(v.q, u.q) + d
            # update distance
            if g < table[u.id]
                haskey(OPEN, u.id) && delete!(OPEN, u.id)
                table[u.id] = g
                enqueue!(OPEN, u.id, g)
            end
        end
    end

    return table
end

"""
    get_distance_tables(
        roadmaps::Vector{Vector{Node{State}}};
        g_func::Function = dist,
        goal_nodes::Vector{Node{State}} = map(rmp -> rmp[2], roadmaps),
    )::Vector{Vector{Float64}} where {State<:AbsState}

compute distance tables for multiple roadmaps by Dijkstra method
"""
function get_distance_tables(
    roadmaps::Vector{Vector{Node{State}}};
    g_func::Function = dist,
    goal_nodes::Vector{Node{State}} = map(rmp -> rmp[2], roadmaps),
)::Vector{Vector{Float64}} where {State<:AbsState}

    return map(
        i -> get_distance_table(roadmaps[i]; g_func = g_func, goal_node = goal_nodes[i]),
        1:length(roadmaps),
    )
end


text:
using ClusterTrees
using StaticArrays
using Test
const P = SVector{3,Float64}

using DelimitedFiles
Q = readdlm(joinpath(@__DIR__,"points.dlm"), Float64)
points = [SVector{3,Float64}(Q[i,:]) for i in axes(Q,1)]

root_center = P(0.5,0.5,0.5)
root_size = 0.5
root_node = ClusterTrees.LevelledTrees.HNode(ClusterTrees.PointerBasedTrees.Node(ClusterTrees.LevelledTrees.Data(0,Int[]), 0, 0, 0, 0), 0)
tree = ClusterTrees.LevelledTrees.LevelledTree([root_node], 1, root_center, root_size, Int[1])

smallest_box_size = 0.1
root_sector = 2
root_sfc_state = 2
root_depth = 1
for i in 1:length(points)
    target = ClusterTrees.LevelledTrees.Router(smallest_box_size, points[i])
    root_state = root(tree), root_sector, root_center, root_size, root_sfc_state, root_depth
    ClusterTrees.update!(tree, root_state, i, target) do tree, node, data
        push!(ClusterTrees.data(tree,node).values, data)
    end
end

bt = ClusterTrees.BlockTrees.BlockTree(tree,tree)
@test ClusterTrees.LevelledTrees.numlevels(tree) == 4
# ClusterTrees.print_tree(bt)

num_levels = length(tree.levels)
num_nodes_per_level = [count(x->true,ClusterTrees.LevelledTrees.LevelIterator(tree,i)) for i in 1:num_levels]

num_nodes_in_block_tree = count(x->true,ClusterTrees.DepthFirstIterator(bt, ClusterTrees.root(bt)))
num_leaves_in_block_tree = count(x->true,ClusterTrees.leaves(bt, ClusterTrees.root(bt)))

@test num_nodes_in_block_tree == sum(m^2 for m in num_nodes_per_level)
@test num_leaves_in_block_tree == num_nodes_per_level[end]^2


text:
using ITensors
using Random

using ITensors: AbstractMPS

function fill_trivial_coefficients(ψ)
  return ψ isa AbstractMPS ? (1, ψ) : ψ
end

function inner_add(α⃗ψ⃗::Tuple{<:Number,<:MPST}...) where {MPST<:AbstractMPS}
  Nₘₚₛ = length(α⃗ψ⃗)
  α⃗ = first.(α⃗ψ⃗)
  ψ⃗ = last.(α⃗ψ⃗)
  N⃡ = (conj(α⃗[i]) * α⃗[j] * inner(ψ⃗[i], ψ⃗[j]) for i in 1:Nₘₚₛ, j in 1:Nₘₚₛ)
  return sum(N⃡)
end

inner_add(ψ⃗...) = inner_add(fill_trivial_coefficients.(ψ⃗)...)

# TODO: this is no longer needed, use randomMPS
function makeRandomMPS(sites; chi::Int=4)::MPS
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(chi, "Link,l=$n") for n in 1:(N - 1)]
  for n in 1:N
    s = sites[n]
    if n == 1
      v[n] = randomITensor(l[n], s)
    elseif n == N
      v[n] = randomITensor(l[n - 1], s)
    else
      v[n] = randomITensor(l[n - 1], l[n], s)
    end
    normalize!(v[n])
  end
  return MPS(v, 0, N + 1)
end

function makeRandomMPO(sites; chi::Int=4)::MPO
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(chi, "Link,l=$n") for n in 1:(N - 1)]
  for n in 1:N
    s = sites[n]
    if n == 1
      v[n] = ITensor(l[n], s, s')
    elseif n == N
      v[n] = ITensor(l[n - 1], s, s')
    else
      v[n] = ITensor(l[n - 1], s, s', l[n])
    end
    randn!(v[n])
    normalize!(v[n])
  end
  return MPO(v, 0, N + 1)
end

# Based on https://discourse.julialang.org/t/lapackexception-1-while-svd-but-not-svdvals/23787
function make_illconditioned_matrix(T=5000)
  t = 0:(T - 1)
  f = LinRange(0, 0.5 - 1 / length(t) / 2, length(t) ÷ 2)
  y = sin.(t)
  function check_freq(f)
    zerofreq = findfirst(iszero, f)
    zerofreq !== nothing &&
      zerofreq != 1 &&
      throw(ArgumentError("If zero frequency is included it must be the first frequency"))
    return zerofreq
  end
  function get_fourier_regressor(t, f)
    zerofreq = check_freq(f)
    N = length(t)
    Nf = length(f)
    Nreg = zerofreq === nothing ? 2Nf : 2Nf - 1
    N >= Nreg || throw(ArgumentError("Too many frequency components $Nreg > $N"))
    A = zeros(N, Nreg)
    sinoffset = Nf
    for fn in 1:Nf
      if fn == zerofreq
        sinoffset = Nf - 1
      end
      for n in 1:N
        phi = 2π * f[fn] * t[n]
        A[n, fn] = cos(phi)
        if fn != zerofreq
          A[n, fn + sinoffset] = -sin(phi)
        end
      end
    end
    return A, zerofreq
  end
  A, z = get_fourier_regressor(t, f)
  return [A y]
end


text:
"""
    primal_contour(f::Function, npoints::Int)

Return `npoints` points with equally spaced angles of the 1-sublevel set of the
homogeneous function `f(x, y)`.
"""
function primal_contour(f::Function, npoints::Int)
    x = Vector{Float64}(undef, npoints)
    y = Vector{Float64}(undef, npoints)
    for (i, α) in enumerate(range(0, stop=2π - 2π/npoints, length=npoints))
        x0 = cos(α)
        y0 = sin(α)
        r = f(x0, y0)
        # f is homogeneous so f(x0/r, y0/r) = 1
        x[i] = x0 / r
        y[i] = y0 / r
    end
    return x, y
end

"""
    dual_contour(f::Function, nhalfspaces::Int, T::Type)

Return a polytope of `nhalfspaces` halfspaces defined by normal vectors of
equally spaced angles for the polar of the 1-sublevel set of the homogeneous
function `f(x, y)`.
"""
function dual_contour(f::Function, nhalfspaces::Int, ::Type{T},
                      point::Vector{T} = [0.0, 0.0],
                      x_axis::Vector{T} = [1.0, 0.0],
                      y_axis::Vector{T} = [0.0, 1.0],
                      cone = false) where T
    h = hrep(Polyhedra.HyperPlane{T, Vector{T}}[],
             Polyhedra.HalfSpace{T, Vector{T}}[], d=length(x_axis))
    for α in range(0, stop=2π - 2π/nhalfspaces, length=nhalfspaces)
        ray = x_axis * cos(α) + y_axis * sin(α)
        λ = f(ray...)
        # We have f(ray/λ) = 1 so the halfspace is
        # (point + ray / λ) ⋅ x ≤ 1 for non-cone
        # (point + ray / λ) ⋅ x ≥ 0 for coen
        a = point + ray / λ
        intersect!(h, HalfSpace(cone ? -a : a, cone ? zero(T) : one(T)))
    end
    return polyhedron(h, Polyhedra.DefaultLibrary{T}(Polyhedra.OppositeMockOptimizer))
end

function Polyhedra.planar_contour(sphere::HyperSphere; npoints=64)
    @assert dimension(sphere) == 2
    return primal_contour((x, y) -> sqrt(x^2 + y^2), npoints)
end

function Polyhedra.planar_contour(ell::PerspectiveDualOrPolarOrNot{<:AbstractEllipsoid};
                           kws...)
    return Polyhedra.planar_contour(ellipsoid(ell); kws...)
end

function Polyhedra.planar_contour(set::Union{Ellipsoid, PolySet, ConvexPolySet,
                                             Piecewise};
                                  npoints=64)
    return primal_contour(scaling_function(set), npoints)
end

function Polyhedra.planar_contour(set::PolarOf{<:Union{Piecewise{U}, ConvexPolySet{T, B, U}}};
                                  npoints=64) where {T, B, U}
    return Polyhedra.planar_contour(dual_contour(scaling_function(Polyhedra.polar(set)),
                                                 npoints, U))
end

function Polyhedra.planar_contour(set::PerspectiveDual{T, <:Householder};
                           npoints=64) where T
    @assert dimension(set) == 2
    # z is a halfspace of the primal so a ray of the dual
    z = [1.0, 0.0, 0.0]
    h1, h2 = set.set.h
    # a is a ray of the primal so a halfspace of the dual
    a = [1, h1, h2]
    b = [h1, -1, 0]
    @assert abs(dot(a, b)) < 1e-8
    c = [h2 / (1 + h1^2), h1*h2 / (1 + h1^2), -1]
    @assert abs(dot(b, c)) < 1e-8
    @assert abs(dot(a, c)) < 1e-8
    polyhedron = dual_contour(scaling_function(set), npoints, T,
                              z, b, c, true)
    # We fix z to 1.0 and eliminate it, this is cheap for H-rep
    return Polyhedra.planar_contour(fixandeliminate(polyhedron, 1, 1.0))
end

function Polyhedra.planar_contour(t::Translation; kws...)
    @assert dimension(t) == 2
    x, y = Polyhedra.planar_contour(t.set; kws...)
    return x .+ t.c[1], y .+ t.c[2]
end

@recipe function f(set::AbstractSet; npoints=64)
    @assert dimension(set) == 2
    seriestype --> :shape
    legend --> false
    Polyhedra.planar_contour(set; npoints=npoints)
end


text:
"""
Abstract type for struct that returns real-time or simulated time
when called (functor)
"""
abstract type TimeFunc end

"""
    UTCDateTimeFuncStruct()

Functor that return real-time as DateTime (UTC) when called
"""
struct UTCDateTimeFuncStruct <: TimeFunc
    func::Function
    args

    UTCDateTimeFuncStruct() = new(now, [Dates.UTC])
end
function (timefunc::UTCDateTimeFuncStruct)()
    timefunc.func(timefunc.args...)
end
UTCDateTimeFunc = UTCDateTimeFuncStruct()

"""
    FloatTimeFuncStruct()

Functor that return real-time as Float when called
"""
struct FloatTimeFuncStruct <: TimeFunc
    func::Function

    FloatTimeFuncStruct() = new(time)
end
function (timefunc::FloatTimeFuncStruct)()
    timefunc.func()
end
FloatTimeFunc = FloatTimeFuncStruct()

"""
Default time function
"""
global _time = UTCDateTimeFunc

# Time as Float64
#_time = FloatTimeFunc

# Time as DateTime (UTC)
#_time = UTCDateTimeFunc


text:
using PooksoftOptionsKit
using PyPlot
using Dates

# assetSet -
assetSet = Set{PSAbstractAsset}()

# setup the option object -
callOptionContract = PSCallOptionContract("AAPL", Date(2020,6,25), 100.0, 1.20, 1; sense=:buy)
push!(assetSet, callOptionContract)

# compute -
result_buyer = compute_option_profit_and_loss_at_expiration(assetSet,90.0,110.0)

# plot -
profit_loss_buyer = result_buyer.value
plot(profit_loss_buyer[:,1], profit_loss_buyer[:,2], "b--")     # buyer - P/L   ]
plot(profit_loss_buyer[:,1], zeros(1000))


text:
# Small script to prepare CI testing (see https://github.com/HolyLab/HolyLabRegistry
# and https://github.com/HolyLab/ImagineInterface.jl).
if VERSION ≥ v"0.7.0"
    using Pkg
end
if VERSION < v"1.1.0"
    using LibGit2
    user_regs = joinpath(DEPOT_PATH[1],"registries")
    mkpath(user_regs)
    Base.shred!(LibGit2.CachedCredentials()) do creds
        for (reg, url) in (
            "General" => "https://github.com/JuliaRegistries/General.git",
            "EmmtRegistry" => "https://github.com/emmt/EmmtRegistry")
            path = joinpath(user_regs, reg);
            LibGit2.with(Pkg.GitTools.clone(
                url, path;
                header = "registry $reg from $(repr(url))",
                credentials = creds)) do repo end
        end
    end
else
    Pkg.Registry.add("General")
    Pkg.Registry.add(RegistrySpec(url="https://github.com/emmt/EmmtRegistry"))
end


text:
module SfSolarSystem

using Unitful, UnitfulAstro, UnitfulAngles, Measurements

using ..SfPlanetary, ..SfGeometry

# useful sources:
# - mostly wikipedia
# - https://ssd.jpl.nasa.gov/?sat_phys_par (planetary satellite physical parameters)
# - https://arxiv.org/abs/1604.06129 (pluto and charon albedo)

const sol = Star("Sol", 1u"Msun", 1u"Rsun", 5772u"K", 4.83, "G2V")

const mercury = Body(
			"Mercury", (0.330114±0.000021)*1e24u"kg", Sphere((2439.4±0.1)u"km"), 0.088, 0.142,
			Orbit(sol, 0.387098u"AU", 0.205630, 3.38, 48.331, 23.124),
			Rotation(0.346, 1407.5u"hr", 0.034))

const venus = Body(
			"Venus", 4.8675u"kg" * big(10)^24, Sphere(6051.8u"km"), 0.76, 0.689,
			Orbit(sol, 0.723332u"AU", 0.006772, 3.86, 76.680, 54.884),
			Rotation(missing, -324.025u"d", 177.36))

const earth = Body(
			"Earth", 1u"Mearth", Spheroid(1u"Rearth_e", 1u"Rearth_p"), 0.306, 0.367,
			Orbit(sol, 1.000001018u"AU", 0.0167086, 7.155, -11.26064, 114.20783), # inclination to ecliptic
			Rotation(0.3307, 0.99726968u"d", 23.4392811))
const moon = Body(
			"Moon", 7.34767309u"kg" * big(10)^22, Sphere(1737.1u"km"), 0.11, 0.12,
			Orbit(earth, 384399u"km", 0.0549,  5.145u"°"),
			Rotation(0.3929, 27.321661u"d", 24))

const mars = Body(
			"Mars", 6.4171u"kg" * big(10)^23, Spheroid(3396.2u"km", 3376.2u"km"), 0.25, 0.17,
			Orbit(sol, 227939200u"km", 0.0934, 5.65, 49.558, 286.502),
			Rotation(0.3644, 1.025957u"d", 25.19))
const phobos = Body(
			"Phobos", 1.0659u"kg" * big(10)^16, Ellipsoid(13.5u"km", 11u"km", 9u"km"), missing, 0.071±0.012,
			Orbit(mars, 9376u"km", 0.0151, 1.093u"°"),
			Rotation(missing, 0.31891023u"d", 0))
const deimos = Body(
			"Deimos", 1.4762u"kg" * big(10)^15, Ellipsoid(7.5u"km", 6.1u"km", 5.5u"km"), missing, 0.068±0.007,
			Orbit(mars, 23463.2u"km", 0.00033, .93u"°"),
			Rotation(missing, 1.263u"d", 0))

const ceres = Body(
			"Ceres", 9.3835u"kg" * big(10)^20, Sphere(469.73u"km"), missing, 0.090±0.0033,
			Orbit(sol, 2.7691651545u"AU", 0.0760090291, 10.59406704, 80.3055316, 73.5976941), # inclination to ecliptic
			Rotation(0.36, 9.07417u"hr", 4))

const jupiter = Body(
			"Jupiter", 1u"Mjup", Spheroid(1u"Rjup_e", 1u"Rjup_p"), 0.503, 0.538,
			Orbit(sol, 5.2044u"AU", 0.0489, 6.09, 100.464, 273.867),
			Rotation(0.2756, 9.925u"hr", 3.13)) # tilt to orbit

const io = Body(
			"Io", 8.931938e22u"kg", Sphere(18621.6u"km"), missing, 0.63±0.02,
			Orbit(jupiter, 421700u"km", 0.0041, 0.05), # not quite an elliptical orbit; has resonance issues
			Rotation(0.37824, 1.769137786u"d", 0))

const europa = Body(
			"Europa", 4.799844e22u"kg", Sphere(1560.8u"km"), missing, 0.67±0.03,
			Orbit(jupiter, 670900u"km", 0.009, 0.47), # not quite an elliptical orbit
			Rotation(0.346, 3.551181u"d", 0.1))

const ganymede = Body(
			"Ganymede", 1.4819e23u"kg", Sphere(2634.1u"km"), missing, 0.43±0.02,
			Orbit(jupiter, 1070400u"km", 0.0013, 0.2), # this one is elliptical
			Rotation(0.3115, 7.15455296u"d", 0.33)),

const callisto = Body(
			"Callisto", 1.075938e23u"kg", Sphere(2410.3u"km"), missing, 0.17±0.02,
			Orbit(jupiter, 1882700u"km", 0.0074, 0.192), # also actually elliptical
			Rotation(0.3549, 16.6890184u"d", 0))

const saturn = Body(
			"Saturn", 5.6834u"kg" * big(10)^26, Spheroid(60268u"km", 54364u"km"), 0.342, 0.499,
			Orbit(sol, 9.5826u"AU", 0.0565, 5.51, 113.665, 339.392),
			Rotation(0.22, 38018u"s", 26.73)) # tilt to orbit)

const titan = Body(
			"Titan", 1.3452e23u"kg", Sphere(2574.73u"km"), missing, 0.22,
			Orbit(saturn, 1221870u"km", 0.0288, 0.34854),# inclination to Saturn's equator
			Rotation(0.3414, 15.945u"d", 0))

const uranus = Body(
			"Uranus", 8.681u"kg" * big(10)^25, Spheroid(25559u"km", 24973u"km"), 0.3, 0.488,
			Orbit(sol, 19.2184u"AU", 0.046381, 6.48, 74.006, 96.998857),
			Rotation(0.23, -0.71833u"d", 97.77)) # tilt to orbit

const neptune = Body(
			"Neptune", 1.02413u"kg" * big(10)^26, Spheroid(24764u"km", 24341u"km"), 0.29, 0.422,
			Orbit(sol, 30.07u"AU", 0.008678, 6.43, 131.784, 276.336),
			Rotation(0.23, 0.6713u"d", 28.32)) # tilt to orbit

const pluto = Body(
			"Pluto", 1.303u"kg" * big(10)^22, Sphere(1188.3u"km"), 0.72±0.07, 0.575, # average geometric
			Orbit(sol, 39.482u"AU", 0.2488, 11.88, 110.299, 113.834),
			Rotation(missing, 6.3872304u"d", 122.53)) # tilt to orbit
const charon = Body(
			"Charon", 1.586u"kg" * big(10)^21, Sphere(606u"km"), 0.25±0.03 , 0.35, # average geometric
			Orbit(pluto, 19591.4u"km", 0.0002, 0.08, 223.046, missing), # SMA to centre of Pluto, not to common barycenter
			Rotation(missing, 6.3872304u"d", 0))

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
# With Robin boundary conditions

# Manufactured Robin function
rfun(x) = 1.0 + ufun(x)

# Construct the FEspace
diritags = [1,2,3,4,5,7]
fespace = H1ConformingFESpace(Float64,model,order,diritags)

# Define test and trial
V = TestFESpace(fespace)
U = TrialFESpace(fespace,ufun)

# Setup integration on Robin boundary
robintags = [6,]
rtrian = BoundaryTriangulation(model,robintags)
rquad = CellQuadrature(rtrian,degree=2)

# Object describing Robin function
rfield = CellField(rtrian,rfun)

# Integrands of the Robin BC
r(v) = inner(v,rfield)
m(v,u) = inner(v,u)

# Define Robin terms
t_ΓR = AffineFETerm(m,r,rtrian,rquad)

# Dummy term that includes jumps.
# Only for testing purposes since the shape
# functions are continuous in this example.
# Thus, including this term should not change the solution
tags = [9,]
strian = SkeletonTriangulation(model,tags)
squad = CellQuadrature(strian,degree=2)
j(v,u) = inner(jump(v),jump(u))
t_ΓS = LinearFETerm(j,strian,squad)

# Define Assembler
assem = SparseMatrixAssembler(V,U)

# Define FE problem
op = LinearFEOperator(V,U,assem,t_Ω,t_ΓN,t_ΓR,t_ΓS)

# Solve!
uh = solve(solver,op)

# Define exact solution and error
e = u - uh

# Compute errors
el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

#writevtk(trian,"trian",nref=4,cellfields=["uh"=>uh,"u"=>u,"e"=>e])

@test el2 < 1.e-8
@test eh1 < 1.e-8


text:
# Julia wrapper for header: /usr/include/scip/pub_conflict.h
# Automatically generated using Clang.jl wrap_c

function SCIPconflicthdlrComp(elem1, elem2)
    ccall((:SCIPconflicthdlrComp, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPconflicthdlrCompName(elem1, elem2)
    ccall((:SCIPconflicthdlrCompName, libscip), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), elem1, elem2)
end

function SCIPconflicthdlrGetData(conflicthdlr)
    ccall((:SCIPconflicthdlrGetData, libscip), Ptr{SCIP_CONFLICTHDLRDATA}, (Ptr{SCIP_CONFLICTHDLR},), conflicthdlr)
end

function SCIPconflicthdlrSetData(conflicthdlr, conflicthdlrdata)
    ccall((:SCIPconflicthdlrSetData, libscip), Cvoid, (Ptr{SCIP_CONFLICTHDLR}, Ptr{SCIP_CONFLICTHDLRDATA}), conflicthdlr, conflicthdlrdata)
end

function SCIPconflicthdlrGetName(conflicthdlr)
    ccall((:SCIPconflicthdlrGetName, libscip), Cstring, (Ptr{SCIP_CONFLICTHDLR},), conflicthdlr)
end

function SCIPconflicthdlrGetDesc(conflicthdlr)
    ccall((:SCIPconflicthdlrGetDesc, libscip), Cstring, (Ptr{SCIP_CONFLICTHDLR},), conflicthdlr)
end

function SCIPconflicthdlrGetPriority(conflicthdlr)
    ccall((:SCIPconflicthdlrGetPriority, libscip), Cint, (Ptr{SCIP_CONFLICTHDLR},), conflicthdlr)
end

function SCIPconflicthdlrIsInitialized(conflicthdlr)
    ccall((:SCIPconflicthdlrIsInitialized, libscip), UInt32, (Ptr{SCIP_CONFLICTHDLR},), conflicthdlr)
end

function SCIPconflicthdlrGetSetupTime(conflicthdlr)
    ccall((:SCIPconflicthdlrGetSetupTime, libscip), Cdouble, (Ptr{SCIP_CONFLICTHDLR},), conflicthdlr)
end

function SCIPconflicthdlrGetTime(conflicthdlr)
    ccall((:SCIPconflicthdlrGetTime, libscip), Cdouble, (Ptr{SCIP_CONFLICTHDLR},), conflicthdlr)
end


text:
module RoMEPlotting

using Cairo, Fontconfig

using Reexport
@reexport using Gadfly
@reexport using Colors

using Statistics, LinearAlgebra
using StatsBase
using Compose
using Dates
using DistributedFactorGraphs
using KernelDensityEstimate, KernelDensityEstimatePlotting
using IncrementalInference, RoME
using DocStringExtensions
using ApproxManifoldProducts
using TensorCast
using Requires

# import ApproxManifoldProducts: mmd! # future dependency

import Gadfly: plot
import KernelDensityEstimatePlotting: plot, drawHorDens, plotKDE
import KernelDensityEstimatePlotting: getColorsByLength
import KernelDensityEstimatePlotting: plotKDE

# assuming this is a good size for everybody
@info "For larger plots when using a browswer, run Gadfly.set_default_plot_size(30cm,20cm)"
# Gadfly.set_default_plot_size(30cm,20cm)

include("ExportAPI.jl")

# EXPERIMENTAL
const AbstractMatrix__{T} = Union{AbstractArray{T,2}, Adjoint{T,<:AbstractArray{T,2}}}

# will be overwritten if flux is present (dont make const)
PlotTypesPose2 = Union{Type{Pose2Pose2}, Type{Pose2Point2BearingRange}, Type{Pose2Point2Range}, Type{Pose2Point2Bearing}}
ExtendedPose2Pose2Types = Pose2Pose2

include("services/PlotBelief.jl")
include("SolverVisualization.jl")
include("services/PlotEllipseUtils.jl")
include("RobotViz.jl")
include("services/PlotPose.jl")
include("services/PlotTree.jl")
include("CurrentlyUnmaintained.jl")
include("PlotHexUtils.jl")
include("PlotFactors.jl")
include("PlotVariables.jl")
include("PlotFactorsReload.jl")
include("Deprecated.jl")
include("NeedsFixing.jl")

function __init__()
  @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" include("FluxSpecificFeatures.jl")
  @require ImageMagick="6218d12a-5da1-5696-b52f-db25d2ecc6d1" include("imagemagick.jl")
end

end