text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ bf40c8ca-5c8c-11eb-2e48-7b906e8e6c5e
using DrWatson

# ╔═╡ 1ef22dae-5c88-11eb-1276-231730d144f6
begin
	@quickactivate

	using Images
	using Random
	using PlutoUI
	using Statistics
	using Rotations
	using Serialization

	pd = projectdir()
	USER_SETTINGS = deserialize("$pd/tmp/user_selection.jls")

	md"""
	# Carve custom faces for your cards!
	"""
end

# ╔═╡ 9d2e6708-88ea-11eb-3852-5b09b6960506
md"""
Mage.jl Deck directory path: $(@bind DECK_DIR TextField(default="/home/dusty/Shared/PlaymatProjects/PlaymatGames/Mage.jl/EDH/decks/Vannifar's Circus"))
"""

# ╔═╡ 9d987b58-88e9-11eb-1b8b-e733e3b39391
DECK_NAME = split(DECK_DIR, '/')[end]

# ╔═╡ c0a76bea-88e9-11eb-144e-65dd27acc30b
begin
	custom_imgs_dir = "$DECK_DIR/custom_images/"
	custom_img_names = [ fn for fn in readdir(custom_imgs_dir) if occursin("png", fn) ]
	push!(custom_img_names, "Backside.png")
	[ try load("$custom_imgs_dir/$img_fn") catch end for img_fn in custom_img_names ]
end

# ╔═╡ 8823075c-88e9-11eb-3aa1-e9a9af738720
deck = deserialize("$DECK_DIR/$DECK_NAME.jls")

# ╔═╡ b978978c-5ee0-11eb-1443-33825d99e295
md"""
#### $(@bind carve_custom_card CheckBox()) Carve cards in custom images dir $(@bind custom_img_select Slider(1:length(custom_img_names), show_value=true))
"""

# ╔═╡ c4bb34e4-5ede-11eb-0871-e16410f0c048
if carve_custom_card
	md"""Shrink custom card $(@bind custom_img_ratio Slider(0:0.01:1, default=0.5, show_value=true)) by percentage or $(@bind specify_custom_size CheckBox(default=true)) specify card width $(@bind custom_width NumberField(200:10:300, default=270))"""
end

# ╔═╡ b92971e6-5d70-11eb-1713-7327a0bca668
md"""
#### $(@bind carve_deck_card CheckBox()) Carve cards in deck $(@bind deck_img_select Slider(1:length(deck[:CARD_FACES]), show_value=true))
"""

# ╔═╡ 6ff626d4-5d72-11eb-110e-49b0464e95a7
if carve_deck_card
	deck_img = deck[:CARD_FACES][deck_img_select][end][end]
	save("blot_img.png", deck_img)
	deck_img
end

# ╔═╡ 501cb384-5d73-11eb-20cb-d3b96fc930d0
begin
	custom_img = load("$custom_imgs_dir/$(custom_img_names[custom_img_select])")
	imgs = carve_deck_card ? [ deck_img ] : [ custom_img ]
end;

# ╔═╡ 4c43b0ba-5ed7-11eb-3fa3-7d9b1ba5e434
if carve_custom_card
	card_ratio = height(custom_img) / width(custom_img)
	imgs[end] = specify_custom_size ? imresize(custom_img, ceil(Int, custom_width * card_ratio), custom_width) : imresize(custom_img, ratio=custom_img_ratio)
end

# ╔═╡ 4ecd1258-5d11-11eb-1b53-71740a55b68d
md"""
##### Vertically shrink image up to : $(@bind vert_shrink_by NumberField(0:200, default=160)) pixels $(@bind shrink_greedy_v CheckBox())
"""

# ╔═╡ 72c7f7b6-5c8e-11eb-26a8-61e826cea075
md"""
##### Horizontally shrink image up to : $(@bind hor_shrink_by NumberField(0:40, default=20)) pixels $(@bind shrink_greedy_h CheckBox())
"""

# ╔═╡ 4a156ba2-5d1b-11eb-07a5-730b3f65033f
md"""
 **$(@bind export_img CheckBox()) Export card image to $DECK_DIR?**
"""

# ╔═╡ 148ab572-5d15-11eb-2a54-89eaf1d95ca4
if export_img
	save("$DECK_DIR/carved_$(randstring(5)).png", imgs[end])
	if carve_deck_card
		deck[:CARD_FACES][deck_img_select] = deck[:CARD_FACES][deck_img_select][begin] => [ imgs[end] ]
	end
end

# ╔═╡ 26ba7e48-5ecb-11eb-22f1-b5013af4ec24
md"""
**Functions**
"""

# ╔═╡ deffb654-5c8d-11eb-3138-a5c7634793bc
begin
	brightness(c::RGB) = mean((c.r, c.g, c.b))
	brightness(c::RGBA) = mean((c.r, c.g, c.b))
end

# ╔═╡ e5dad786-5c8d-11eb-1da6-e5fd55842485
convolve(img, k) = imfilter(img, reflect(k))

# ╔═╡ f334877c-5ce4-11eb-1d51-83452cede3d8
energy(∇x, ∇y) = sqrt.(∇x.^2 .+ ∇y.^2)

# ╔═╡ ed9f0b4c-5c8d-11eb-0bda-4d2fe4942325
function energy(img)
	∇y = convolve(brightness.(img), Kernel.sobel()[1])
	∇x = convolve(brightness.(img), Kernel.sobel()[2])
	energy(∇x, ∇y)
end

# ╔═╡ 88d7b806-5c88-11eb-220a-9bae6e26693b
function mark_path(img, path)
	img′ = copy(img)
	m = size(img, 2)
	for (i, j) in enumerate(path)
		for j′ in j-1:j+1
			img′[i, clamp(j′, 1, m)] = RGB(1,0,1)
		end
	end
	img′
end

# ╔═╡ 93f5643a-5ba3-11eb-3bfa-3f1e453ec23b
function remove_in_each_row_views(img, column_numbers)
	@assert size(img, 1) == length(column_numbers)
	m, n = size(img)
	local img′ = similar(img, m, n-1)

	for (i, j) in enumerate(column_numbers)
		img′[i, 1:j-1] .= @view img[i, 1:j-1]
		img′[i, j:end] .= @view img[i, j+1:end]
	end

	img′
end

# ╔═╡ a1b7250c-5ba3-11eb-0e28-2b8e62c2007a
function shrink_n(img, n, min_seam, imgs=[]; show_lightning=false)
	n==0 && return push!(imgs, img)
	e = energy(img)

	seam_energy(seam) = sum(e[i, seam[i]] for i in 1:size(img, 1))
	_, min_j = findmin(map(j->seam_energy(min_seam(e, j)), 1:size(e, 2)))
	min_seam_vec = min_seam(e, min_j)
	img′ = remove_in_each_row_views(img, min_seam_vec)

	if show_lightning
		push!(imgs, mark_path(img, min_seam_vec))
	else
		push!(imgs, img′)
	end

	shrink_n(img′, n-1, min_seam, imgs)
end

# ╔═╡ 26500512-5d00-11eb-37a8-458570a946fd
function shrink!(img, min_seam, cols=true)
	img = cols ? img : rotr90(img)
	e = energy(img)

	seam_energy(seam) = sum(e[i, seam[i]] for i in 1:size(img, 1))
	_, min_j = findmin(map(j->seam_energy(min_seam(e, j)), 1:size(e, 2)))
	min_seam_vec = min_seam(e, min_j)

	img = remove_in_each_row_views(img, min_seam_vec)
	img = cols ? img : rotl90(img)

	return img
end

# ╔═╡ a7b310bc-5ba3-11eb-38e5-3d025712d57e
function greedy_seam(energies, starting_pixel::Int)
	is = [ starting_pixel ]
	m, n = size(energies)

	for k in 2:m
		es = energies[k, clamp(is[end]-1,1,n):clamp(is[end]+1,1,n)]
		push!(is, clamp(is[end] + argmin(es) - clamp(is[end],0,2), 1, n))
	end

	return is
end

# ╔═╡ d1dd1162-5d10-11eb-3d7f-8f10be30efa7
if shrink_greedy_v
	greedy_carved_v = shrink_n(rotr90(imgs[end]), vert_shrink_by, greedy_seam, show_lightning=false)
	md"""Shrink by: $(@bind greedy_v Slider(1:vert_shrink_by, show_value=true)) Save changes to img? $(@bind save_v CheckBox())"""
end

# ╔═╡ 38bc0512-5ee4-11eb-1576-1fac86b5e112
if shrink_greedy_v
	size(rotl90(greedy_carved_v[greedy_v]))
end

# ╔═╡ e061aca2-5d12-11eb-3661-59430c9f0ce4
if shrink_greedy_v && (@isdefined greedy_carved_v)
	rotl90(greedy_carved_v[greedy_v])
end

# ╔═╡ 0d2b79d4-5d13-11eb-03c5-f90edc2114e0
if (@isdefined save_v) && save_v
	imgs[end] = deepcopy(rotl90(greedy_carved_v[greedy_v]))
	"Changes saved!"
end

# ╔═╡ 3fe076c8-5c8e-11eb-3222-ad52ddd23256
if shrink_greedy_h
	greedy_carved_h = shrink_n(imgs[end], hor_shrink_by, greedy_seam, show_lightning=false)

	md"""
	Shrink by: $(@bind greedy_h Slider(1:hor_shrink_by, show_value=true))
	Save changes to img? $(@bind save_h CheckBox())
	"""
end

# ╔═╡ e1c05af0-5ee4-11eb-153f-af5abf6bd5a9
if shrink_greedy_h
	size(greedy_carved_h[greedy_h])
end

# ╔═╡ a3c7f5f8-5c8e-11eb-3807-290228dc3ea8
if shrink_greedy_h && (@isdefined greedy_carved_h)
	greedy_carved_h[greedy_h]
end

# ╔═╡ b4d3252a-5d10-11eb-16c2-1330d81f5f7c
if (@isdefined save_h) && save_h
	imgs[end]=deepcopy(greedy_carved_h[greedy_h])
	"Changes saved!"
end

# ╔═╡ 2c06e12e-5d7d-11eb-0da4-c71de86166e8
if ((@isdefined save_v) && save_v) || ((@isdefined save_h) && save_h)
	serialize("$DECK_DIR/$DECK_NAME.jls", deck)
end

# ╔═╡ Cell order:
# ╟─bf40c8ca-5c8c-11eb-2e48-7b906e8e6c5e
# ╠═1ef22dae-5c88-11eb-1276-231730d144f6
# ╠═9d2e6708-88ea-11eb-3852-5b09b6960506
# ╠═9d987b58-88e9-11eb-1b8b-e733e3b39391
# ╟─c0a76bea-88e9-11eb-144e-65dd27acc30b
# ╟─8823075c-88e9-11eb-3aa1-e9a9af738720
# ╟─501cb384-5d73-11eb-20cb-d3b96fc930d0
# ╟─b978978c-5ee0-11eb-1443-33825d99e295
# ╟─c4bb34e4-5ede-11eb-0871-e16410f0c048
# ╟─4c43b0ba-5ed7-11eb-3fa3-7d9b1ba5e434
# ╟─b92971e6-5d70-11eb-1713-7327a0bca668
# ╟─6ff626d4-5d72-11eb-110e-49b0464e95a7
# ╟─4ecd1258-5d11-11eb-1b53-71740a55b68d
# ╟─38bc0512-5ee4-11eb-1576-1fac86b5e112
# ╟─d1dd1162-5d10-11eb-3d7f-8f10be30efa7
# ╟─e061aca2-5d12-11eb-3661-59430c9f0ce4
# ╟─0d2b79d4-5d13-11eb-03c5-f90edc2114e0
# ╟─72c7f7b6-5c8e-11eb-26a8-61e826cea075
# ╟─3fe076c8-5c8e-11eb-3222-ad52ddd23256
# ╟─e1c05af0-5ee4-11eb-153f-af5abf6bd5a9
# ╟─a3c7f5f8-5c8e-11eb-3807-290228dc3ea8
# ╟─b4d3252a-5d10-11eb-16c2-1330d81f5f7c
# ╟─4a156ba2-5d1b-11eb-07a5-730b3f65033f
# ╟─148ab572-5d15-11eb-2a54-89eaf1d95ca4
# ╟─2c06e12e-5d7d-11eb-0da4-c71de86166e8
# ╟─26ba7e48-5ecb-11eb-22f1-b5013af4ec24
# ╟─deffb654-5c8d-11eb-3138-a5c7634793bc
# ╟─e5dad786-5c8d-11eb-1da6-e5fd55842485
# ╟─f334877c-5ce4-11eb-1d51-83452cede3d8
# ╟─ed9f0b4c-5c8d-11eb-0bda-4d2fe4942325
# ╟─88d7b806-5c88-11eb-220a-9bae6e26693b
# ╟─93f5643a-5ba3-11eb-3bfa-3f1e453ec23b
# ╟─a1b7250c-5ba3-11eb-0e28-2b8e62c2007a
# ╟─26500512-5d00-11eb-37a8-458570a946fd
# ╟─a7b310bc-5ba3-11eb-38e5-3d025712d57e


text:
struct GridworldAgent
    gridworld::Gridworld
    planner::GridworldPlanner
end

function run(agent::GridworldAgent; max_steps = 1e6)
    state = init(agent.gridworld)
    num_steps = 0
    while !checkGoal(agent.gridworld, state) && num_steps < max_steps
        num_steps += 1
        render(agent.gridworld, state)
        best_action, info = act(agent.planner, state)
        updateResiduals!(agent.planner, info)
        state, cost = step(agent.gridworld, state, best_action)
    end
    if num_steps < max_steps
        println("Reached goal in ", num_steps, " steps")
    else
        println("Did not reach goal")
    end
    render(agent.gridworld, state)
end

ext: jl


text:
## No * `QUIZ_BOWL`:

"""
    question(io::IO, question_type::Symbol, q, answers...)

Write a question as a line in a tab separated file specified to `io`

* `io`: IO object specifying the file
* `question_type`: a symbol to dispatch on.
* qanswers...`: Each question type has a different structure for answers

Currently these question types are supported

* `MC`: multiple choice: Answer is tuple of pairs, e.g.: `("ans1" => true, "ans2" => false, ...)` -- exactly one is `true`
* `MA`  multiple answer: Answer is tuple of pairs, e.g.: `("ans1" => true, "ans2" => false, "ans3" => true, ...)` -- one or more is `true`
* `TF`  true false: answer is either `true` or `false`
* `ESS` Essay: no answer given
* `ORD` Ordered list.  Answer  is a tuple of options, e.g.: `("option 1", "option 2", ...)`
* `MAT` Matching. Answer is a tuple of answer-match pairs `("answer1"=>"match1", "answer2"=>"match2")`
* `FIB` fill in blank: question uses `____` for blank;  answer is a tuple of possible answers, e.g. `("ans1", "ans2", ...)`.
* `FIB_PLUS`  fill in blank plus. Question has variables marked as [var1], .... Answer is collection of pairs of the form `"var1" => ("possible", "exact", "matches")`
* `FIL`  File upload: no answer
* `NUM`: numeric answer. Answer specified as `value, [tolerance]`
* `SR` short reponse "answer" is a prompt
* `OP` Likert scale. No answer
* `JUMBLED_SENTENCE`: Question marks blanks with [var1], [var2] (variable names enclosed in []); answer is a collection of pairs of the form `"var"=>"choice"` with a special pair possible of the type `nothing=>("distractor 1", "distractor 2", ...)` for 1 or more distractors

"""
function question(io::IO, qt::Symbol, q, answers...)

    qq = create_html(q)
    print(io, string(qt))
    print(io, "\t")
    print(io, qq)
    _write_answer(io, Val(qt), answers...)
    println(io)

end

## helpers
## XXX Deprecate these
@deprecate lquestion(io::IO, qt::Symbol, q, answers...)  question(io, qt, LaTeX(q), answers...)
@deprecate mdquestion(io::IO, qt::Symbol, q, answers...) question(io, qt, LaTeX′(q), answers...)

## https://help.blackboard.com/Learn/Instructor/Tests_Pools_Surveys/Reuse_Questions/Upload_Questions
## upload question types
## https://help.blackboard.com/Learn/Instructor/Tests_Pools_Surveys/Question_Types

const MC = :MC   # multiple choice: ("ans1" => true, "ans2" => false, ...)
const MA = :MA   # multiple answer: ("ans1" => true, "ans2" => false, "ans3" => true, ...)
const TF = :TF   # true false: answer::Bool
const ESS = :ESS # Essay: no answer given
const ORD = :ORD # answer = answers, in order, in container, `("one", "two", "three")`
const MAT = :MAT # answer-match pairs, e.g. `("answer"=>"match", "answer"=>"match"...)`
const FIB = :FIB # fill in blank: user answer ∈ ("ans1", "ans2", ...)
const FIB_PLUS = :FIB_PLUS # variable list answers: variable, answer as tuples
const FIL = :FIL # File upload: no answer
const SR = :SR   # short reponse: answer = short_response
const OP = :OP   # opinion
const NUM = :NUM # numeric: answer, [tolerance]
const JUMBLED_SENTENCE = :JUMBLED_SENTENCE # choices, variables as tuples
const QUIZ_BOWL = :QUIZ_BOWL # question_words, phrases as tuple

# answers = ("ans1"=>true, "ans2"=>false, ...)
function _write_answer(io::IO, ::Val{:MC}, answers)
    for (k,v) in answers
        print(io, "\t")
        print(io, create_html(k, strip=true))
        print(io, "\t")
        print(io, v ? "correct" : "incorrect")
    end
end

# answers = ("ans1"=>true, "ans2"=>false, ...)
function _write_answer(io::IO, ::Val{:MA}, answers)
    for (k,v) in answers
        print(io, "\t")
        print(io, create_html(k, strip=true))
        print(io, "\t")
        print(io, v ? "correct" : "incorrect")
    end
end

# answer [true | false]
function _write_answer(io::IO, ::Val{:TF}, answer::Bool)
    print(io, "\t")
    print(io, string(answer))
end

# no answer
function _write_answer(io::IO, ::Val{:ESS}, answers...)
    # nothing
end

# answers = (ans1, ..., ansn)
function _write_answer(io::IO, ::Val{:ORD}, answers)
    print(io, "\t")
    print(io, join(answers, "\t"))
end

## Matching answer is answer-match pairs
function _write_answer(io::IO, ::Val{:MAT}, answers)
    for (a,m) in answers
        print(io, "\t")
        print(io, create_html(a, strip=false))
        print(io, "\t")
        print(io, create_html(m,strip=false))
    end
end

# question has ____ marking blank
# answers = (ans1, ..., ansn)
function _write_answer(io::IO, ::Val{:FIB}, answers)
    print(io, "\t")
    print(io, join(answers, "\t"))
end

# quesstin has [var1], [var2], ...
# answers = ("var1" => ("possible", "exact", matches"), "var2"=> ...)
function _write_answer(io::IO, ::Val{:FIB_PLUS}, answers)
    n = length(answers)
    for i in eachindex(answers)
        (a, choices) = answers[i]
        print(io, "\t")
        print(io, a)
        for choice in choices
            print(io,"\t")
            print(io, choice)
        end
        i < n && print(io, "\t")
    end
    out = String(take!(io))
    out
end

function _write_answer(io::IO, ::Val{:FIL}, answers...)
    # nothing
end

function _write_answer(io::IO, ::Val{:SR}, answer = "A short response...")
    print(io, "\t")
    print(io, create_html(answer, strip=true))
end

function _write_answer(io::IO, ::Val{:OP}, answers...)
    # nothing
end

# answer, [tolerance]
function _write_answer(io::IO, ::Val{:NUM},  answer, tolerance=nothing)
    print(io, "\t")
    print(io, answer)
    if !isnothing(tolerance)
        print(io, "\t")
        print(io, tolerance)
    end
end

# ("var"=>"choice", "var1"=>"choice", nothing=>(distractors,))
# no markup for "choice"
function _write_answer(io::IO, ::Val{:JUMBLED_SENTENCE},  answers)
    do_last = ()
    for (var, choice) in answers
        if var == nothing
            do_last = choice
            continue
        end
        print(io, "\t")
        print(io, choice)
        print(io, "\t")
        print(io, var)
        print(io, "\t")
    end
    for choice in do_last
        print(io, "\t")
        print(io, choice)
        print(io, "\t")
    end
end

## Student gets answers, they provide question
# function _write_answer(io::IO, ::Val{:QUIZ_BOWL}, answers)
#     for (q,ph) in answers
#         print(io, "\t")

#         print(io, q)
#     end
#     print(io, "\t") #???
#     for (q,ph) in answers
#         print(io, "\t")
#         print(io, ph)
#     end
#     out = String(take!(io))
#     out
# end


text:
struct Information
    f_optimum::Float64
    x_optimum::Array{Float64}
end

"""
    Information Structure

`Information` can be used to store the true optimum in order to stop a metaheuristic early.

Properties:

- `f_optimum` known minimum.
- `x_optimum` known minimizer.

If `Options` is provided, then `optimize` will stop when `|f(x) - f(x_optimum)| < Options.f_tol`
or `‖ x - x_optimum ‖ < Options.x_tol` (euclidean distance).

# Example

If you want an approximation to the minimum with accuracy of `1e-3` (|f(x) - f(x*)| < 1e-3),
then you may use `Information`.

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> bounds = [  -10.0 -10 -10; # lower bounds
                    10.0  10 10 ] # upper bounds
2×3 Array{Float64,2}:
 -10.0  -10.0  -10.0
  10.0   10.0   10.0

julia> information = Information(f_optimum = 0.0)
Information(0.0, Float64[])

julia> options = Options(f_tol = 1e-3)
Options(0.0, 0.001, 0.0, 0.0, 1000.0, 0.0, 0.0, 0, false, true, false, :minimize)

julia> state = optimize(f, bounds, ECA(information=information, options=options))
+=========== RESULT ==========+
| Iter.: 22
| f(x) = 0.000650243
| solution.x = [0.022811671589729583, 0.007052331140376011, -0.008951836265056107]
| f calls: 474
| Total time: 0.0106 s
+============================+
```
"""
function Information(;#
    f_optimum = NaN,
    x_optimum::Array{Float64} = Float64[],
)

    Information(Float64(f_optimum), x_optimum)

end


text:
"""
Get the curve data from MSH file by breaking triangles and test into
segments.
"""
function read_curve_msh(path::AbstractString)
    # First of all this should only work with gmsh meshes
    last(splitext(path)) == ".msh" || throw("$(path) is not an XML file")
    # For now we will only handle the following tags, others are ignored
    # FIXME: some of the intersting unhandled ones are periodic ...
    tags = ("\$MeshFormat", "\$Nodes", "\$Elements")
    # We shall always store nodes, i.e. all the coordinates used to define
    # entities of the mesh
    nodes = Vector{SVector{3, Float64}}()

    entities = Dict{Int64, Vector}()

    stream = open(path)
    # Node are always 2d but we take 2d as a z = 0 plane
    is_2d = true
    used_tag = ""
    for ln in eachline(stream)
        line = split(strip(ln), " ")
        # This is a tag. Open tags get to set used_tag, EndTag is just skipped
        if startswith(first(line), "\$")
            tag = first(line)
            if in(tag, tags)
                used_tag = tag
            end
            continue
        end

        # We don't do anatyhing with the line that follows line with meshformat
        if used_tag == "\$MeshFormat"
            continue
        end

        # Id Nodes tag is used
        if used_tag == "\$Nodes"
            # The array for vertices is allocated and then we move on
            if isempty(nodes)
                num_nodes = parse(Int64, first(line))
                # NOTE: GMSH stores xyz even in 2D
                resize!(nodes, num_nodes)
                continue
            end
            # The remining line fill in the coordinates
            node_index = parse(Int64, first(line))
            # Update z plane
            z = last(line)
            is_2d = (z == "0") && is_2d
            # One based indexing
            nodes[node_index] = SVector(1+parse(Float64, line[2]),
                                        1+parse(Float64, line[3]),
                                        1+parse(Float64, z))
        end

        # Lines following Element
        if used_tag == "\$Elements"
            # Are skipped with the element is point or number of element info
            if length(line) == 1 || parse(Int64, line[2]) == 15
                continue
            end

            # The rest builds gmsh elm-type
            row = map(s -> parse(Int64, s), line)
            index, elm_type, tag = row[1:3]
            elm_nodes = tuple(row[3+tag+1:end]...)

            if haskey(entities, elm_type)
                push!(entities[elm_type], elm_nodes)
            else
                entities[elm_type] = [elm_nodes]
            end
        end
    end
    close(stream)

    # At this point we have cells. Only look at those with largest dim
    celltype = maximum(keys(entities))
    # Assert simplices tris or tets
    @assert celltype == 2 || celltype == 4

    # Break cells into segments
    # Cell local index numbering
    range = (celltype == 2) ? (1, 2, 3) : (1, 2, 3, 4)
    edge_indices = collect(Combinatorics.combinations(range, 2))

    # Now we break each cell into segments
    segments = Set{Tuple{Int, Int}}()
    for cell in entities[celltype]
        # Break cell to edges
        for (i0, i1) in edge_indices
            # One based indexing
            v0, v1 = cell[i0], cell[i1]
            push!(segments, (v0 < v1) ? (v0, v1) : (v1, v0))
        end
    end
    # Unique can be converted
    segments = collect(segments)
    # For tets or embedded tris we're done
    !is_2d && return (nodes, segments)

    # Project to 2d
    points2d = Vector{SVector{2, Float64}}(length(nodes))
    for (row, point) in enumerate(nodes)
        points2d[row] = SVector(point[1], point[2])
    end
    (points2d, segments)
end


text:
using PolynomialFactors
using AbstractAlgebra
P = PolynomialFactors

## Wilkinson (good out to 50 or so, then issues)

@testset "Test factor over Poly{Int}" begin

## test of factoring over Z[x]
R, x = ZZ["x"]
f = (x-1)*(x-2)*(x-3)
U = P.poly_factor(f)
@test prod(collect(keys(U))) - f == zero(x)

f = (x-1)*(x-2)^2*(x-3)^3
U = P.poly_factor(f)
@test U[x-3] == 3
@test f - prod([k^v for (k,v) in U]) == zero(f)

f = (x-2)*(3x-4)^2*(6x-7)^2
U = P.poly_factor(f)
@test U[6x-7] == 2

f = (x^5 - x - 1)
U = P.poly_factor(f)
@test U[f] == 1

f = x^4
U = P.poly_factor(f)
@test U[x] == 4

d = Dict(x-1 => 3, 2x-3=>4, x^2+x+1=>3)
f = prod([k^v for (k,v) in  d])
U = P.poly_factor(f)
for (k,val) in d
   @test U[k] == val
end

end

@testset "Test factor over Poly{BigInt}" begin
## BigInt
## issue #40 in Roots.jl
R,x = ZZ["x"]
f = x^2 - big(2)^256
U = P.poly_factor(f)
@test U[x - big(2)^128] == 1

p = x^15 - 1
@test length(P.poly_factor(p)) == 4

p = 1 + x^3 + x^6 + x^9 + x^12
@test length(P.poly_factor(p)) == 2
@test p - prod([k^v for (k,v) in P.poly_factor(p)]) == zero(p)

W(n) = prod([x-i for i in 1:20])
U = P.poly_factor(W(20))
@test U[x-5] == 1

## Swinnerton-Dyer Polys are slow to resolve, as over `p` they factor into linears, but over Z are irreducible.
## so the "fish out" step exhausts all possibilities.
S1 = x^2 - 2
U = P.poly_factor(S1)
@test U[S1] == 1

S2 = x^4 - 10*x^2 + 1
U = P.poly_factor(S2)
@test U[S2] == 1

S3 = x^8 - 40x^6 + 352x^4 - 960x^2 + 576
U = P.poly_factor(S3)
@test length(U) == 1

## Cyclotomic polys are irreducible over Z[x] too (https://en.wikipedia.org/wiki/Cyclotomic_polynomial)
C5 = x^4 + x^3 + x^2 + x + 1
U = P.poly_factor(C5)
@test U[C5] == 1

C10 = x^4 - x^3 + x^2 -x + 1
U = P.poly_factor(C10)
@test U[C10] == 1

C15 = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
U = P.poly_factor(C15)
@test U[C15] == 1

C25 = x^20 + x^15 + x^10 + x^5 + 1
U = P.poly_factor(C25)
@test U[C25] == 1

println("Test factor over Poly{Rational{BigInt}}")
## Rational
R,x = QQ["x"]
f = -17 * (x - 1//2)^3 * (x-3//4)^4
U = P.poly_factor(f)
@test U[-1 + 2x] == 3
@test f - prod([k^v for (k,v) in P.poly_factor(f)]) == zero(f)

end

# @testset "Test rational_roots" begin
# ### Rational roots
# R,x = ZZ["x"]
# W(n) = prod([x-i for i in 1:20])
# V = rational_roots(W(20))
# @test all(V .== 1:20)

# f = (2x-3)^4 * (5x-6)^7
# V = rational_roots(f)
# @test 3//2 in V
# @test 6//5 in V

# end

@testset "Test factormod" begin
    ## factormod
    R,x = ZZ["x"]

    # factormod has elements in GF(q)
    C10 = x^4 - x^3 + x^2 -x + 1
    U = P.factormod(C10, 5)
    # (1+x)^4
    @test length(U) == 1

    C25 = x^20 + x^15 + x^10 + x^5 + 1
    U = P.factormod(C25, 5)
    # (x+4)^20
    @test length(U) == 1

    U = P.factormod(x^4 + 1, 5)
    # (x^2+2)*(x^2+3)
    @test length(U) == 2

    U = P.factormod(x^4 + 1, 2)
    # (x+1)^4
    @test length(U) == 1

    p = 5x^5 - x^4 - 1
    U = P.factormod(p, 7)
    # 5 * (x^4+3*x^3+4*x^2+3*x+4) * (x + 1)
    V = 5 * (x^4+3*x^3+4*x^2+3*x+4) * (x + 1)
    vs = PolynomialFactors.poly_coeffs(V)
    ps = PolynomialFactors.poly_coeffs(p)
    ## v-p is not zero, v-pmod 7 should be:
    @test !all(iszero(vs - ps))
    @test all(iszero.(mod.(vs-ps,7)))

end


text:
"""
    get_SOA_table(id)
    get_SOA_table(table_name)

Given the id or name of a `mort.SOA.org` table, grab it and return it as a `MortalityTable`.

!!! Remember that not all tables have been tested to work.
"""
function get_SOA_table(id::Int)
    readXTbML(joinpath(artifact"mort.soa.org","t$id.xml"))
end

function get_SOA_table(table_name::String;source_map=table_source_map)
    entry = findfirst(x-> x.name == table_name, source_map)
    readXTbML(joinpath(artifact"mort.soa.org","t$(source_map[entry].id).xml"))
end

"""
    table(id)
    table(name)

Given the id or name of a `mort.SOA.org` table, grab it and return it as a `MortalityTable`.

!!! Remember that not all tables have been tested to work.
"""
table = get_SOA_table


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

function SCIPcreateVar(scip, var, name, lb, ub, obj, vartype, initial, removable, vardelorig, vartrans, vardeltrans, varcopy, vardata)
    ccall((:SCIPcreateVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Cstring, Cdouble, Cdouble, Cdouble, SCIP_VARTYPE, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_VARDATA}), scip, var, name, lb, ub, obj, vartype, initial, removable, vardelorig, vartrans, vardeltrans, varcopy, vardata)
end

function SCIPcreateVarBasic(scip, var, name, lb, ub, obj, vartype)
    ccall((:SCIPcreateVarBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Cstring, Cdouble, Cdouble, Cdouble, SCIP_VARTYPE), scip, var, name, lb, ub, obj, vartype)
end

function SCIPwriteVarName(scip, file, var, type)
    ccall((:SCIPwriteVarName, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}, Ptr{SCIP_VAR}, UInt32), scip, file, var, type)
end

function SCIPwriteVarsList(scip, file, vars, nvars, type, delimiter)
    ccall((:SCIPwriteVarsList, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}, Ptr{Ptr{SCIP_VAR}}, Cint, UInt32, UInt8), scip, file, vars, nvars, type, delimiter)
end

function SCIPwriteVarsLinearsum(scip, file, vars, vals, nvars, type)
    ccall((:SCIPwriteVarsLinearsum, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint, UInt32), scip, file, vars, vals, nvars, type)
end

function SCIPwriteVarsPolynomial(scip, file, monomialvars, monomialexps, monomialcoefs, monomialnvars, nmonomials, type)
    ccall((:SCIPwriteVarsPolynomial, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Cint}, Cint, UInt32), scip, file, monomialvars, monomialexps, monomialcoefs, monomialnvars, nmonomials, type)
end

function SCIPparseVar(scip, var, str, initial, removable, varcopy, vardelorig, vartrans, vardeltrans, vardata, endptr, success)
    ccall((:SCIPparseVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Cstring, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_VARDATA}, Ptr{Cstring}, Ptr{UInt32}), scip, var, str, initial, removable, varcopy, vardelorig, vartrans, vardeltrans, vardata, endptr, success)
end

function SCIPparseVarName(scip, str, var, endptr)
    ccall((:SCIPparseVarName, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Ptr{SCIP_VAR}}, Ptr{Cstring}), scip, str, var, endptr)
end

function SCIPparseVarsList(scip, str, vars, nvars, varssize, requiredsize, endptr, delimiter, success)
    ccall((:SCIPparseVarsList, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Ptr{SCIP_VAR}}, Ptr{Cint}, Cint, Ptr{Cint}, Ptr{Cstring}, UInt8, Ptr{UInt32}), scip, str, vars, nvars, varssize, requiredsize, endptr, delimiter, success)
end

function SCIPparseVarsLinearsum(scip, str, vars, vals, nvars, varssize, requiredsize, endptr, success)
    ccall((:SCIPparseVarsLinearsum, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{Cint}, Cint, Ptr{Cint}, Ptr{Cstring}, Ptr{UInt32}), scip, str, vars, vals, nvars, varssize, requiredsize, endptr, success)
end

function SCIPparseVarsPolynomial(scip, str, monomialvars, monomialexps, monomialcoefs, monomialnvars, nmonomials, endptr, success)
    ccall((:SCIPparseVarsPolynomial, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Ptr{Ptr{Ptr{SCIP_VAR}}}}, Ptr{Ptr{Ptr{Cdouble}}}, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cint}}, Ptr{Cint}, Ptr{Cstring}, Ptr{UInt32}), scip, str, monomialvars, monomialexps, monomialcoefs, monomialnvars, nmonomials, endptr, success)
end

function SCIPfreeParseVarsPolynomialData(scip, monomialvars, monomialexps, monomialcoefs, monomialnvars, nmonomials)
    ccall((:SCIPfreeParseVarsPolynomialData, libscip), Cvoid, (Ptr{SCIP_}, Ptr{Ptr{Ptr{Ptr{SCIP_VAR}}}}, Ptr{Ptr{Ptr{Cdouble}}}, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cint}}, Cint), scip, monomialvars, monomialexps, monomialcoefs, monomialnvars, nmonomials)
end

function SCIPcaptureVar(scip, var)
    ccall((:SCIPcaptureVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPreleaseVar(scip, var)
    ccall((:SCIPreleaseVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}), scip, var)
end

function SCIPchgVarName(scip, var, name)
    ccall((:SCIPchgVarName, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cstring), scip, var, name)
end

function SCIPtransformVar(scip, var, transvar)
    ccall((:SCIPtransformVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{Ptr{SCIP_VAR}}), scip, var, transvar)
end

function SCIPtransformVars(scip, nvars, vars, transvars)
    ccall((:SCIPtransformVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}), scip, nvars, vars, transvars)
end

function SCIPgetTransformedVar(scip, var, transvar)
    ccall((:SCIPgetTransformedVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{Ptr{SCIP_VAR}}), scip, var, transvar)
end

function SCIPgetTransformedVars(scip, nvars, vars, transvars)
    ccall((:SCIPgetTransformedVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}), scip, nvars, vars, transvars)
end

function SCIPgetNegatedVar(scip, var, negvar)
    ccall((:SCIPgetNegatedVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{Ptr{SCIP_VAR}}), scip, var, negvar)
end

function SCIPgetNegatedVars(scip, nvars, vars, negvars)
    ccall((:SCIPgetNegatedVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}), scip, nvars, vars, negvars)
end

function SCIPgetBinvarRepresentative(scip, var, repvar, negated)
    ccall((:SCIPgetBinvarRepresentative, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{Ptr{SCIP_VAR}}, Ptr{UInt32}), scip, var, repvar, negated)
end

function SCIPgetBinvarRepresentatives(scip, nvars, vars, repvars, negated)
    ccall((:SCIPgetBinvarRepresentatives, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Ptr{SCIP_VAR}}, Ptr{UInt32}), scip, nvars, vars, repvars, negated)
end

function SCIPflattenVarAggregationGraph(scip, var)
    ccall((:SCIPflattenVarAggregationGraph, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetProbvarLinearSum(scip, vars, scalars, nvars, varssize, constant, requiredsize, mergemultiples)
    ccall((:SCIPgetProbvarLinearSum, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{Cint}, Cint, Ptr{Cdouble}, Ptr{Cint}, UInt32), scip, vars, scalars, nvars, varssize, constant, requiredsize, mergemultiples)
end

function SCIPgetProbvarSum(scip, var, scalar, constant)
    ccall((:SCIPgetProbvarSum, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{Cdouble}), scip, var, scalar, constant)
end

function SCIPgetActiveVars(scip, vars, nvars, varssize, requiredsize)
    ccall((:SCIPgetActiveVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Ptr{Cint}, Cint, Ptr{Cint}), scip, vars, nvars, varssize, requiredsize)
end

function SCIPgetVarRedcost(scip, var)
    ccall((:SCIPgetVarRedcost, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarImplRedcost(scip, var, varfixing)
    ccall((:SCIPgetVarImplRedcost, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, UInt32), scip, var, varfixing)
end

function SCIPgetVarFarkasCoef(scip, var)
    ccall((:SCIPgetVarFarkasCoef, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarLbAtIndex(scip, var, bdchgidx, after)
    ccall((:SCIPgetVarLbAtIndex, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), scip, var, bdchgidx, after)
end

function SCIPgetVarUbAtIndex(scip, var, bdchgidx, after)
    ccall((:SCIPgetVarUbAtIndex, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), scip, var, bdchgidx, after)
end

function SCIPgetVarBdAtIndex(scip, var, boundtype, bdchgidx, after)
    ccall((:SCIPgetVarBdAtIndex, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BOUNDTYPE, Ptr{SCIP_BDCHGIDX}, UInt32), scip, var, boundtype, bdchgidx, after)
end

function SCIPgetVarWasFixedAtIndex(scip, var, bdchgidx, after)
    ccall((:SCIPgetVarWasFixedAtIndex, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, UInt32), scip, var, bdchgidx, after)
end

function SCIPgetVarSol(scip, var)
    ccall((:SCIPgetVarSol, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarSols(scip, nvars, vars, vals)
    ccall((:SCIPgetVarSols, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}), scip, nvars, vars, vals)
end

function SCIPclearRelaxSolVals(scip)
    ccall((:SCIPclearRelaxSolVals, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPsetRelaxSolVal(scip, var, val)
    ccall((:SCIPsetRelaxSolVal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, val)
end

function SCIPsetRelaxSolVals(scip, nvars, vars, vals, includeslp)
    ccall((:SCIPsetRelaxSolVals, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, UInt32), scip, nvars, vars, vals, includeslp)
end

function SCIPsetRelaxSolValsSol(scip, sol, includeslp)
    ccall((:SCIPsetRelaxSolValsSol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SOL}, UInt32), scip, sol, includeslp)
end

function SCIPisRelaxSolValid(scip)
    ccall((:SCIPisRelaxSolValid, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPmarkRelaxSolValid(scip, includeslp)
    ccall((:SCIPmarkRelaxSolValid, libscip), SCIP_RETCODE, (Ptr{SCIP_}, UInt32), scip, includeslp)
end

function SCIPmarkRelaxSolInvalid(scip)
    ccall((:SCIPmarkRelaxSolInvalid, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPgetRelaxSolVal(scip, var)
    ccall((:SCIPgetRelaxSolVal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetRelaxSolObj(scip)
    ccall((:SCIPgetRelaxSolObj, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPisStrongbranchDownFirst(scip, var)
    ccall((:SCIPisStrongbranchDownFirst, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPstartStrongbranch(scip, enablepropagation)
    ccall((:SCIPstartStrongbranch, libscip), SCIP_RETCODE, (Ptr{SCIP_}, UInt32), scip, enablepropagation)
end

function SCIPendStrongbranch(scip)
    ccall((:SCIPendStrongbranch, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPgetVarStrongbranchFrac(scip, var, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
    ccall((:SCIPgetVarStrongbranchFrac, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), scip, var, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
end

function SCIPgetVarStrongbranchWithPropagation(scip, var, solval, lpobjval, itlim, maxproprounds, down, up, downvalid, upvalid, ndomredsdown, ndomredsup, downinf, upinf, downconflict, upconflict, lperror, newlbs, newubs)
    ccall((:SCIPgetVarStrongbranchWithPropagation, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}, Ptr{UInt32}, Ptr{Clonglong}, Ptr{Clonglong}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{Cdouble}, Ptr{Cdouble}), scip, var, solval, lpobjval, itlim, maxproprounds, down, up, downvalid, upvalid, ndomredsdown, ndomredsup, downinf, upinf, downconflict, upconflict, lperror, newlbs, newubs)
end

function SCIPgetVarStrongbranchInt(scip, var, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
    ccall((:SCIPgetVarStrongbranchInt, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), scip, var, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
end

function SCIPgetVarsStrongbranchesFrac(scip, vars, nvars, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
    ccall((:SCIPgetVarsStrongbranchesFrac, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), scip, vars, nvars, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
end

function SCIPgetVarsStrongbranchesInt(scip, vars, nvars, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
    ccall((:SCIPgetVarsStrongbranchesInt, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), scip, vars, nvars, itlim, down, up, downvalid, upvalid, downinf, upinf, downconflict, upconflict, lperror)
end

function SCIPgetLastStrongbranchLPSolStat(scip, branchdir)
    ccall((:SCIPgetLastStrongbranchLPSolStat, libscip), SCIP_LPSOLSTAT, (Ptr{SCIP_}, SCIP_BRANCHDIR), scip, branchdir)
end

function SCIPgetVarStrongbranchLast(scip, var, down, up, downvalid, upvalid, solval, lpobjval)
    ccall((:SCIPgetVarStrongbranchLast, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}, Ptr{UInt32}, Ptr{Cdouble}, Ptr{Cdouble}), scip, var, down, up, downvalid, upvalid, solval, lpobjval)
end

function SCIPsetVarStrongbranchData(scip, var, lpobjval, primsol, down, up, downvalid, upvalid, iter, itlim)
    ccall((:SCIPsetVarStrongbranchData, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble, Cdouble, UInt32, UInt32, Clonglong, Cint), scip, var, lpobjval, primsol, down, up, downvalid, upvalid, iter, itlim)
end

function SCIPtryStrongbranchLPSol(scip, foundsol, cutoff)
    ccall((:SCIPtryStrongbranchLPSol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{UInt32}, Ptr{UInt32}), scip, foundsol, cutoff)
end

function SCIPgetVarStrongbranchNode(scip, var)
    ccall((:SCIPgetVarStrongbranchNode, libscip), Clonglong, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarStrongbranchLPAge(scip, var)
    ccall((:SCIPgetVarStrongbranchLPAge, libscip), Clonglong, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarNStrongbranchs(scip, var)
    ccall((:SCIPgetVarNStrongbranchs, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPaddVarLocksType(scip, var, locktype, nlocksdown, nlocksup)
    ccall((:SCIPaddVarLocksType, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_LOCKTYPE, Cint, Cint), scip, var, locktype, nlocksdown, nlocksup)
end

function SCIPaddVarLocks(scip, var, nlocksdown, nlocksup)
    ccall((:SCIPaddVarLocks, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint, Cint), scip, var, nlocksdown, nlocksup)
end

function SCIPlockVarCons(scip, var, cons, lockdown, lockup)
    ccall((:SCIPlockVarCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_CONS}, UInt32, UInt32), scip, var, cons, lockdown, lockup)
end

function SCIPunlockVarCons(scip, var, cons, lockdown, lockup)
    ccall((:SCIPunlockVarCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_CONS}, UInt32, UInt32), scip, var, cons, lockdown, lockup)
end

function SCIPchgVarObj(scip, var, newobj)
    ccall((:SCIPchgVarObj, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newobj)
end

function SCIPaddVarObj(scip, var, addobj)
    ccall((:SCIPaddVarObj, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, addobj)
end

function SCIPadjustedVarLb(scip, var, lb)
    ccall((:SCIPadjustedVarLb, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, lb)
end

function SCIPadjustedVarUb(scip, var, ub)
    ccall((:SCIPadjustedVarUb, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, ub)
end

function SCIPchgVarLb(scip, var, newbound)
    ccall((:SCIPchgVarLb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newbound)
end

function SCIPchgVarUb(scip, var, newbound)
    ccall((:SCIPchgVarUb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newbound)
end

function SCIPchgVarLbNode(scip, node, var, newbound)
    ccall((:SCIPchgVarLbNode, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Ptr{SCIP_VAR}, Cdouble), scip, node, var, newbound)
end

function SCIPchgVarUbNode(scip, node, var, newbound)
    ccall((:SCIPchgVarUbNode, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Ptr{SCIP_VAR}, Cdouble), scip, node, var, newbound)
end

function SCIPchgVarLbGlobal(scip, var, newbound)
    ccall((:SCIPchgVarLbGlobal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newbound)
end

function SCIPchgVarUbGlobal(scip, var, newbound)
    ccall((:SCIPchgVarUbGlobal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, newbound)
end

function SCIPchgVarLbLazy(scip, var, lazylb)
    ccall((:SCIPchgVarLbLazy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, lazylb)
end

function SCIPchgVarUbLazy(scip, var, lazyub)
    ccall((:SCIPchgVarUbLazy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, lazyub)
end

function SCIPtightenVarLb(scip, var, newbound, force, infeasible, tightened)
    ccall((:SCIPtightenVarLb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, force, infeasible, tightened)
end

function SCIPtightenVarUb(scip, var, newbound, force, infeasible, tightened)
    ccall((:SCIPtightenVarUb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, force, infeasible, tightened)
end

function SCIPinferVarFixCons(scip, var, fixedval, infercons, inferinfo, force, infeasible, tightened)
    ccall((:SCIPinferVarFixCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{SCIP_CONS}, Cint, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, fixedval, infercons, inferinfo, force, infeasible, tightened)
end

function SCIPinferVarLbCons(scip, var, newbound, infercons, inferinfo, force, infeasible, tightened)
    ccall((:SCIPinferVarLbCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{SCIP_CONS}, Cint, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, infercons, inferinfo, force, infeasible, tightened)
end

function SCIPinferVarUbCons(scip, var, newbound, infercons, inferinfo, force, infeasible, tightened)
    ccall((:SCIPinferVarUbCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{SCIP_CONS}, Cint, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, infercons, inferinfo, force, infeasible, tightened)
end

function SCIPinferBinvarCons(scip, var, fixedval, infercons, inferinfo, infeasible, tightened)
    ccall((:SCIPinferBinvarCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, UInt32, Ptr{SCIP_CONS}, Cint, Ptr{UInt32}, Ptr{UInt32}), scip, var, fixedval, infercons, inferinfo, infeasible, tightened)
end

function SCIPinferVarFixProp(scip, var, fixedval, inferprop, inferinfo, force, infeasible, tightened)
    ccall((:SCIPinferVarFixProp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{SCIP_PROP}, Cint, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, fixedval, inferprop, inferinfo, force, infeasible, tightened)
end

function SCIPinferVarLbProp(scip, var, newbound, inferprop, inferinfo, force, infeasible, tightened)
    ccall((:SCIPinferVarLbProp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{SCIP_PROP}, Cint, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, inferprop, inferinfo, force, infeasible, tightened)
end

function SCIPinferVarUbProp(scip, var, newbound, inferprop, inferinfo, force, infeasible, tightened)
    ccall((:SCIPinferVarUbProp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{SCIP_PROP}, Cint, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, inferprop, inferinfo, force, infeasible, tightened)
end

function SCIPinferBinvarProp(scip, var, fixedval, inferprop, inferinfo, infeasible, tightened)
    ccall((:SCIPinferBinvarProp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, UInt32, Ptr{SCIP_PROP}, Cint, Ptr{UInt32}, Ptr{UInt32}), scip, var, fixedval, inferprop, inferinfo, infeasible, tightened)
end

function SCIPtightenVarLbGlobal(scip, var, newbound, force, infeasible, tightened)
    ccall((:SCIPtightenVarLbGlobal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, force, infeasible, tightened)
end

function SCIPtightenVarUbGlobal(scip, var, newbound, force, infeasible, tightened)
    ccall((:SCIPtightenVarUbGlobal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, var, newbound, force, infeasible, tightened)
end

function SCIPcomputeVarLbGlobal(scip, var)
    ccall((:SCIPcomputeVarLbGlobal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPcomputeVarUbGlobal(scip, var)
    ccall((:SCIPcomputeVarUbGlobal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPcomputeVarLbLocal(scip, var)
    ccall((:SCIPcomputeVarLbLocal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPcomputeVarUbLocal(scip, var)
    ccall((:SCIPcomputeVarUbLocal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarMultaggrLbGlobal(scip, var)
    ccall((:SCIPgetVarMultaggrLbGlobal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarMultaggrUbGlobal(scip, var)
    ccall((:SCIPgetVarMultaggrUbGlobal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarMultaggrLbLocal(scip, var)
    ccall((:SCIPgetVarMultaggrLbLocal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarMultaggrUbLocal(scip, var)
    ccall((:SCIPgetVarMultaggrUbLocal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarClosestVlb(scip, var, sol, closestvlb, closestvlbidx)
    ccall((:SCIPgetVarClosestVlb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_SOL}, Ptr{Cdouble}, Ptr{Cint}), scip, var, sol, closestvlb, closestvlbidx)
end

function SCIPgetVarClosestVub(scip, var, sol, closestvub, closestvubidx)
    ccall((:SCIPgetVarClosestVub, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_SOL}, Ptr{Cdouble}, Ptr{Cint}), scip, var, sol, closestvub, closestvubidx)
end

function SCIPaddVarVlb(scip, var, vlbvar, vlbcoef, vlbconstant, infeasible, nbdchgs)
    ccall((:SCIPaddVarVlb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Ptr{UInt32}, Ptr{Cint}), scip, var, vlbvar, vlbcoef, vlbconstant, infeasible, nbdchgs)
end

function SCIPaddVarVub(scip, var, vubvar, vubcoef, vubconstant, infeasible, nbdchgs)
    ccall((:SCIPaddVarVub, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Ptr{UInt32}, Ptr{Cint}), scip, var, vubvar, vubcoef, vubconstant, infeasible, nbdchgs)
end

function SCIPaddVarImplication(scip, var, varfixing, implvar, impltype, implbound, infeasible, nbdchgs)
    ccall((:SCIPaddVarImplication, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, UInt32, Ptr{SCIP_VAR}, SCIP_BOUNDTYPE, Cdouble, Ptr{UInt32}, Ptr{Cint}), scip, var, varfixing, implvar, impltype, implbound, infeasible, nbdchgs)
end

function SCIPaddClique(scip, vars, values, nvars, isequation, infeasible, nbdchgs)
    ccall((:SCIPaddClique, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Ptr{UInt32}, Cint, UInt32, Ptr{UInt32}, Ptr{Cint}), scip, vars, values, nvars, isequation, infeasible, nbdchgs)
end

function SCIPcalcCliquePartition(scip, vars, nvars, cliquepartition, ncliques)
    ccall((:SCIPcalcCliquePartition, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Cint}, Ptr{Cint}), scip, vars, nvars, cliquepartition, ncliques)
end

function SCIPcalcNegatedCliquePartition(scip, vars, nvars, cliquepartition, ncliques)
    ccall((:SCIPcalcNegatedCliquePartition, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Cint}, Ptr{Cint}), scip, vars, nvars, cliquepartition, ncliques)
end

function SCIPcleanupCliques(scip, infeasible)
    ccall((:SCIPcleanupCliques, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{UInt32}), scip, infeasible)
end

function SCIPgetNCliques(scip)
    ccall((:SCIPgetNCliques, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNCliquesCreated(scip)
    ccall((:SCIPgetNCliquesCreated, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetCliques(scip)
    ccall((:SCIPgetCliques, libscip), Ptr{Ptr{SCIP_CLIQUE}}, (Ptr{SCIP_},), scip)
end

function SCIPhaveVarsCommonClique(scip, var1, value1, var2, value2, regardimplics)
    ccall((:SCIPhaveVarsCommonClique, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}, UInt32, Ptr{SCIP_VAR}, UInt32, UInt32), scip, var1, value1, var2, value2, regardimplics)
end

function SCIPwriteCliqueGraph(scip, fname, writenodeweights)
    ccall((:SCIPwriteCliqueGraph, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, UInt32), scip, fname, writenodeweights)
end

function SCIPremoveVarFromGlobalStructures(scip, var)
    ccall((:SCIPremoveVarFromGlobalStructures, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPchgVarBranchFactor(scip, var, branchfactor)
    ccall((:SCIPchgVarBranchFactor, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, branchfactor)
end

function SCIPscaleVarBranchFactor(scip, var, scale)
    ccall((:SCIPscaleVarBranchFactor, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, scale)
end

function SCIPaddVarBranchFactor(scip, var, addfactor)
    ccall((:SCIPaddVarBranchFactor, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, addfactor)
end

function SCIPchgVarBranchPriority(scip, var, branchpriority)
    ccall((:SCIPchgVarBranchPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint), scip, var, branchpriority)
end

function SCIPupdateVarBranchPriority(scip, var, branchpriority)
    ccall((:SCIPupdateVarBranchPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint), scip, var, branchpriority)
end

function SCIPaddVarBranchPriority(scip, var, addpriority)
    ccall((:SCIPaddVarBranchPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint), scip, var, addpriority)
end

function SCIPchgVarBranchDirection(scip, var, branchdirection)
    ccall((:SCIPchgVarBranchDirection, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, branchdirection)
end

function SCIPchgVarType(scip, var, vartype, infeasible)
    ccall((:SCIPchgVarType, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_VARTYPE, Ptr{UInt32}), scip, var, vartype, infeasible)
end

function SCIPfixVar(scip, var, fixedval, infeasible, fixed)
    ccall((:SCIPfixVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{UInt32}, Ptr{UInt32}), scip, var, fixedval, infeasible, fixed)
end

function SCIPaggregateVars(scip, varx, vary, scalarx, scalary, rhs, infeasible, redundant, aggregated)
    ccall((:SCIPaggregateVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt32}), scip, varx, vary, scalarx, scalary, rhs, infeasible, redundant, aggregated)
end

function SCIPmultiaggregateVar(scip, var, naggvars, aggvars, scalars, constant, infeasible, aggregated)
    ccall((:SCIPmultiaggregateVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, Ptr{UInt32}, Ptr{UInt32}), scip, var, naggvars, aggvars, scalars, constant, infeasible, aggregated)
end

function SCIPdoNotAggr(scip)
    ccall((:SCIPdoNotAggr, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPdoNotMultaggr(scip)
    ccall((:SCIPdoNotMultaggr, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPdoNotMultaggrVar(scip, var)
    ccall((:SCIPdoNotMultaggrVar, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPallowDualReds(scip)
    ccall((:SCIPallowDualReds, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPallowObjProp(scip)
    ccall((:SCIPallowObjProp, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPmarkDoNotMultaggrVar(scip, var)
    ccall((:SCIPmarkDoNotMultaggrVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPenableVarHistory(scip)
    ccall((:SCIPenableVarHistory, libscip), Cvoid, (Ptr{SCIP_},), scip)
end

function SCIPdisableVarHistory(scip)
    ccall((:SCIPdisableVarHistory, libscip), Cvoid, (Ptr{SCIP_},), scip)
end

function SCIPupdateVarPseudocost(scip, var, solvaldelta, objdelta, weight)
    ccall((:SCIPupdateVarPseudocost, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble), scip, var, solvaldelta, objdelta, weight)
end

function SCIPgetVarPseudocostVal(scip, var, solvaldelta)
    ccall((:SCIPgetVarPseudocostVal, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, solvaldelta)
end

function SCIPgetVarPseudocostValCurrentRun(scip, var, solvaldelta)
    ccall((:SCIPgetVarPseudocostValCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, solvaldelta)
end

function SCIPgetVarPseudocost(scip, var, dir)
    ccall((:SCIPgetVarPseudocost, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarPseudocostCurrentRun(scip, var, dir)
    ccall((:SCIPgetVarPseudocostCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarPseudocostCount(scip, var, dir)
    ccall((:SCIPgetVarPseudocostCount, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarPseudocostCountCurrentRun(scip, var, dir)
    ccall((:SCIPgetVarPseudocostCountCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarPseudocostVariance(scip, var, dir, onlycurrentrun)
    ccall((:SCIPgetVarPseudocostVariance, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR, UInt32), scip, var, dir, onlycurrentrun)
end

function SCIPcalculatePscostConfidenceBound(scip, var, dir, onlycurrentrun, clevel)
    ccall((:SCIPcalculatePscostConfidenceBound, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR, UInt32, SCIP_CONFIDENCELEVEL), scip, var, dir, onlycurrentrun, clevel)
end

function SCIPsignificantVarPscostDifference(scip, varx, fracx, vary, fracy, dir, clevel, onesided)
    ccall((:SCIPsignificantVarPscostDifference, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{SCIP_VAR}, Cdouble, SCIP_BRANCHDIR, SCIP_CONFIDENCELEVEL, UInt32), scip, varx, fracx, vary, fracy, dir, clevel, onesided)
end

function SCIPpscostThresholdProbabilityTest(scip, var, frac, threshold, dir, clevel)
    ccall((:SCIPpscostThresholdProbabilityTest, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble, SCIP_BRANCHDIR, SCIP_CONFIDENCELEVEL), scip, var, frac, threshold, dir, clevel)
end

function SCIPisVarPscostRelerrorReliable(scip, var, threshold, clevel)
    ccall((:SCIPisVarPscostRelerrorReliable, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, SCIP_CONFIDENCELEVEL), scip, var, threshold, clevel)
end

function SCIPgetVarPseudocostScore(scip, var, solval)
    ccall((:SCIPgetVarPseudocostScore, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, solval)
end

function SCIPgetVarPseudocostScoreCurrentRun(scip, var, solval)
    ccall((:SCIPgetVarPseudocostScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, solval)
end

function SCIPgetVarVSIDS(scip, var, dir)
    ccall((:SCIPgetVarVSIDS, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarVSIDSCurrentRun(scip, var, dir)
    ccall((:SCIPgetVarVSIDSCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarConflictScore(scip, var)
    ccall((:SCIPgetVarConflictScore, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarConflictScoreCurrentRun(scip, var)
    ccall((:SCIPgetVarConflictScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarConflictlengthScore(scip, var)
    ccall((:SCIPgetVarConflictlengthScore, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarConflictlengthScoreCurrentRun(scip, var)
    ccall((:SCIPgetVarConflictlengthScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarAvgConflictlength(scip, var, dir)
    ccall((:SCIPgetVarAvgConflictlength, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarAvgConflictlengthCurrentRun(scip, var, dir)
    ccall((:SCIPgetVarAvgConflictlengthCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarAvgInferences(scip, var, dir)
    ccall((:SCIPgetVarAvgInferences, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarAvgInferencesCurrentRun(scip, var, dir)
    ccall((:SCIPgetVarAvgInferencesCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarAvgInferenceScore(scip, var)
    ccall((:SCIPgetVarAvgInferenceScore, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarAvgInferenceScoreCurrentRun(scip, var)
    ccall((:SCIPgetVarAvgInferenceScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPinitVarBranchStats(scip, var, downpscost, uppscost, downvsids, upvsids, downconflen, upconflen, downinfer, upinfer, downcutoff, upcutoff)
    ccall((:SCIPinitVarBranchStats, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), scip, var, downpscost, uppscost, downvsids, upvsids, downconflen, upconflen, downinfer, upinfer, downcutoff, upcutoff)
end

function SCIPinitVarValueBranchStats(scip, var, value, downvsids, upvsids, downconflen, upconflen, downinfer, upinfer, downcutoff, upcutoff)
    ccall((:SCIPinitVarValueBranchStats, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), scip, var, value, downvsids, upvsids, downconflen, upconflen, downinfer, upinfer, downcutoff, upcutoff)
end

function SCIPgetVarAvgCutoffs(scip, var, dir)
    ccall((:SCIPgetVarAvgCutoffs, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarAvgCutoffsCurrentRun(scip, var, dir)
    ccall((:SCIPgetVarAvgCutoffsCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR), scip, var, dir)
end

function SCIPgetVarAvgCutoffScore(scip, var)
    ccall((:SCIPgetVarAvgCutoffScore, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarAvgCutoffScoreCurrentRun(scip, var)
    ccall((:SCIPgetVarAvgCutoffScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetVarAvgInferenceCutoffScore(scip, var, cutoffweight)
    ccall((:SCIPgetVarAvgInferenceCutoffScore, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, cutoffweight)
end

function SCIPgetVarAvgInferenceCutoffScoreCurrentRun(scip, var, cutoffweight)
    ccall((:SCIPgetVarAvgInferenceCutoffScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, cutoffweight)
end

function SCIPprintVar(scip, var, file)
    ccall((:SCIPprintVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{FILE}), scip, var, file)
end


text:
mutable struct Music <: AbstractSound{sfMusic}
    ptr::Ptr{sfMusic}
    function Music(file::AbstractString)
        ptr = sfMusic_createFromFile(file)
        ptr == C_NULL && error("Couldn't create Music object from file '$file'.")
        this = new(ptr)
        finalizer(this) do obj
            try
                destroy!(obj)
            finally
                return
            end
        end
        return this
    end
end

destroy!(this::Music) = sfMusic_destroy(this)

timepos!(music::Music, offset::Integer)  = sfMusic_setPlayingOffset(music, sfTime(offset))
volume!(music::Music, volume::Real)      = sfMusic_setVolume(music, volume)
loop!(music::Music, loop::Bool)          = sfMusic_setLoop(music, loop ? sfTrue : sfFalse)
pitch!(music::Music, pitch::Real)        = sfMusic_setPitch(music, pitch)

volume(music::Music)    = sfMusic_getVolume(music)
getstatus(music::Music) = Status[ sfMusic_getStatus(music) ]
pitch(music::Music)     = sfMusic_getPitch(music)
timepos(music::Music)   = sfMusic_getPlayingOffset(music).microseconds
islooping(music::Music) = sfMusic_getLoop(music) == sfTrue

play!(music::Music)  = sfMusic_play(music)
stop!(music::Music)  = sfMusic_stop(music)
pause!(music::Music) = sfMusic_pause(music)

isplaying(music::Music) = getstatus(music) == :playing
ispaused(music::Music)  = getstatus(music) == :paused
isstopped(music::Music) = getstatus(music) == :stopped

nchannels(music::Music)  = Int(sfMusic_getChannelCount(music))
duration(music::Music)   = sfMusic_getDuration(music).microseconds
samplerate(music::Music) = Int(sfMusic_getSampleRate(music))

getprops(music::Music) = (
    loop = islooping(music),
    status = getstatus(music),
    timepos = timepos(music),
    volume = volume(music),
    pitch = pitch(music),
    channels = nchannels(music),
    duration = duration(music),
    samplerate = samplerate(music),
)

function Base.show(io::IO, music::Music)
    d = round(1e-6 * duration(music), digits = 2)
    n = nchannels(music)
    sr = round(1e-3 * samplerate(music), digits = 1)
    print(io, "Music ♫ $(d)s ($(sr)kHz, $(n)-channels)")
end


text:
# Combined 1D/2D plots across multiple frames, multi-processing version.
#
# To run on a single node,
# julia -p $ncores demo_1d2d_mp_pyplot.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Distributed, Glob
@everywhere using Vlasiator, PyPlot, Printf, LaTeXStrings

@assert matplotlib.__version__ ≥ "3.4" "Require Matplotlib version 3.4+ to use subfigure!"

@everywhere struct Varminmax{T}
   "Density, [amu/cc]"
   ρmin::T
   ρmax::T
   "Velocity, [km/s]"
   vmin::T
   vmax::T
   "Pressure, [nPa]"
   pmin::T
   pmax::T
   "Magnetic field, [nT]"
   bmin::T
   bmax::T
   "Electric field, [nT]"
   emin::T
   emax::T
end

@everywhere function init_figure(loc, norms, ticks, pArgs1, varminmax)
   fig = plt.figure(myid(), constrained_layout=true, figsize=(12, 7.2))
   subfigs = fig.subfigures(1, 2, wspace=0.01, width_ratios=[2,1])

   axsL = subfigs[1].subplots(5, 1, sharex=true)
   axsR = subfigs[2].subplots(2, 1, sharex=true)

   # Set line plots' axes
   axsL[end].set_xlim(loc[1], loc[end])

   (;ρmin, ρmax, vmin, vmax, pmin, pmax, bmin, bmax, emin, emax) = varminmax

   axsL[1].set_ylim(ρmin, ρmax)
   axsL[2].set_ylim(vmin, vmax)
   axsL[3].set_ylim(pmin, pmax)
   axsL[4].set_ylim(bmin, bmax)
   axsL[5].set_ylim(emin, emax)
   for ax in axsL
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.grid(true)
   end

   axsL[end].set_xlabel(L"x [$R_E$]"; fontsize=14)
   axsL[1].set_ylabel("n [amu/cc]";   fontsize=14)
   axsL[2].set_ylabel("V [km/s]";     fontsize=14)
   axsL[3].set_ylabel("P [nPa]";      fontsize=14)
   axsL[4].set_ylabel("B [nT]";       fontsize=14)
   axsL[5].set_ylabel("E [mV/m]";     fontsize=14)

   fakeline = loc
   l1 = axsL[1].plot(loc, fakeline, label="Proton density", color="#1f77b4")
   l2 = axsL[2].plot(loc, fakeline, label="Vx",             color="#1f77b4")
   l3 = axsL[2].plot(loc, fakeline, label="Vy",             color="#ff7f0e")
   l4 = axsL[2].plot(loc, fakeline, label="Vz",             color="#2ca02c")
   l5 = axsL[3].plot(loc, fakeline, label="Ram",            color="#1f77b4")
   l6 = axsL[3].plot(loc, fakeline, label="Thermal",        color="#ff7f0e")
   l7 = axsL[4].plot(loc, fakeline, label="Bx",             color="#1f77b4")
   l8 = axsL[4].plot(loc, fakeline, label="By",             color="#ff7f0e")
   l9 = axsL[4].plot(loc, fakeline, label="Bz",             color="#2ca02c")
   l10= axsL[5].plot(loc, fakeline, label="Ex",             color="#1f77b4")
   l11= axsL[5].plot(loc, fakeline, label="Ey",             color="#ff7f0e")
   l12= axsL[5].plot(loc, fakeline, label="Ez",             color="#2ca02c")

   ls = (l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12)

   axsL[2].legend(;loc="lower left",  ncol=3, frameon=false, fontsize=12)
   axsL[3].legend(;loc="upper right", ncol=2, frameon=false, fontsize=12)
   axsL[4].legend(;loc="upper right", ncol=3, frameon=false, fontsize=12)
   axsL[5].legend(;loc="lower right", ncol=3, frameon=false, fontsize=12)

   vl1 = axsL[1].vlines(loc[1], ρmin, ρmax; colors="r", linestyle="dashed", alpha=0.5)
   vl2 = axsL[2].vlines(loc[1], vmin, vmax; colors="r", linestyle="dashed", alpha=0.5)
   vl3 = axsL[3].vlines(loc[1], pmin, pmax; colors="r", linestyle="dashed", alpha=0.5)
   vl4 = axsL[4].vlines(loc[1], bmin, bmax; colors="r", linestyle="dashed", alpha=0.5)
   vl5 = axsL[5].vlines(loc[1], emin, emax; colors="r", linestyle="dashed", alpha=0.5)

   hl2 = axsL[2].hlines(0.0, loc[1], loc[end]; colors="k", linestyle="dashed", alpha=0.2)
   hl4 = axsL[4].hlines(0.0, loc[1], loc[end]; colors="k", linestyle="dashed", alpha=0.2)
   hl5 = axsL[5].hlines(0.0, loc[1], loc[end]; colors="k", linestyle="dashed", alpha=0.2)

   vlines = (vl1, vl2, vl3, vl4, vl5)

   for ax in axsR
      ax.set_aspect("equal")

      # Set border line widths
      for loc in ("left", "bottom", "right", "top")
         edge = get(ax.spines, loc, nothing)
         edge.set_linewidth(2.0)
      end
      ax.xaxis.set_tick_params(width=2.0, length=3)
      ax.yaxis.set_tick_params(width=2.0, length=3)
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())

      ax.set_ylabel(pArgs1.stry; fontsize=14)
   end

   axsR[2].set_xlabel(pArgs1.strx; fontsize=14)

   axsR[1].set_title("Alfvén speed", fontsize=14)
   axsR[2].set_title("Sound speed", fontsize=14)

   x, y = Vlasiator.get_axis(pArgs1)
   fakedata = fill(NaN32, length(y), length(x))

   c1 = axsR[1].pcolormesh(x, y, fakedata, norm=norms[1], cmap=matplotlib.cm.turbo)
   c2 = axsR[2].pcolormesh(x, y, fakedata, norm=norms[2], cmap=matplotlib.cm.turbo)

   rInner = 31.8e6 # [m]
   circle1 = plt.Circle((0, 0), rInner/Vlasiator.RE, facecolor="w", edgecolor="tab:purple")
   circle2 = plt.Circle((0, 0), rInner/Vlasiator.RE, facecolor="w", edgecolor="tab:purple")
   axsR[1].add_patch(circle1)
   axsR[2].add_patch(circle2)

   im_ratio = length(y)/length(x)
   fraction = 0.046 * im_ratio

   cb1 = colorbar(c1; ax=axsR[1], ticks=ticks[1], fraction, pad=0.02, extend="max")
   cb1.ax.set_ylabel("[km/s]"; fontsize=14)

   cb2 = colorbar(c2; ax=axsR[2], ticks=ticks[2], fraction, pad=0.02, extend="max")
   cb2.ax.set_ylabel("[km/s]"; fontsize=14)

   fig.suptitle("Density Pulse Run", fontsize="x-large")

   cs = (c1, c2)

   return fig, subfigs, ls, vlines, cs
end

@everywhere function update_vline(h, x)
   seg_old = h.get_segments()
   ymin = seg_old[1][1, 2]
   ymax = seg_old[1][2, 2]

   seg_new = [[x ymin; x ymax]]

   h.set_segments(seg_new)
end

@everywhere function update_plot!(subfigs, ls, vlines, cs, outdir, file, cellids, loc)
   isfile(outdir*file[end-8:end-5]*".png") && return

   println("file = $(basename(file))")
   meta = load(file)

   rho = readvariable(meta, "proton/vg_rho", cellids) |> vec
   v = readvariable(meta, "proton/vg_v", cellids)
   p = readvariable(meta, "vg_pressure", cellids) .* 1f9 |> vec # [nPa]

   vmag2 = sum(x -> x*x, v, dims=1) |> vec
   pram = rho .* Vlasiator.mᵢ .* vmag2 .* 1f9 # [nPa]

   b = readvariable(meta, "vg_b_vol", cellids) .* 1f9 #[nT]
   e = readvariable(meta, "vg_e_vol", cellids) .* 1f3 #[mV/m]

   ls[1][1].set_ydata(rho ./ 1f6)
   ls[2][1].set_ydata(@views v[1,:] ./ 1f3)
   ls[3][1].set_ydata(@views v[2,:] ./ 1f3)
   ls[4][1].set_ydata(@views v[3,:] ./ 1f3)
   ls[5][1].set_ydata(pram)
   ls[6][1].set_ydata(p)
   ls[7][1].set_ydata(@view b[1,:])
   ls[8][1].set_ydata(@view b[2,:])
   ls[9][1].set_ydata(@view b[3,:])
   ls[10][1].set_ydata(@view e[1,:])
   ls[11][1].set_ydata(@view e[2,:])
   ls[12][1].set_ydata(@view e[3,:])

   imagnetopause_ = findfirst(<(0.0), @views b[3,:])
   for vline in vlines
      update_vline(vline, loc[imagnetopause_])
   end

   str_title = @sprintf "Sun-Earth line, t= %4.1fs" meta.time
   subfigs[1].suptitle(str_title, fontsize="x-large")

   data = Vlasiator.prep2d(meta, "VA", :z)'
   cs[1].set_array(data ./ 1f3)

   data = Vlasiator.prep2d(meta, "VS", :z)'
   cs[2].set_array(data ./ 1f3)

   savefig(outdir*file[end-8:end-5]*".png", bbox_inches="tight")
   return
end

function make_jobs(files)
   for f in files
       put!(jobs, f)
   end
end

@everywhere function do_work(jobs, status,
   outdir, loc, norms, ticks, pArgs1, cellids, varminmax)

   fig, subfigs, ls, vlines, cs = init_figure(loc, norms, ticks, pArgs1, varminmax)
   while true
      file = take!(jobs)
      update_plot!(subfigs, ls, vlines, cs, outdir, file, cellids, loc)
      put!(status, true)
   end
   close(fig)
end

############################################################################################
files = glob("bulk*.vlsv", ".")

nfile = length(files)
# Set output directory
outdir = "1d2d/"

# Set contour plots' axes
axisunit = EARTH
# Upper/lower limits for each variable
ρmin, ρmax = 0.0, 14.0     # [amu/cc]
vmin, vmax = -620.0, 150.0 # [km/s]
pmin, pmax = 0.0, 2.8      # [nPa]
bmin, bmax = -60.0, 60.0   # [nT]
emin, emax = -8.0, 8.0     # [mV/m]
vamin, vamax = 50.0, 600.0 # [km/s]
vsmin, vsmax = 50.0, 600.0 # [km/s]

varminmax = Varminmax(ρmin, ρmax, vmin, vmax, pmin, pmax, bmin, bmax, emin, emax)

meta = load(files[1])

pArgs1 = Vlasiator.set_args(meta, "VA", axisunit; normal=:none)
norm1, ticks1 = Vlasiator.set_colorbar(Linear, vamin, vamax)
norm2, ticks2 = Vlasiator.set_colorbar(Linear, vsmin, vsmax)

norms = (norm1, norm2)
ticks = (ticks1, ticks2)

const jobs   = RemoteChannel(()->Channel{String}(nfile))
const status = RemoteChannel(()->Channel{Bool}(nworkers()))

xmin, xmax = 7.0, 17.0 # Earth radii
point1 = [xmin, 0, 0] .* Vlasiator.RE
point2 = [xmax, 0, 0] .* Vlasiator.RE

meta = load(files[1])
cellids, _, _ = getcellinline(meta, point1, point2)
loc = range(xmin, xmax, length=length(cellids))

println("Total number of files: $nfile")
println("Running with $(nworkers()) workers...")

@async make_jobs(files) # Feed the jobs channel with all files to process.

@sync for p in workers()
   @async remote_do(do_work, p, jobs, status,
      outdir, loc, norms, ticks, pArgs1, cellids, varminmax)
end

let n = nfile
   t = @elapsed while n > 0 # wait for all jobs to complete
      take!(status)
      n -= 1
   end
   println("Finished in $(round(t, digits=2))s.")
end


text:
using Test
using FlightMechanics

# TESTS for TAS CAS EAS based on values from
# http://www.hochwarth.com/misc/AviationCalculator.html

# --- height = 0 m ---
h = 0.0  # m
T, p, rho, a = atmosphere_isa(h)
# From TAS
tas = 100.0 * KT2MS  # m/s
expected_eas = 100.0 * KT2MS  # m/s
expected_cas = 100.0 * KT2MS  # m/s

eas = tas2eas(tas, rho)
cas = tas2cas(tas, rho, p)
@test isapprox(eas, expected_eas)
@test isapprox(cas, expected_cas)

# From CAS
cas = 100.0 * KT2MS  # m/s
expected_eas = 100.0 * KT2MS  # m/s
expected_tas = 100.0 * KT2MS  # m/s

eas = cas2eas(cas, rho, p)
tas = cas2tas(cas, rho, p)
@test isapprox(eas, expected_eas)
@test isapprox(tas, expected_tas)

# From EAS
eas = 100.0 * KT2MS  # m/s
expected_cas = 100.0 * KT2MS  # m/s
expected_tas = 100.0 * KT2MS  # m/s

cas = eas2cas(eas, rho, p)
tas = eas2tas(eas, rho)
@test isapprox(cas, expected_cas)
@test isapprox(tas, expected_tas)

# --- height = 5000 m ---
h = 5000.0  # m
T, p, rho, a = atmosphere_isa(h)
# From TAS
tas = 200.0 * KT2MS  # m/s
expected_eas = 155.03684 * KT2MS  # m/s
expected_cas = 155.95551 * KT2MS  # m/s

eas = tas2eas(tas, rho)
cas = tas2cas(tas, rho, p)
@test isapprox(eas, expected_eas, atol=1e-4)
@test isapprox(cas, expected_cas, atol=1e-4)

# From CAS
cas = 200.0 * KT2MS  # m/s
expected_eas = 198.101308 * KT2MS  # m/s
expected_tas = 255.553851 * KT2MS  # m/s

eas = cas2eas(cas, rho, p)
tas = cas2tas(cas, rho, p)
@test isapprox(eas, expected_eas, atol=1e-4)
@test isapprox(tas, expected_tas, atol=1e-4)

# From EAS
eas = 200.0 * KT2MS  # m/s
expected_cas = 201.95290 * KT2MS  # m/s
expected_tas = 258.00319 * KT2MS  # m/s

cas = eas2cas(eas, rho, p)
tas = eas2tas(eas, rho)
@test isapprox(cas, expected_cas, atol=1e-4)
@test isapprox(tas, expected_tas, atol=1e-4)

# --- Velocities from qc ---
# Values from http://www.aerospaceweb.org/design/scripts/atmosphere/
# geometric altitude = 8010.0807  # m  --> geopotential = 8000.0  # m
# velocity = 100  # m/s
h = 8000.0  # m
T, p, rho, a = atmosphere_isa(h)
tas = 100  # m/s

qinf = 38295.5172  # Pa (Total Head)
qc = qinf - p
exp_cas = 66.0304  # m/s
exp_tas = 100  # m/s
exp_eas = 65.4758  # m/s

eas = qc2eas(qc, p)
tas_ = qc2tas(qc, rho, p)
cas = qc2cas(qc)

@test isapprox(cas, exp_cas, atol=1e-3)
@test isapprox(eas, exp_eas, atol=1e-3)
@test isapprox(tas_, exp_tas, atol=1e-3)

# --- Dynamic pressure ---
# Values from http://www.aerospaceweb.org/design/scripts/atmosphere/
# geometric altitude = 8010.0807  # m  --> geopotential = 8000.0  # m
# velocity = 100  # m/s
h = 8000.0  # m
T, p, rho, a = atmosphere_isa(h)
tas = 100  # m/s

exp_qinf_inc = 2625.8356  # Pa (True Dynamic Pressure)
qinf_inc = incompressible_qinf(tas, rho)
@test isapprox(qinf_inc, exp_qinf_inc, atol=0.0001)

M = tas / a
exp_qinf_comp = 38295.5172  # Pa (Total Head)
qinf_comp = compressible_qinf(M, p)
@test isapprox(qinf_comp, exp_qinf_comp, atol=0.01)

M = [0.6, 1.3, 3.0]         # Ref=> Aerodynamics, Anderson, page 550
p = 1
exp_qinf_comp = [1.276, 2.714, 12.06]
qinf_comp = []
for i = 1:length(exp_qinf_comp)
    push!(qinf_comp, compressible_qinf(M[i], p))
end
@test isapprox(qinf_comp, exp_qinf_comp, atol=0.01)

# --- tas, alpha, beta from body ---
u, v, w = 100, 0, 0  # m/s
rv = tas_alpha_beta_from_uvw(u, v, w)
@test isapprox(rv, [100, 0, 0])

u, v, w = 10, 0, 10  # m/s
rv = tas_alpha_beta_from_uvw(u, v, w)
@test isapprox(rv, [sqrt(200), atan(1), 0])

u, v, w = 10, 10, 0  # m/s
rv = tas_alpha_beta_from_uvw(u, v, w)
@test isapprox(rv, [sqrt(200), 0, asin(10/sqrt(200))])


text: my_f(x, y) = 2x+y

max_stars_repo_name: csquires/MyExample

max_issues_repo_name: csquires/MyExample

max_forks_repo_name: csquires/MyExample


text:
#Choose the action based on the policy and belief
#Performs search of Tree returning best action from top node
#Can return tree info as second paramter
function action_info(p::POMCPPlanner, b; tree_in_info=false)
    local a::action_type(p.problem)
    info = Dict{Symbol, Any}()
    try
        tree = POMCPTree(p.problem, p.solver.tree_queries)
        a = search(p, b, tree, info)
        p._tree = Nullable(tree)
        if p.solver.tree_in_info || tree_in_info
            info[:tree] = tree
        end
    catch ex
        # Note: this might not be type stable, but it shouldn't matter too much here
        a = convert(action_type(p.problem), default_action(p.solver.default_action, p.problem, b, ex))
        info[:exception] = ex
    end
    return a, info
end

#Generic action selection based on policy and current belief
#returns best action, after doing search of tree (first is because could return tree info as second parameter)
action(p::POMCPPlanner, b) = first(action_info(p, b))

#Takes a policy and searches the tree based on initial belief
#Returns best action
#This calls simulate() for each node open to search in the tree.
function search(p::POMCPPlanner, b, t::POMCPTree, info::Dict)
    all_terminal = true
    i = 0
    start_us = CPUtime_us()
    #Cycle for defined number of loops of sample+simulate to update policy
    for i in 1:p.solver.tree_queries
        #Check if search time has expired
        if CPUtime_us() - start_us >= 1e6*p.solver.max_time
            break
        end
        s = rand(p.rng, b) #SAMPLE State from belief <-----------!!!!!!!!!!!!!!!!!!!!!
        #Check for terminal state in the problem;
        if !POMDPs.isterminal(p.problem, s)
            #Not terminal, so simulate the random state selected, and add to tree;
            simulate(p, s, POMCPObsNode(t, 1), p.solver.max_depth)
            all_terminal = false
        end
    end
    #Record search information
    info[:search_time_us] = CPUtime_us() - start_us
    info[:tree_queries] = i

    #Account for nothing to search.
    if all_terminal
        throw(AllSamplesTerminal(b))
    end

    #Search the roots children nodes (actions) for the best value
    h = 1
    best_node = first(t.children[h])
    best_v = t.v[best_node]
    @assert !isnan(best_v)
    for node in t.children[h][2:end]
        if t.v[node] >= best_v
            best_v = t.v[node]
            best_node = node
        end
    end

    #Return best value's action
    return t.a_labels[best_node]
end

solve(solver::POMCPSolver, pomdp::POMDP) = POMCPPlanner(solver, pomdp)

#I ADDED
#struct Obs_Counts{S,A,O}
#    M::Dict{Tuple{S,A,O}, Int}
#end

#struct Trans_Counts{S,A}
#    M::Dict{Tuple{S,A,S}, Int}
#end

#Get the observation probability
#Base BAPOMDP state
mutable struct BAPOMDPState{S}
    s::S                        #State
    oc::Vector{Int}             #Multi dimensional array: Array{Int}(S,A,O)
    tc::Vector{Int}             #Multi dimensional array: Array{Int}(S,A,S)
end
#TODO: Improve to handle continuous, or infinite number of states
#    oc::Dict{Tuple{S,A,O},Int}  #Observation counts
#    tc::Dict{Tuple{S,A,S},Int}  #Transition counts

function obs_count_prob(p::POMDP,s::BAPOMDPState,a,o)
    s_idx = find(x -> x==s.s,states(p))
    a_idx = find(x -> x==a,actions(p))
    o_idx = find(x -> x==o,observations(p))
    return first(s.oc[s_idx,a_idx,o_idx]/sum(s.oc[s_idx,a_idx,:])) #convert to scalar
end

#Get the transition probability
function trans_count_prob(p::POMDP,s::BAPOMDPState,a,sp::BAPOMDPState)
    s_idx = find(x -> x==s.s,states(p))
    a_idx = find(x -> x==a,actions(p))
    sp_idx = find(x -> x==sp.s,states(p))
    return first(s.tc[s_idx,a_idx,sp_idx]/sum(s.tc[s_idx,a_idx,:])) #convert to scalar
end

#Get the distribution of transition probabilities
function trans_count_dist(p::POMDP,s::BAPOMDPState,a)
    s_idx = find(x -> x==s.s,states(p))
    a_idx = find(x -> x==a,actions(p))
    return reshape(s.tc[s_idx,a_idx,:]/sum(s.tc[s_idx,a_idx,:]),(n_states(p)))
end

#Initialize the BAPOMDP state
function initiate_state(p::POMDP,s)
    tc = ones(Int,POMDPs.n_states(p),POMDPs.n_actions(p),POMDPs.n_states(p))
    oc = ones(Int,POMDPs.n_states(p),POMDPs.n_actions(p),POMDPs.n_observations(p))
    return BAPOMDPState{state_type(p)}(s,oc,tc)
end

#Copy BAPOMDP state
function copy(s::BAPOMDPState)
    return BAPOMDPState(s.s,Base.copy(s.oc),Base.copy(s.tc))
end

#Increment the counts for the new sp
function increment_trans_obs_counts(p::POMDP,s::BAPOMDPState,a,o,sp::BAPOMDPState)
    s_idx = find(x -> x==s.s,states(p))
    a_idx = find(x -> x==a,actions(p))
    sp_idx = find(x -> x==sp.s,states(p))
    o_idx = find(x -> x==o,observations(p))
    sp.tc[s_idx,a_idx,sp_idx] += 1 #only can do transition at this time, observation will be updated later
    sp.oc[s_idx,a_idx,o_idx] += 1 #only can do transition at this time, observation will be updated later
end

#Rollout probably needs the current counts to do the updater()
#Actually do a generate sp from s and action for particle filter
# See https://github.com/JuliaPOMDP/ParticleFilters.jl
#See the SimpleParticleFilter update function calls these two functions
#https://github.com/JuliaPOMDP/ParticleFilters.jl/blob/master/src/ParticleFilters.jl
#Add all probability to the obs_weight, and the generate_s should just be a random state
function ParticleFilters.generate_s(model::POMDP,s::BAPOMDPState,a,rng::AbstractRNG)
    sp = copy(s)
    s_index = rand(Categorical(trans_count_dist(model,s,a))) #Handle states that are not integers
    sp.s = states(model)[s_index]
    return sp
end

#OUtput of this is put into a WeightedParticleBelief, paired with the output of the generate_s() above
#only just the probability of sp + o, given a + s; P(s',o | a,s)
#NOTE: Modifying the sp counts.
function ParticleFilters.obs_weight(model::POMDP,a,s::BAPOMDPState,sp::BAPOMDPState,o)
    increment_trans_obs_counts(model,s,a,o,sp)
    return obs_count_prob(model,s,a,o) * trans_count_prob(model,s,a,sp) #merge the observation and transition probabilities together
end

#Also do a generate observation probability distribution from new state sp and action taken.
#Don't need, just need obs_weight; unless it is unweighted
#function ParticleFilters.observation(model::,a,sp)
#    obsDist = Array{Float}(n_obs(O)) #?????
#    return obsDist
#end
#Where is the transition and observation probabilities used??

#struct SA_Model{S}
#    s::Vector{Tuple{S,Trans_Counts,Obs_Counts}}
#end

#function generate_sa_model(t::Trans_Counts,o::Obs_Counts)
#end
#END ADDED

#Called per node visited, and with a specific state
#Return the roll-up value of the an observation reward....., plus node at each
function simulate(p::POMCPPlanner, s, hnode::POMCPObsNode, steps::Int)
    #Checks for final step depth to expand, or if state is terminal for the problem
    #Then end value accumulation
    if steps == 0 || isterminal(p.problem, s)
        return 0.0
    end

    t = hnode.tree
    h = hnode.node

    #Get log of total number of times visited this observation node (h)
    ltn = log(t.total_n[h])

    #Keep track of the equal value best nodes.
    best_nodes = empty!(p._best_node_mem)
    best_criterion_val = -Inf
    #Search each action node of this observation node for best criterion value
    for node in t.children[h]
        #get total number of times visited this action node
        n = t.n[node]

        #Decide if should absolutely pick this node, or pick it based on its value, or pick it based on exploration
        if n == 0 && ltn <= 0.0
            criterion_value = t.v[node]
        elseif n == 0 && t.v[node] == -Inf
            criterion_value = Inf
        else
            criterion_value = t.v[node] + p.solver.c*sqrt(ltn/n)
        end

        #See if better node value, or same node value
        if criterion_value > best_criterion_val
            best_criterion_val = criterion_value
            empty!(best_nodes)
            push!(best_nodes, node)
        elseif criterion_value == best_criterion_val
            push!(best_nodes, node)
        end
    end
    #Sample the best node (if multiple of them with same value) randomly
    ha = rand(p.rng, best_nodes)
    a = t.a_labels[ha]

    #Given the best action a and the sampled state in the belief state, then get the s', o, r
    sp, o, r = generate_sor(p.problem, s, a, p.rng)

    #Look-up the observation node under the current action node, else return 0
    hao = get(t.o_lookup, (ha, o), 0)
    if hao == 0
        #No observation node under the action node, so add new observation node
        hao = insert_obs_node!(t, p.problem, ha, o)
        #Estimate the value based on a state and observation node.
        #This is where a random rollout is performed to determine value, with the current policy
        #Dont think I need to add BA information in the rollout estimation
        v = estimate_value(p.solved_estimator,
                           p.problem,
                           sp,
                           POMCPObsNode(t, hao),
                           steps-1)
        R = r + discount(p.problem)*v
    else
        #Has observation node under the action node, so
        R = r + discount(p.problem)*simulate(p, sp, POMCPObsNode(t, hao), steps-1)
    end

    #Increment number of times visited for observation and best action node picked.
    t.total_n[h] += 1
    t.n[ha] += 1

    #Add the current reward, plus rollout of a new observation node, or the value of the next best action after this observation node
    # minus the current value, divided by total times visited.
    t.v[ha] += (R-t.v[ha])/t.n[ha]

    #Return future reward + current reward
    return R
end


text:
using Stats
# using Logging

function estimateNGFEpsilon(referenceImage::regImage;cutoffPercent = 80)
  m = size(referenceImage.data)
  h = referenceImage.voxelsize
  Rc = referenceImage.data

  shortFiniteDiffX = spdiagm((-ones(m[1]-1,1),ones(m[1]-2,1)),[0, 1],m[1]-1,m[1])/(2*h[1])
  #shortFiniteDiffX[1] = shortFiniteDiffX[2]
  #shortFiniteDiffX[end] = shortFiniteDiffX[end-1]

  shortFiniteDiffY = spdiagm((-ones(m[2]-1,1),ones(m[2]-2,1)),[0, 1],m[2]-1,m[2])/(2*h[2])
  #shortFiniteDiffY[1] = shortFiniteDiffY[2]
  #shortFiniteDiffY[end] = shortFiniteDiffY[end-1]

  averageX = spdiagm((0.5*ones(m[1]-2,1),0.5*ones(m[1]-1,1)),[-1, 0],m[1],m[1]-1)
  averageY = spdiagm((ones(m[2]-2,1)*0.5,ones(m[2]-1,1)*0.5),[-1, 0],m[2],m[2]-1)

  G1 = sparse(kron(speye(m[2]),shortFiniteDiffX));
  G2 = sparse(kron(shortFiniteDiffY,speye(m[1])));

  AvgX = sparse(kron(speye(m[2]),averageX));
  AvgY = sparse(kron(averageY,speye(m[1])));

  gradTx  = G1*referenceImage.data[:]
  gradTy  = G2*referenceImage.data[:]

  quantileX = percentile(AvgX * (gradTx.^2),cutoffPercent)
  quantileY = percentile(AvgY * (gradTy.^2),cutoffPercent)

  estimatedEps = sqrt(0.5*(quantileX + quantileY))

  lengthGT = sqrt(sum(AvgX * (gradTx.^2) + AvgY * (gradTy.^2), 2) + estimatedEps^2) #epsilon-norm

  vectorlength = reshape(lengthGT,(size(referenceImage)))

  try
    figure(51)
    clf()
    showImage(vectorlength)
    s = @sprintf("estimated epsilon: %3.3e", estimatedEps)
    title(s)
    colorbar()
    figure(52)
    clf()
    plt.hist(lengthGT,bins=100)
  catch err
    warn("Error while plotting. This is normal if it happens in the automated test framework.")
    println(err)
  end

  return (estimatedEps,vectorlength)

end


text:
using PyCall

np = pyimport("numpy")

P = pyimport("numpy.polynomial")
@pydef mutable struct Doubler <: P.Polynomial
    function __init__(self, x=10)
        self.x = x
    end
    my_method(self, arg1::Number) = arg1 + 20
    x2.get(self) = self.x * 2
    function x2.set!(self, new_val)
        self.x = new_val / 2
    end
end

@show Doubler().x2

@pydef mutable struct MyClass
    function __init__(self, x)
        self.x = x
    end
end

#=
py"""
class MyClass():
    def __init__(self, x):
        self.x = x
    def edit(self, v):
        self.x[0] = v
        return self.x
"""
=#

jo = [1 2 3; 4 5 6]
mc = MyClass(jo)
x = PyArray(PyObject(mc.edit(-10)))
@show x
x[2]=-200
@show x
@show mc.x

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_nq_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_nq")
JLLWrappers.@generate_main_file("GAP_pkg_nq", UUID("00bf61e5-c0f9-5a3f-ba8b-68a621695b53"))
end  # module GAP_pkg_nq_jll


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

@testset "init" begin
    prob = DDEProblemLibrary.prob_dde_constant_1delay_ip
    prob_scalar = DDEProblemLibrary.prob_dde_constant_1delay_scalar

    inferred = [BS3(), Tsit5(), RK4(), Vern6()]
    for alg in inferred
        ddealg = MethodOfSteps(alg)

        @test_broken @inferred init(prob, ddealg)
        @test_broken @inferred init(prob_scalar, ddealg)
    end

    notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()]
    for alg in notinferred
        ddealg = MethodOfSteps(alg)

        @test_broken @inferred init(prob, ddealg)
        @test_broken @inferred init(prob_scalar, ddealg)
    end
end

@testset "discontinuity_time" begin
    prob_inplace = DDEProblemLibrary.prob_dde_constant_1delay_ip
    prob_scalar = DDEProblemLibrary.prob_dde_constant_1delay_scalar

    for prob in (prob_inplace, prob_scalar)
        int = init(prob, MethodOfSteps(Tsit5()))
        @inferred DelayDiffEq.discontinuity_time(int, (u, p, t) -> 1.0, 0.0, (0.5, 1.5))
    end
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
export hard, soft, soft_intern, soft2, soft2_intern, soft_unrolled_intern, soft_unrolled
export soft_fcm, sof_fcm_intern
export exchange!, substitute!, substitute

# hard coded function calls in inner loop
function hard(s::String)
    io = IOBuffer()
    for c in s
        c = lowercase(c)
        c = c + 20
        c = uppercase(c)
        c = c - 20
        write(io, c)
    end
    String(take!(io))
end

# the same functionality, more flexible with a vararg list of
# functions as arguments
# performance is 4 times worse - many allocations
function soft_intern(s::String, funlist::Function...)
    io = IOBuffer()
    for c in s
        for fun in funlist
            c = fun(c)::Char
        end
        write(io, c)
    end
    String(take!(io))
end

# make use of the flexible implementation
function soft(s::String)
    soft_intern(s, lowercase, c->c+20, uppercase, c->c-20)
end

# the same functionality, more flexible with a vector of
# functions as arguments
# performance is 4 times worse - many allocations
function soft2_intern(s::String, funlist::Vector{Function})
    io = IOBuffer()
    for c in s
        for fun in funlist
            c = fun(c)::Char
        end
        write(io, c)
    end
    String(take!(io))
end

# make use of the flexible implementation
function soft2(s::String)
    funlist = [lowercase, c->c+20, uppercase, c->c-20]
    soft2_intern(s, funlist)
end

# an attempt to improve perfomance using `eval`
# works!
function soft_unrolled_intern(funlist::Function...)

    function code_for_funcall(fun::Function)
        :( c = $(fun)(c) )
    end
    function unrolled()
        Expr(:block, code_for_funcall.(funlist)...)
    end

    ex = :(
       (s::String) ->
       let c::Char, io = IOBuffer()
           for c in s
               $(unrolled())
               write(io, c)
           end
           String(take!(io))
       end
    )
    f = eval(ex)
    (x) -> Base.invokelatest(f, x)
end

"""

    soft_unrolled(s::String)

The following data show the timing, which is achieved by the different implementations.
When the function is generated, and executed separately, perfomance is
as good as the hard-coded version!

```jldoctest
julia> srand(0); s = randstring(10^7);

julia> @time hard(s);
  0.497349 seconds (13.58 k allocations: 16.785 MiB)

julia> @time hard(s);
  0.469577 seconds (28 allocations: 16.003 MiB)

julia> @time soft(s);
  2.087998 seconds (9.70 M allocations: 165.476 MiB, 2.98% gc time)

julia> @time soft(s);
  2.015281 seconds (9.67 M allocations: 171.554 MiB, 2.53% gc time)

julia> @time soft_unrolled(s);
  0.568373 seconds (38.56 k allocations: 18.327 MiB)

julia> @time soft_unrolled(s);
  0.506165 seconds (19.58 k allocations: 17.188 MiB, 4.41% gc time)

julia> @time ex = soft_unrolled_intern(lowercase, c->c+20, uppercase, c->c-20)
  0.020691 seconds (6.41 k allocations: 413.712 KiB)
#5 (generic function with 1 method)

julia> @time ex(s);
  0.565612 seconds (19.58 k allocations: 17.129 MiB, 9.19% gc time)

julia> @time ex(s);
  0.477627 seconds (37 allocations: 16.003 MiB)

julia> @time ex(s);
  0.479458 seconds (37 allocations: 16.003 MiB)
"""
function soft_unrolled(s::String)
    soft_unrolled_intern(lowercase, c->c+20, uppercase, c->c-20)(s)
end

# make that generic
function substitute!(outer::Expr, critical, list)
    function code_for_funcall(element)
        exchange!(deepcopy(critical), :( :REPLACE ), :( $element ))
    end

    function unrolled()
        Expr(:block, code_for_funcall.(list)...)
    end

    exchange(outer, :( :BLOCK ), unrolled())
end

substitute(outer::Expr, critical, list) = substitute!(deepcopy(outer), critical, list)

function exchange!(code::Expr, dummy, replacement)
    if dummy == code
        replacement
    else
        for i = 1:length(code.args)
            code.args[i] = exchange(code.args[i], dummy, replacement)
        end
    end
    code
end

function exchange!(code::Any, dummy::Any, replacement::Any)
    code == dummy ? replacement : code
end

"""
    An attempt to use 'first-class-methods'
"""

struct FCM{T<:Function} <: Function
    f::T
end

function (f::FCM{T})(a::Char)::Char where T<:Function
    ff = f.f
    ff(a)::Char
end

@inline function soft_fcm_intern(s::String, funlist::FCM...)
    io = IOBuffer()
    for c in s
        for fun in funlist
            c = fun(c)::Char
        end
        write(io, c)
    end
    String(take!(io))
end

# make use of the flexible implementation
function soft_fcm(s::String)
    soft_fcm_intern(s, FCM.((lowercase, c->c+20, uppercase, c->c-20))...)
end


text:
using DifferentialEquations

function main()
    t0, tf = 0, 10
    x0 = [1.0, 2, 3]
    """
    dx: next x
    """
    @Loggable function dynamics!(dx, x, p, t; u)
        @log x
        dx .= 0.99*x + u
        @onlylog u_next = dx
    end
    simulator = Simulator(x0, apply_inputs(dynamics!; u=zeros(3));
                          Problem=:Discrete,
                          tf=tf,  # default step length = 1 for Discrete problem
                         )
    df = solve(simulator)
end


text:
# important: set d=1 in BridgeLandmarks.jl
using BridgeLandmarks
using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD2

workdir = @__DIR__
println(workdir)
cd(workdir)

n = 3
nshapes = 1
T = 1.0

xobs0 = [PointF(-0.5), PointF(0.0), PointF(0.1)]
#xobsT = [exp.(q0[i])  for i in 1:n ] .- PointF(3.0)
xobsT = [PointF(-0.5), PointF(0.2), PointF(1.0)]

JLD2.@save "data_exp1-1D.jld2" xobs0 xobsT n nshapes


text:
using Seaborn
using Test

# Let's just make sure this doesn't error
rugplot([1,3,5])

# Issue #11
Seaborn.set(style="ticks")