text:
function optimize(tuning_run::Run)
    stop       = RemoteChannel(()->Channel{Bool}(1))
    report     = RemoteChannel(()->Channel{Bool}(1))

    @spawn tuning_run.stopping_criterion(tuning_run.duration, stop)
    @spawn tuning_run.reporting_criterion(tuning_run.report_after, report)

    results    = initialize_search_tasks!(tuning_run)
    best       = get_new_best(results)
    iteration  = 1
    start_time = time()
    put!(tuning_run.channel, best)

    while !take!(stop)
        best                   = get_new_best(results, best)
        iteration             += 1
        best.current_iteration = iteration

        if take!(report)
            best.current_time = time() - start_time
            put!(tuning_run.channel, best)
        end
    end

    best.is_final     = true
    best.current_time = time() - start_time
    put!(tuning_run.channel, best)

    end_search_tasks!(results)
    close(tuning_run.channel)
    close(stop)
    close(report)
end


text:
const  √ = sqrt
const  ∑ = sum
const  ⊗ = kron

ext: jl


text:
using Base.Broadcast: broadcastable
using ChainRules
using ChainRulesCore
using ChainRulesTestUtils
using ChainRulesTestUtils: rand_tangent, _fdm
using Compat: hasproperty, only, cispi, eachcol
using FiniteDifferences
using LinearAlgebra
using LinearAlgebra.BLAS
using LinearAlgebra: dot
using Random
using StaticArrays
using Statistics
using Test

Random.seed!(1) # Set seed that all testsets should reset to.

function include_test(path)
    println("Testing $path:")  # print so TravisCI doesn't timeout due to no output
    @time include(path)  # show basic timing, (this will print a newline at end)
end

println("Testing ChainRules.jl")
@testset "ChainRules" begin
    include_test("test_helpers.jl")
    println()
    @testset "rulesets" begin
        @testset "Core" begin
            include_test("rulesets/Core/core.jl")
        end

        @testset "Base" begin
            include_test("rulesets/Base/base.jl")
            include_test("rulesets/Base/fastmath_able.jl")
            include_test("rulesets/Base/evalpoly.jl")
            include_test("rulesets/Base/array.jl")
            include_test("rulesets/Base/arraymath.jl")
            include_test("rulesets/Base/indexing.jl")
            include_test("rulesets/Base/mapreduce.jl")
            include_test("rulesets/Base/sort.jl")
        end
        println()

        @testset "Statistics" begin
            include_test("rulesets/Statistics/statistics.jl")
        end
        println()

        @testset "LinearAlgebra" begin
            include_test("rulesets/LinearAlgebra/dense.jl")
            include_test("rulesets/LinearAlgebra/norm.jl")
            include_test("rulesets/LinearAlgebra/matfun.jl")
            include_test("rulesets/LinearAlgebra/structured.jl")
            include_test("rulesets/LinearAlgebra/symmetric.jl")
            include_test("rulesets/LinearAlgebra/factorization.jl")
            include_test("rulesets/LinearAlgebra/blas.jl")
            include_test("rulesets/LinearAlgebra/lapack.jl")
        end
        println()

        @testset "Random" begin
            include_test("rulesets/Random/random.jl")
        end
        println()
    end
end


text:
# Extract the bow shock location from 2D equatorial run.
# Courtesy of Vertti Tarvus, reimplemented from Python.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot

# Upstream solar wind temperature
Tsw = 0.5e6 #[K]

file = "bulk.0001347.vlsv"
meta = load(file)

x = LinRange{Float32}(meta.coordmin[1], meta.coordmax[1], meta.ncells[1])
y = LinRange{Float32}(meta.coordmin[2], meta.coordmax[2], meta.ncells[2])

# Obtain thermal temperature
T = meta["T"]
T = reshape(T, meta.ncells[1], meta.ncells[2])

x_crossing = zeros(Float32, meta.ncells[2])
y_crossing = y

# Extract bow shock location from the 1st point which fulfills the threshold: T > 4 * Tsw
# Note: this won't work if the bow shock does not cover the whole domain y extension!
for j in 1:meta.ncells[2] # scan in y direction
   ind_ = findlast(>(4*Tsw), @view T[:,j]) # count from upstream
   x_crossing[j] = x[ind_]
end

# Julia is column-major, while Python is row-major
imshow(T', extent=(x[1], x[end], y[1], y[end]), origin="lower", cmap=plt.get_cmap("ocean"))
plot(x_crossing, y_crossing, "r")
axis("scaled")
savefig("bs_temp_test.png")


text:
using FymEnvs
using LinearAlgebra

if !isdefined(Main, :TestEnvs)
    include("custom_env.jl")  # to avoid conflict
    using .TestEnvs
end

function test_custom_fym()
    fym = TestEnv(; max_t=0.1)

    env = fym.env
    log_dir = "data/test"
    file_name = "custom.jld2"
    reset!(env)
    _sample(fym, nothing, log_dir, file_name)
    data, info = load(joinpath(log_dir, file_name),
                      with_info=true)
    print(info)

    reset!(env)
    _sample2(fym, nothing, log_dir, file_name)
    data, info = load(joinpath(log_dir, file_name),
                      with_info=true)
    print(info)
end

function _sample(fym::FymEnv, agent, log_dir, file_name)
    env = fym.env
    logger = Logger(log_dir=log_dir, file_name=file_name)
    obs = observe_flat(env)
    i = 0
    @time while true
        if agent == nothing
            action = nothing
            if time(env.clock) > 0.5
                fym.controller.K = - 3*Matrix(I, 3, 3)
            end
        end
        next_obs, reward, done, info = env.step(action=action)
        record(logger, info)
        obs = next_obs
        i += 1
        if done
            config = Dict(
                          "name" => env.name,
                       )
            set_info!(logger, config)
            break
        end
    end
    close!(env)
    close!(logger)
end

function _sample2(fym::FymEnv, agent, log_dir, file_name)
    env = fym.env
    logger = Logger(log_dir=log_dir, file_name=file_name)
    obs = observe_flat(env)
    i = 0
    @time while true
        if agent == nothing
            action = nothing
            if time(env.clock) > 0.5
                fym.controller.K = - 3*Matrix(I, 3, 3)
            end
        end
        next_obs, reward, done, info = env.step(action=action)
        record(logger, info)
        obs = next_obs
        i += 1
        if done
            config = Dict(
                          "name" => Dict(),
                       )
            set_info!(logger, config)
            break
        end
    end
    close!(env)
    close!(logger)
end

test_custom_fym()


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")
include("test_problem_and_solver.jl")

function error_for_plane_interface(
    x0,
    normal,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
    eta,
)
    L = 1.0
    W = 1.0
    lambda, mu = 1.0, 2.0
    dx = 1.0 / nelmts
    penalty = penaltyfactor / dx * (lambda + mu)
    alpha = 0.1
    stiffness = CutCellDG.HookeStiffness(lambda, mu, lambda, mu)

    elasticitybasis = LagrangeTensorProductBasis(2, polyorder)

    mergedmesh, cellquads, facequads, interfacequads, levelset =
        construct_mesh_and_quadratures(
            [L, W],
            nelmts,
            elasticitybasis,
            x -> plane_distance_function(x, normal, x0),
            numqp,
        )

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        elasticitybasis,
        cellquads,
        stiffness,
        mergedmesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        elasticitybasis,
        interfacequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        x -> displacement(alpha, x),
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        x -> onboundary(x, L, W),
        penalty,
    )
    CutCellDG.assemble_body_force!(
        sysrhs,
        x -> body_force(lambda, mu, alpha, x),
        elasticitybasis,
        cellquads,
        mergedmesh,
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mergedmesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mergedmesh)

    solution = matrix \ rhs
    nodaldisplacement = reshape(solution, 2, :)

    err = mesh_L2_error(
        nodaldisplacement,
        x -> displacement(alpha, x),
        elasticitybasis,
        cellquads,
        mergedmesh,
    )
end

function test_plane_interface_linear_elements()

    x0 = [0.8, 0.0]
    interfaceangle = 30.0
    normal = [cosd(interfaceangle), sind(interfaceangle)]
    powers = [3, 4, 5]
    nelmts = [2^p + 1 for p in powers]
    polyorder = 1
    numqp = required_quadrature_order(polyorder)
    penaltyfactor = 1e2
    eta = 1

    err = [
        error_for_plane_interface(
            x0,
            normal,
            ne,
            polyorder,
            numqp,
            penaltyfactor,
            eta,
        ) for ne in nelmts
    ]

    u1err = [er[1] for er in err]
    u2err = [er[2] for er in err]

    dx = 1.0 ./ nelmts
    u1rate = convergence_rate(dx, u1err)
    u2rate = convergence_rate(dx, u2err)

    @test all(u1rate .> 1.95)
    @test all(u2rate .> 1.95)

end

function test_plane_interface_quadratic_elements()
    x0 = [0.1, 0.0]
    interfaceangle = -28.0
    normal = [cosd(interfaceangle), sind(interfaceangle)]
    powers = [3, 4, 5]
    nelmts = [2^p + 1 for p in powers]
    polyorder = 2
    numqp = required_quadrature_order(polyorder)
    penaltyfactor = 1e2
    eta = 1

    err = [
        error_for_plane_interface(
            x0,
            normal,
            ne,
            polyorder,
            numqp,
            penaltyfactor,
            eta,
        ) for ne in nelmts
    ]

    u1err = [er[1] for er in err]
    u2err = [er[2] for er in err]

    dx = 1.0 ./ nelmts
    u1rate = convergence_rate(dx, u1err)
    u2rate = convergence_rate(dx, u2err)

    @test all(u1rate .> 2.95)
    @test all(u2rate .> 2.95)
end

test_plane_interface_linear_elements()
test_plane_interface_quadratic_elements()


text:
using Revise
using PyPlot
using PyCall
using Statistics
using LinearAlgebra
using Infiltrator
using IterativeSolvers
using Bem2d

"""
    plotbcsUeff(els, titlestring)

Plot boundary conditions and Ueff
"""
function plotbcsUeff(els, bcs, Ueff, suptitlestring)
    figure(figsize=(20, 20))
    subplot(2, 2, 1)
    quiver(els.xcenter[1:els.endidx], els.ycenter[1:els.endidx],
           bcs[1:2:end], bcs[2:2:end])
    gca().set_aspect("equal")
    title("bcs")

    subplot(2, 2, 2)
    plot(bcs[1:2:end], "rx", label="x")
    plot(bcs[2:2:end], "b+", label="y")
    legend()
    title("bcs")

    subplot(2, 2, 3)
    quiver(els.xcenter[1:els.endidx], els.ycenter[1:els.endidx],
           Ueff[1:2:end], Ueff[2:2:end])
    gca().set_aspect("equal")
    title("Ueff")

    subplot(2, 2, 4)
    plot(Ueff[1:2:end], "rx", label="x")
    plot(Ueff[2:2:end], "b+", label="y")
    legend()
    title("Ueff")

    suptitle(suptitlestring)
end

"""
    plotgeometry(els, titlestring)

els structure geometry plotting for diagnostics
"""
function plotgeometry(els, titlestring)
    scale = 500
    figure(figsize=(15,10))
    for i in 1:els.endidx
        plot([els.x1[i], els.x2[i]], [els.y1[i], els.y2[i]], "-g")
        quiver(els.xcenter[i], els.ycenter[i],
               els.xnormal[i], els.ynormal[i],
               width=1e-3, scale=1e2)
        text(els.xcenter[i] + scale * els.xnormal[i],
             els.ycenter[i] + scale * els.ynormal[i],
             string(i), fontsize=10)
    end
    title(titlestring)
    gca().set_aspect("equal")
end

"""
    weldedcircle()

Comparing homogeneous and welded circle BEM solutions
"""
function weldedcircle()
    close("all")
    PLOTGEOMETRY = false
    mu1 = 3e10
    nu1 = 0.25
    mu2 = .5 * mu1
    nu2 = 0.25

    npts = 100
    offset = 100 # meters

    # Element geometries and data structures for the homogeneous circle case
    els1 = Elements(Int(1e5))
    nels = 20
    r = 10e3
    x1, y1, x2, y2 = discretizedarc(deg2rad(0), deg2rad(180), r, nels)
    addelsez!(els1, x1, y1, x2, y2, "T")
    x1, y1, x2, y2 = discretizedarc(deg2rad(180), deg2rad(360), r, nels)
    addelsez!(els1, x1, y1, x2, y2, "B")
    idx1 = getidxdict(els1)
    # PLOTGEOMETRY && plotgeometry(els1, "Circle boundaries and normals")

    T_B_BT, _ = PUTC(slip2dispstress, els1, idx1["B"], 1:els1.endidx, mu1, nu1)
    _, H_T_BT = PUTC(slip2dispstress, els1, idx1["T"], 1:els1.endidx, mu1, nu1)
    bcs1 = zeros(2*els1.endidx)
    bcs1[2*10] = 1.0
    bcs1[2*11] = 1.0
    Ueff1 = [H_T_BT ; T_B_BT] \ bcs1
    # plotbcsUeff(els1, bcs1, Ueff1, "homogeneous circle")
    xgrid1, ygrid1 = obsgrid(-r, -r, r, r, npts)
    U1, S1 = constdispstress(slip2dispstress, xgrid1, ygrid1, els1, 1:els1.endidx,
                             Ueff1[1:2:end], Ueff1[2:2:end], mu1, nu1)
    # plotfields(els1, reshape(xgrid1, npts, npts), reshape(ygrid1, npts, npts),
    #            U1, S1, "homogeneous circle")

    # Element geometries and data structures for the welded circle case
    els2 = Elements(Int(1e5))
    nels = 20
    r = 10e3
    x1, y1, x2, y2 = discretizedarc(deg2rad(0), deg2rad(180), r, nels)
    addelsez!(els2, x1, y1, x2, y2, "T")
    x1, y1, x2, y2 = discretizedline(-r, 0, r, 0, nels)
    addelsez!(els2, x1, y1, x2, y2, "midT")
    x1, y1, x2, y2 = discretizedarc(deg2rad(180), deg2rad(360), r, nels)
    addelsez!(els2, x1, y1, x2, y2, "B")
    x1, y1, x2, y2 = discretizedline(-r, 0, r, 0, nels) # Peculiar
    addelsez!(els2, x2, y2, x1, y1, "midB")
    idx2 = getidxdict(els2)
    PLOTGEOMETRY && plotgeometry(els2, "Welded circle boundaries and normals")

    TH = zeros(2 * els2.endidx, 2 * els2.endidx)
    T_midT_TmidT, H_midT_TmidT = PUTC(slip2dispstress, els2, idx2["midT"],
                                      [idx2["T"]; idx2["midT"]],
                                      mu1, nu1)
    T_midB_BmidB, H_midB_BmidB = PUTC(slip2dispstress, els2,
                                      idx2["midB"], [idx2["B"]; idx2["midB"]],
                                      mu2, nu2)
    _, H_T_TmidT = PUTC(slip2dispstress, els2, idx2["T"],
                        [idx2["T"]; idx2["midT"]], mu1, nu1)
    T_B_BmidB, _ = PUTC(slip2dispstress, els2, idx2["B"],
                        [idx2["B"]; idx2["midB"]], mu2, nu2)

    TH[1:40, 1:80] = T_midT_TmidT
    TH[1:40, 81:160] = -T_midB_BmidB
    TH[41:80, 1:80] = H_midT_TmidT
    TH[41:80, 81:160] = H_midB_BmidB
    TH[81:120, 1:80] = H_T_TmidT
    TH[121:160, 81:160] = T_B_BmidB
    bcs2 = zeros(2*els2.endidx)
    bcs2[2*50] = 1.0 # These are the indices after welding
    bcs2[2*51] = 1.0

    # Direct solve
    Ueff2 = TH \ bcs2
    UeffT2 = Ueff2[1:80]
    UeffB2 = Ueff2[81:160]

    # Volume visualization for two domain case
    Tidx = findall(x -> x > 0, ygrid1)
    Bidx = findall(x -> x < 0, ygrid1)
    Tx = xgrid1[Tidx]
    Ty = ygrid1[Tidx]
    Bx = xgrid1[Bidx]
    By = ygrid1[Bidx]
    UT2, ST2 = constdispstress(slip2dispstress, Tx, Ty, els2,
                               [idx2["T"]; idx2["midT"]],
                               UeffT2[1:2:end], UeffT2[2:2:end], mu1, nu1)
    UB2, SB2 = constdispstress(slip2dispstress, Bx, By, els2,
                               [idx2["B"]; idx2["midB"]],
                               UeffB2[1:2:end], UeffB2[2:2:end], mu2, nu2)
    U2 = [UB2 ; UT2] # Note B, T ordering
    S2 = [SB2 ; ST2] # Note B, T ordering
    plotfields(els2, reshape(xgrid1, npts, npts), reshape(ygrid1, npts, npts),
               U2, S2, "welded circle")
    # plotfields(els2, reshape(xgrid1, npts, npts), reshape(ygrid1, npts, npts),
    #            U2 .- U1, S2 .- S1, "residuals")
end
weldedcircle()


text:
using Flux, XConv, BSON, IterTools
using Flux: onecold, logitcrossentropy, @epochs
using CUDA
using Parameters: @with_kw
using ProgressMeter

CUDA.unsafe_free!(h::Flux.OneHotMatrix) = CUDA.unsafe_free!(h.data)

cd(dirname(@__FILE__))

if CUDA.has_cuda()
    device = gpu
    CUDA.allowscalar(false)
    iterator = CuIterator
    @info "Training on GPU"
else
    iterator = ()
    device = cpu
    @info "Training on CPU"
end

include("./datautils.jl")
include("./modelutils.jl")

@with_kw mutable struct Args
    name::String = "MNIST"
    batchsize::Int = 128
    η::Float64 = 3e-3
    η_fact::Float64 = .9
    progress_fact::Float64 = .001
    epochs::Int = 20
    splitr_::Float64 = 0.1
    probe_size::Int = 8
    mode::String = "TrueGrad"
    savepath::String = "./$(name)_bench"
end

function accuracy(test_data, m)
    correct, total = 0, 0
    for (x, y) in iterator(test_data)
        correct += sum(onecold(cpu(m(x)), 0:9) .== onecold(cpu(y), 0:9))
        total += size(y, 2)
    end
    test_accuracy = correct / total
    test_accuracy
end

augment(x) = x .+ device(0.1f0*randn(eltype(x), size(x)))
filename(args) = joinpath(args.savepath, "$(args.name)_conv_$(args.mode)_$(args.batchsize)_$(args.probe_size).bson")

function check(fname)
    isfile(fname) ? (BSON.@load fname acc) : (return false)
    acc < .1 ? rm(fname) : (return true)
    false
end

function train(; kws...)
    # Initialize the hyperparameters
    args = Args(; kws...)
    check(filename(args)) && return

    # Load the train, validation data
    train_data, val_data = get_train_data(args)
    test_data = get_test_data(args)

    @info("Constructing Model")
    # Defining the loss and accuracy functions
    m = build_model(args.name) |> device
    loss(x, y) = logitcrossentropy(m(x), y)

    @info("Training $(args.name) with batch size $(args.batchsize)")
    lhist = []

    # Defining the optimizer
    opt = ADAM(args.η)
    ps = Flux.params(m)

    # Initialize gradient mode
    XConv.initXConv(args.probe_size, args.mode)
    # Starting to train models
    p = Progress(length(train_data) * args.epochs)
    best = 0
    for epoch in 1:args.epochs
        Base.flush(Base.stdout)
        local l, acc
	acc = accuracy(val_data, m)
	n = length(train_data)
        for (x, y) in iterator(train_data)
            gs = Flux.gradient(ps) do
                l = loss(x, y)
                l
            end
            Flux.update!(opt, ps, gs)
            push!(lhist, l)
	    ProgressMeter.next!(p; showvalues = [(:loss, l), (:epoch, epoch), (:Mode, args.mode), (:ps, args.probe_size), (:η, args.η), (:accuracy, acc)])
	end

        validation_loss = 0f0
        for (x, y) in iterator(val_data)
            validation_loss += loss(x, y)
        end
        validation_loss /= length(val_data)
        acc_loc = accuracy(val_data, m)
	(acc_loc/acc - 1) < args.progress_fact && (args.η *= args.η_fact)
        @info "Epoch $epoch validation with $(XConv._params) loss = $(validation_loss), accuracy = $(acc_loc)"
    end
    acc = accuracy(test_data, m)
    BSON.@save filename(args) params=cpu.(params(m)) acc lhist
    GC.gc()
end

b_sizes = [2^i for i=5:11]
ps_sizes = [0..., [2^i for i=1:7]...]

for b in b_sizes
    for ps in ps_sizes
        η = 3e-3
        train(;η=η, epochs=10+2*log2(b),batchsize=b, probe_size=ps, name=d, mode= ps>0 ? "EVGrad" : "TrueGrad")
    end
end

nps, nb = length(ps_sizes)+1, length(b_sizes)
accs = Dict("MNIST" => zeros(nps, nb))
d = "MNIST"

for (i,b) in enumerate(b_sizes)
    @show b, d
    BSON.@load joinpath("./$(d)_bench/", "$(d)_conv_TrueGrad_$(b)_0.bson") acc
    accs[d][1, i] = acc
    for (j,ps) in enumerate(ps_sizes)
        BSON.@load joinpath("./$(d)_bench/", "$(d)_conv_EVGrad_$(b)_$(ps).bson") acc
        accs[d][j+1, i] = acc
    end
end

# ppri
header(i) = i==1 ? "| default   |" : "| ``ps=$(ps_sizes[i-1])`` |"

str_mnist = prod([header(j)*prod([" ``$(round(aa, digits=4))`` |" for aa in accs["MNIST"][j,:]])*"\n" for j=1:nps])

top = "|           |"*prod(["``B=$b``   |" for b in b_sizes])*"\n|"*"|:---------:"^(nb+1)*"|\n"

println("MNIST")
println(top*str_mnist)

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    RegularGridData(data, origin, spacing)

Regularly spaced `data` georeferenced with `origin` and `spacing`.
The `data` argument is a dictionary mapping variable names to Julia
arrays with the actual data.

`NaN` or `missing` values in the Julia arrays are interpreted as
non-valid. They can be used to mask the variables on the grid.

## Examples

Given `poro` and `perm` two 2-dimensional Julia arrays containing
values of porosity and permeability, the following code can be used
to georeference the data:

```julia
julia> data = Dict(:porosity => poro, :permeability => perm)
julia> RegularGridData(data, [0.,0.,0.], [1.,1.,1.])
```

Alternatively, one can omit `origin` and `spacing` for default
values of zeros and ones:

```julia
julia> RegularGridData{Float64}(data)
```

See also: [`RegularGrid`](@ref)
"""
struct RegularGridData{T<:Real,N} <: AbstractSpatialData{T,N}
  data::Dict{Symbol,<:AbstractArray}
  origin::NTuple{N,T}
  spacing::NTuple{N,T}

  # state fields
  dims::Dims{N}

  function RegularGridData{T,N}(data, origin, spacing) where {N,T<:Real}
    sizes = [size(array) for array in values(data)]
    @assert length(unique(sizes)) == 1 "data dimensions must be the same for all variables"
    @assert length(sizes[1]) == N "inconsistent number of dimensions for given origin/spacing"
    @assert all(spacing .> 0) "spacing must be positive"
    new(data, origin, spacing, sizes[1])
  end
end

RegularGridData(data::Dict{Symbol,<:AbstractArray}, origin::Vector{T}, spacing::Vector{T}) where {T<:Real} =
  RegularGridData{T,length(origin)}(data, (origin...,), (spacing...,))

RegularGridData{T}(data::Dict{Symbol,<:AbstractArray{<:Any,N}}) where {N,T<:Real} =
  RegularGridData{T,N}(data, (zeros(T,N)...,), (ones(T,N)...,))

coordinates(geodata::RegularGridData{T,N}) where {N,T<:Real} = Dict(Symbol("x$i") => T for i=1:N)

variables(geodata::RegularGridData) = Dict(var => eltype(array) for (var,array) in geodata.data)

npoints(geodata::RegularGridData) = prod(geodata.dims)

function coordinates!(buff::AbstractVector{T}, geodata::RegularGridData{T,N}, ind::Int) where {N,T<:Real}
  intcoords = CartesianIndices(geodata.dims)[ind]
  for i in 1:N
    @inbounds buff[i] = geodata.origin[i] + (intcoords[i] - 1)*geodata.spacing[i]
  end
end

value(geodata::RegularGridData, ind::Int, var::Symbol) = geodata.data[var][ind]

function Base.view(geodata::RegularGridData{T,N}, inds::AbstractVector{Int}) where {N,T<:Real}
  data = Dict(var => view(array, inds) for (var,array) in geodata.data)

  coords = Matrix{T}(N, length(inds))
  for ind in inds
    coords[:,ind] .= coordinates(geodata, ind)
  end

  PointSetData(data, coords)
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, geodata::RegularGridData{T,N}) where {N,T<:Real}
  dims = join(geodata.dims, "×")
  print(io, "$dims RegularGridData{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", geodata::RegularGridData{T,N}) where {N,T<:Real}
  println(io, geodata)
  println(io, "  origin:  ", geodata.origin)
  println(io, "  spacing: ", geodata.spacing)
  println(io, "  variables")
  varlines = ["    └─$var ($(eltype(array)))" for (var, array) in geodata.data]
  print(io, join(varlines, "\n"))
end


text:
"""
	subplot(fim=nothing; kwargs...)

Manage figure subplot configuration and selection.

Full option list at [`subplot`](http://docs.generic-mapping-tools.org/latest/subplot.html)

Parameters
----------

- **grid** :: [Type => Str | Tuple | Array]

    Specifies the number of rows and columns of subplots. Ex grid=(2,3)
- **F** | **dims** | **dimensions** | **size** | **sizes** :: [Type => Str | Tuple, NamedTuple]

    Specify the dimensions of the figure.
    [`-F`](http://docs.generic-mapping-tools.org/latest/subplot.html#f)
- **A** | **autolabel** | **fixedlabel** :: [Type => Str | number]

    Specify automatic tagging of each subplot. This sets the tag of the first, top-left subplot and others follow sequentially.
    [`-A`](http://docs.generic-mapping-tools.org/latest/subplot.html#a)
- $(GMT.opt_B)
- **C** | **clearance** :: [Type => Str | number]

    Reserve a space of dimension clearance between the margin and the subplot on the specified side. Settings specified under **begin** directive apply to all panels.
    [`-C`](http://docs.generic-mapping-tools.org/latest/subplot.html#c)
- $(GMT.opt_J)
- **M** | **margin** | **margins** :: [Type => Str]

    The margin space that is added around each subplot beyond the automatic space allocated for tick marks, annotations, and labels.
    [`-M`](http://docs.generic-mapping-tools.org/latest/subplot.html#m)
- $(GMT.opt_R)
- **SC** | **SR** | **col_axes** | **row_axes** :: [Type => Str | NamedTuple]

    Set subplot layout for shared axes. Set separately for rows (SR) and columns (SC).
    [`-S`](http://docs.generic-mapping-tools.org/latest/subplot.html#s)
- **T** | **title** :: [Type => Str]

    While individual subplots can have titles, the entire figure may also have a overarching title.
    [`-T`](http://docs.generic-mapping-tools.org/latest/subplot.html#t)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
"""
function subplot(fim=nothing; stop=false, kwargs...)

	FirstModern[1] = true			# To know if we need to compute -R in plot. Due to a GMT6.0 BUG

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode

	# In case :title exists we must use and delete it to avoid double parsing
	cmd = ((val = find_in_dict(d, [:T :title])[1]) !== nothing) ? " -T\"" * val * "\"" : ""
	cmd, opt_B, opt_J, opt_R = parse_BJR(d, cmd, "", false, " ")
	cmd, = parse_common_opts(d, cmd, [:params], true)
	cmd  = parse_these_opts(cmd, d, [[:M :margin :margins]])
	cmd  = add_opt(d, cmd, "A", [:A :autolabel],
                  (Anchor=("+J", arg2str), anchor=("+j", arg2str), label="", clearance=("+c", arg2str), fill=("+g", add_opt_fill), pen=("+p", add_opt_pen), offset=("+o", arg2str), roman="_+r", Roman="_+R", vertical="_+v"))
	cmd = add_opt(d, cmd, "SC", [:SC :col_axes],
	              (top=("t", nothing, 1), bott=("b", nothing, 1), bottom=("b", nothing, 1), label="+l", grid=("+w", add_opt_pen)))
	cmd = add_opt(d, cmd, "SR", [:SR :row_axes],
	              (left=("l", nothing, 1), right=("r", nothing, 1), label="+l", parallel="_+p", row_title="_+t", top_row_title="_+tc", grid=("+w", add_opt_pen)))
	opt_C = add_opt(d, "", "", [:C :clearance],
				  (left=(" -Cw", arg2str), right=(" -Ce", arg2str), bott=(" -Cs", arg2str), bottom=(" -Cs", arg2str), top=(" -Cn", arg2str)))
	cmd = add_opt(d, cmd, "Fs", [:Fs :panels_size :panel_size :panel_sizes])
	cmd = add_opt(d, cmd, "Ff", [:Ff :splot_size])

	if ((val = find_in_dict(d, [:F :dims :dimensions :size :sizes], false)[1]) !== nothing || show_kwargs[1])
		if (isa(val, NamedTuple) && haskey(nt2dict(val), :width))	# Preferred way
			cmd *= " -F" * helper_sub_F(val)		# VAL = (width=x, height=x, fwidth=(...), fheight=(...))
			del_from_dict(d, [:F :dims :dimensions :size :sizes])
		else
			cmd = add_opt(d, cmd, "F", [:F :dims :dimensions :size :sizes],
			              (panels=("-s", helper_sub_F, 1), size=("-f", helper_sub_F, 1), frac=("+f", helper_sub_F), fractions=("+f", helper_sub_F), clearance=("+c", arg2str), outine=("+p", add_opt_pen), fill=("+g", add_opt_fill), divlines=("+w", add_opt_pen)))
		end
	end

	do_set = false;		do_show = false
	if (fim !== nothing)
		t = lowercase(string(fim))
		if     (t == "end" || t == "stop")  stop = true
		elseif (t == "show")  stop = true;  do_show = true
		elseif (t == "set")   do_set = true
		end
	else
		IamModern[1] = false;	IamSubplot[1] = false			# Make sure we always start a clean session
		if (!stop && length(kwargs) == 0)  stop = true  end		# To account for the subplot() call case
	end
	# ------------------------------ End parsing inputs --------------------------------

	if (!stop && !do_set)
		((val_ = find_in_dict(d, [:grid])[1]) === nothing) && error("SUBPLOT: 'grid' keyword is mandatory")
		cmd = arg2str(val_, 'x') * " " * cmd * opt_C			# Also add the eventual global -C clearance option
		(dbg_print_cmd(d, cmd) !== nothing) && return cmd		# Vd=2 cause this return

		if (!IamModern[1])			# If we are not in modern mode, issue a gmt("begin") first
			# Default name (GMTplot.ps) is set in gmt_main()
			fname = ((val_ = find_in_dict(d, [:fmt])[1]) !== nothing) ? string("GMTplot ", val_) : ""
			if ((val_ = find_in_dict(d, [:figname :name :savefig])[1]) !== nothing)
				fname = get_format(string(val_), nothing, d)		# Get the fig name and format.
			end
			gmt("begin " * fname)
		end
		gmt("subplot begin " * cmd);
		IamSubplot[1] = true
	elseif (do_set)
		(!IamSubplot[1]) && error("Cannot call subplot(set, ...) before setting dimensions")
		lix, pane = parse_c(d, cmd)
		cmd = pane * cmd				# Here we don't want the "-c" part
		cmd = add_opt(d, cmd, 'A', [:fixedlabel]) * opt_C			# Also add the eventual this panel -C clearance option
		if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end		# Vd=2 cause this return
		gmt("subplot set " * cmd)
	else
		show = (do_show || haskey(d, :show)) ? " show" : ""
		try
			gmt("subplot end");		gmt("end" * show);		catch
		end
		IamModern[1] = false;		IamSubplot[1] = false
	end
	return nothing
end

# --------------------------------------------------------------------------
function helper_sub_F(arg, dumb=nothing)::String
	# dims=(1,2)
	# dims=(panels=(1,2), frac=((2,3),(3,4,5)))
	# dims=(width=xx, height=yy, fwidth=(), fheight=(), fill=:red, outline=(3,:red))
	out = ""
	if (isa(arg, String))
		out = arg2str(arg)
	elseif (isa(arg, NamedTuple) || isa(arg, Dict) || isa(arg, Tuple{Tuple, Tuple}) || isa(arg, Tuple{Tuple, Number}))
		d = mura_arg(arg)
		if ((val = find_in_dict(d, [:panels])[1]) !== nothing)
			if (isa(val, Tuple{Tuple, Tuple}))		# ex: dims=(panels=((2,4),(2.5,5,1.25)),)
				out *= arg2str(val[1], ',') * '/' * arg2str(val[2], ',')
			else
				out = arg2str(val)
			end
		end
		if ((val = find_in_dict(d, [:frac :fractions])[1]) !== nothing)		# ex: dims=(frac=((2,3),(3,4,5)))
			!isa(val, Tuple{Tuple, Tuple}) && error("'frac' option must be a tuple(tuple, tuple)")
			out *= arg2str(val[1], ',') * '/' * arg2str(val[2], ',')
		end
		if (haskey(d, :width))
			out *= string(d[:width], '/')
			out = (!haskey(d, :height)) ? string(out, d[:width]) : string(out, d[:height])
		end
		if (haskey(d, :fwidth))
			out *= "+f" * arg2str(d[:fwidth], ',')
			(!haskey(d, :fheight)) && error("SUBPLOT: when using 'fwidth' must also set 'fheight'")
			out *= '/' * arg2str(d[:fheight], ',')
		end
		if ((val = find_in_dict(d, [:fill], false)[1]) !== nothing)  out *= "+g" * add_opt_fill(val)  end
		if ((val = find_in_dict(d, [:clearance], false)[1]) !== nothing)  out *= "+c" * arg2str(val)  end
		if (haskey(d, :outline))   out *= "+p" * add_opt_pen(d, [:outline])  end
		if (haskey(d, :divlines))  out *= "+w" * add_opt_pen(d, [:divlines]) end
	elseif (isa(arg, Tuple))		# Hopefully only for the "dims=(panels=(xsize, ysize),)"
		out = arg2str(arg)
	end
	(out == "") && error("SUBPLOT: garbage in DIMS option")
	return out
end

# --------------------------------------------------------------------------
function mura_arg(arg)::Dict
	# Barrier function to contain a possible type instability
	if (isa(arg, Tuple{Tuple, Number}))  arg = (arg[1], (arg[2],))  end	# This looks terribly type instable
	# Need first case because for example dims=(panels=((2,4),(2.5,5,1.25)),) shows up here only as
	# arg = ((2, 4), (2.5, 5, 1.25)) because this function was called from within add_opt()
	if (isa(arg, Tuple{Tuple, Tuple}))  d = Dict(:panels => arg)
	else                                d = (isa(arg, NamedTuple)) ? nt2dict(arg) : arg
	end
	d
end


text:
using DiffEqBase, DiffEqJump, OrdinaryDiffEq, StochasticDiffEq, Test
using StableRNGs
rng = StableRNG(12345)

a = ExtendedJumpArray(rand(rng,3),rand(rng,2))
b = ExtendedJumpArray(rand(rng,3),rand(rng,2))

a.=b

@test a.u == b.u
@test a.jump_u == b.jump_u
@test a == b

c = rand(rng,5)
d = 2.0

a .+ d
a .= b .+ d
a .+ c .+ d
a .= b .+ c .+ d

rate = (u,p,t) -> u[1]
affect! = (integrator) -> (integrator.u[1] = integrator.u[1]/2)
jump = VariableRateJump(rate,affect!,interp_points=1000)
jump2 = deepcopy(jump)

f = function (du,u,p,t)
  du[1] = u[1]
end

prob = ODEProblem(f,[0.2],(0.0,10.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2; rng=rng)

integrator = init(jump_prob,Tsit5())

sol = solve(jump_prob,Tsit5())
sol = solve(jump_prob,Rosenbrock23(autodiff=false))
sol = solve(jump_prob,Rosenbrock23())

# @show sol[end]
# display(sol[end])

@test maximum([sol[i][2] for i in 1:length(sol)]) <= 1e-12
@test maximum([sol[i][3] for i in 1:length(sol)]) <= 1e-12

g = function (du,u,p,t)
  du[1] = u[1]
end

prob = SDEProblem(f,g,[0.2],(0.0,10.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2; rng=rng)

sol = solve(jump_prob,SRIW1())

@test maximum([sol[i][2] for i in 1:length(sol)]) <= 1e-12
@test maximum([sol[i][3] for i in 1:length(sol)]) <= 1e-12

function ff(du,u,p,t)
    if p == 0
        du .= 1.01u
    else
        du .= 2.01u
    end
end

function gg(du,u,p,t)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
end

rate_switch(u,p,t) = u[1]*1.0

function affect_switch!(integrator)
    integrator.p = 1
end

jump_switch = VariableRateJump(rate_switch,affect_switch!)

prob = SDEProblem(ff,gg,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,2))
jump_prob = JumpProblem(prob, Direct(), jump_switch; rng=rng)
solve(jump_prob, SRA1(), dt = 1.0)

## Some integration tests

function f2(du,u,p,t)
  du[1] = u[1]
end

prob = ODEProblem(f2,[0.2],(0.0,10.0))
rate2(u,p,t) = 2
affect2!(integrator) = (integrator.u[1] = integrator.u[1]/2)
jump = ConstantRateJump(rate2,affect2!)
jump_prob = JumpProblem(prob,Direct(),jump; rng=rng)
sol = solve(jump_prob,Tsit5())
sol(4.0)
sol[4]

rate2(u,p,t) = u[1]
affect2!(integrator) = (integrator.u[1] = integrator.u[1]/2)
jump = VariableRateJump(rate2,affect2!)
jump2 = deepcopy(jump)
jump_prob = JumpProblem(prob,Direct(),jump,jump2; rng=rng)
sol = solve(jump_prob,Tsit5())
sol(4.0)
sol[4]

function g2(du,u,p,t)
  du[1] = u[1]
end

prob = SDEProblem(f2,g2,[0.2],(0.0,10.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2; rng=rng)
sol = solve(jump_prob,SRIW1())
sol(4.0)
sol[4]

function f3(du, u, p, t)
    du .= u
end

prob = ODEProblem(f3, [1.0 2.0; 3.0 4.0], (0.0, 1.0))
rate3(u,p,t) = u[1] + u[2]
affect3!(integrator) = (integrator.u[1] = 0.25; integrator.u[2] = 0.5;
                        integrator.u[3] = 0.75; integrator.u[4] = 1)
jump = VariableRateJump(rate3, affect3!)
jump_prob = JumpProblem(prob, Direct(), jump; rng=rng)
sol = solve(jump_prob,Tsit5())


text:
mutable struct WordVectors{S<:AbstractString, T<:Real, H<:Integer}
    vocab::Vector{S} # vocabulary
    vectors::Array{T, 2} # the vectors computed from word2vec
    vocab_hash::Dict{S, H}
end

function WordVectors(vocab::AbstractArray{S,1},
                     vectors::AbstractArray{T,2}) where {S<: AbstractString, T<:Real}
    length(vocab) == size(vectors, 2) ||
        throw(DimensionMismatch("Dimension of vocab and vectors are inconsistent."))
    vocab_hash = Dict{S, Int}()
    for (i, word) in enumerate(vocab)
        vocab_hash[word] = i
    end
    WordVectors(vocab, vectors, vocab_hash)
end

function Base.show(io::IO, wv::WordVectors{S,T}) where {S,T}
    len_vecs, num_words = size(wv.vectors)
    print(io, "WordVectors $(num_words) words, $(len_vecs)-element $(T) vectors")
end

"""
    vocabulary(wv)

Return the vocabulary as a vector of words of the WordVectors `wv`.
"""
vocabulary(wv::WordVectors) = wv.vocab

"""
    in_vocabulary(wv, word)

Return `true` if `word` is part of the vocabulary of the WordVector `wv` and
`false` otherwise.
"""
in_vocabulary(wv::WordVectors, word::AbstractString) = word in wv.vocab

"""
    size(wv)

Return the word vector length and the number of words as a tuple.
"""
size(wv::WordVectors) = size(wv.vectors)

"""
    index(wv, word)

Return the index of `word` from the WordVectors `wv`.
"""
index(wv::WordVectors, word) = wv.vocab_hash[word]

"""
    get_vector(wv, word)

Return the vector representation of `word` from the WordVectors `wv`.
"""
get_vector(wv::WordVectors, word) =
      (idx = wv.vocab_hash[word]; wv.vectors[:,idx])

"""
    cosine(wv, word, n=10)

Return the position of `n` (by default `n = 10`) neighbors of `word` and their
cosine similarities.
"""
function cosine(wv::WordVectors, word, n=10)
    metrics = wv.vectors'*get_vector(wv, word)
    topn_positions = sortperm(metrics[:], rev = true)[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    similarity(wv, word1, word2)

Return the cosine similarity value between two words `word1` and `word2`.
"""
function similarity(wv::WordVectors, word1, word2)
    return get_vector(wv, word1)'*get_vector(wv, word2)
end

"""
    cosine_similar_words(wv, word, n=10)

Return the top `n` (by default `n = 10`) most similar words to `word`
from the WordVectors `wv`.
"""
function cosine_similar_words(wv::WordVectors, word, n=10)
    indx, metr = cosine(wv, word, n)
    return vocabulary(wv)[indx]
end

"""
    analogy(wv, pos, neg, n=5)

Compute the analogy similarity between two lists of words. The positions
and the similarity values of the top `n` similar words will be returned.
For example,
`king - man + woman = queen` will be
`pos=[\"king\", \"woman\"], neg=[\"man\"]`.
"""
function analogy(wv::WordVectors{S,T,H}, pos::AbstractArray, neg::AbstractArray, n= 5
                ) where {S<:AbstractString, T<:Real, H<:Integer}
    m, n_vocab = size(wv)
    n_pos = length(pos)
    n_neg = length(neg)
    anal_vecs = Matrix{T}(undef, m, n_pos + n_neg)

    for (i, word) in enumerate(pos)
        anal_vecs[:,i] = get_vector(wv, word)
    end
    for (i, word) in enumerate(neg)
        anal_vecs[:,i+n_pos] = -get_vector(wv, word)
    end
    mean_vec = mean(anal_vecs, dims=2)
    metrics = wv.vectors'*mean_vec
    top_positions = sortperm(metrics[:], rev = true)[1:n+n_pos+n_neg]
    for word in [pos;neg]
        idx = index(wv, word)
        loc = findfirst(x->x==idx, top_positions)
        if loc != nothing
            splice!(top_positions, loc)
        end
    end
    topn_positions = top_positions[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    analogy_words(wv, pos, neg, n=5)

Return the top `n` words computed by analogy similarity between
positive words `pos` and negaive words `neg`. from the WordVectors `wv`.
"""
function analogy_words(wv::WordVectors, pos, neg, n=5)
    indx, metr = analogy(wv, pos, neg, n)
    return vocabulary(wv)[indx]
end

"""
    wordvectors(filename [,type=Float64][; kind=:text, skip=false, normalize=true, fasttext=false])

Generate a WordVectors type object from a file.

# Arguments
  * `filename::AbstractString` the embeddings file name
  * `type::Type` type of the embedding vector elements; default `Float64`

# Keyword arguments
  * `kind::Symbol` specifies whether the embeddings file is textual (`:text`)
or binary (`:binary`); default `:text`
  * `skip::Bool` in binary embeddings files specifies whether the newline
byte is missing or not (use `true` for Google pre-trained models); default `false`
  * `normalize:Bool` specifies whether to normalize the embedding vectors
i.e. return unit vectors; default `true`
  * `delim:Char` delimiter; default `' '`
  * `fasttext::Bool` setting to `true` fixes a bug that occurs when loading
FastText generated text embeddings by not stripping the line before parsing.
The option is ignored for the binary format; default `false`
"""
function wordvectors(filename::AbstractString,
                     ::Type{T};
                     kind::Symbol=:text,
                     skip::Bool=false,
                     normalize::Bool=true,
                     delim::Char=' ',
                     fasttext::Bool=false) where T <: Real
    if kind == :binary
        return _from_binary(T, filename, skip, normalize, delim)
    elseif kind == :text
        return _from_text(T, filename, normalize, delim, fasttext)
    else
        throw(ArgumentError("Unknown embedding file kind $(kind)"))
    end
end

wordvectors(filename::AbstractString;
            kind::Symbol=:text,
            skip::Bool=false,
            normalize::Bool=true,
            delim::Char=' ',
            fasttext::Bool=false) =
    wordvectors(filename, Float64, kind=kind,
                normalize=normalize, delim=delim, fasttext=fasttext)

# Generate a WordVectors object from binary file
function _from_binary(::Type{T},
                      filename::AbstractString,
                      skip::Bool=true,
                      normalize::Bool=true,
                      delim::Char=' ') where T<:Real
    sb = ifelse(skip, 0, 1)
    open(filename) do f
        header = strip(readline(f))
        vocab_size,vector_size = map(x -> parse(Int, x), split(header, delim))
        vocab = Vector{String}(undef, vocab_size)
        vectors = zeros(T, vector_size, vocab_size)
        binary_length = sizeof(Float32) * vector_size
        for i in 1:vocab_size
            vocab[i] = strip(readuntil(f, delim))
            vector = reinterpret(Float32, read(f, binary_length))
            if normalize
                vector = vector ./ norm(vector)  # unit vector
            end
            vectors[:,i] = T.(vector)
            read(f, sb) # new line
        end
        return WordVectors(vocab, vectors)
    end
end

# Generate a WordVectors object from text file
function _from_text(::Type{T},
                    filename::AbstractString,
                    normalize::Bool=true,
                    delim::Char=' ',
                    fasttext::Bool=false) where T<:Real
    open(filename) do f
        header = strip(readline(f))
        vocab_size,vector_size = map(x -> parse(Int, x), split(header, delim))
        vocab = Vector{String}(undef, vocab_size)
        vectors = Matrix{T}(undef, vector_size, vocab_size)
        for (i, line) in enumerate(readlines(f))
            !fasttext && (line = strip(line))
            parts = split(line, delim)
            word = parts[1]
            vector = map(x-> parse(T, x), parts[2:end])
            vocab[i] = word
            if normalize
                vector = vector ./ norm(vector)  # unit vector
            end
            vectors[:, i] = vector
        end
       return WordVectors(vocab, vectors)
    end
end


text:
#=
 Copyright 2015-2016 Oliver Heimlich

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, see <http://www.gnu.org/licenses/>.
=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
# setrounding(Interval, :correct)
@format full

@testset "minimal.absRevBin_test" begin
    @test abs_rev(∅, entireinterval(Float64))[2] == ∅
    @test abs_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test abs_rev(Interval(0.0, 1.0), Interval(7.0, 9.0))[2] == ∅
    @test abs_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test abs_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2] == ∅
    @test abs_rev(Interval(1.0, 1.0), entireinterval(Float64))[2] == Interval(-1.0, 1.0)
    @test abs_rev(Interval(0.0, 0.0), entireinterval(Float64))[2] == Interval(0.0, 0.0)
    @test abs_rev(Interval(-1.0, -1.0), entireinterval(Float64))[2] == ∅
    @test abs_rev(Interval(0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023), entireinterval(Float64))[2] == Interval(-0x1.fffffffffffffp1023, 0x1.fffffffffffffp1023)
    @test abs_rev(Interval(0x1p-1022, 0x1p-1022), entireinterval(Float64))[2] == Interval(-0x1p-1022, 0x1p-1022)
    @test abs_rev(Interval(-0x1p-1022, -0x1p-1022), entireinterval(Float64))[2] == ∅
    @test abs_rev(Interval(-0x1.fffffffffffffp1023, -0x1.fffffffffffffp1023), entireinterval(Float64))[2] == ∅
    @test abs_rev(Interval(1.0, 2.0), entireinterval(Float64))[2] == Interval(-2.0, 2.0)
    @test abs_rev(Interval(1.0, 2.0), Interval(0.0, 2.0))[2] == Interval(1.0, 2.0)
    @test abs_rev(Interval(0.0, 1.0), Interval(-0.5, 2.0))[2] == Interval(-0.5, 1.0)
    @test abs_rev(Interval(-1.0, 1.0), entireinterval(Float64))[2] == Interval(-1.0, 1.0)
    @test abs_rev(Interval(-1.0, 0.0), entireinterval(Float64))[2] == Interval(0.0, 0.0)
    @test abs_rev(Interval(0.0, Inf), entireinterval(Float64))[2] == entireinterval(Float64)
    @test abs_rev(entireinterval(Float64), entireinterval(Float64))[2] == entireinterval(Float64)
    @test abs_rev(Interval(-Inf, 0.0), entireinterval(Float64))[2] == Interval(0.0, 0.0)
    @test abs_rev(Interval(1.0, Inf), Interval(-Inf, 0.0))[2] == Interval(-Inf, -1.0)
    @test abs_rev(Interval(-1.0, Inf), entireinterval(Float64))[2] == entireinterval(Float64)
    @test abs_rev(Interval(-Inf, -1.0), entireinterval(Float64))[2] == ∅
    @test abs_rev(Interval(-Inf, 1.0), entireinterval(Float64))[2] == Interval(-1.0, 1.0)
end
# FactCheck.exitstatus()


text:
############################################################################################
"""
$(SIGNATURES)
Sample from constraint if 'prior'.

# Examples
```julia
```

"""
function sample_constraint(_rng::Random.AbstractRNG, prior, val)
    return val
end
function sample_constraint(_rng::Random.AbstractRNG, prior::Distributions.Distribution, val)
    return rand(_rng, prior)
end
function sample_constraint(_rng::Random.AbstractRNG, priorᵥ::Vector{<:Distributions.Distribution}, val::AbstractVector)
    return rand.(_rng, priorᵥ)
end
function sample_constraint(_rng::Random.AbstractRNG, priorᵥ::AbstractArray, val::AbstractArray)
    return map(iter -> sample_constraint(_rng, priorᵥ[iter], val[iter]), eachindex(priorᵥ))
end
function sample_constraint(_rng::Random.AbstractRNG, priorᵥ::NamedTuple{names}, val::NamedTuple) where {names}
    return NamedTuple{names}(Tuple(map(iter -> sample_constraint(_rng, priorᵥ[iter], val[iter]), names)))
end

############################################################################################
"""
$(SIGNATURES)
Evaluate Log density of 'prior' at 'θ'.

# Examples
```julia
```

"""
function log_prior(prior, θ::T) where {T}
    return 0.0
end
function log_prior(prior::Distributions.Distribution, θ::T) where {T}
    return Distributions.logpdf(prior, θ)
end
function log_prior(priorᵥ::Vector{<:Distributions.Distribution}, θ::AbstractVector)
    return sum(Distributions.logpdf.(priorᵥ, θ))
end
function log_prior(priorᵥ::AbstractArray, θ::AbstractArray) where {T}
    return sum(map(log_prior, priorᵥ, θ))
end
function log_prior(priorᵥ::A, θ::B) where {A<:NamedTuple,B<:NamedTuple}
    return sum(map(log_prior, priorᵥ, θ))
end

############################################################################################
"""
$(SIGNATURES)
Evaluate Log density and eventual Jacobian adjustments from transformation of 'prior' at 'θ'.

# Examples
```julia
```

"""
function log_prior_with_transform(prior, θ::T) where {T}
    return 0.0
end
function log_prior_with_transform(prior::Distributions.Distribution, θ::T) where {T}
    return Bijectors.logpdf_with_trans(prior, θ, true)
end
function log_prior_with_transform(
    priorᵥ::Vector{<:Distributions.Distribution}, θ::AbstractVector
)
    return sum(Bijectors.logpdf_with_trans.(priorᵥ, θ, true))
end
function log_prior_with_transform(priorᵥ::AbstractArray, θ::AbstractArray) where {T}
    return sum(map(log_prior_with_transform, priorᵥ, θ))
end
function log_prior_with_transform(priorᵥ::A, θ::B) where {A<:NamedTuple,B<:NamedTuple}
    return sum(map(log_prior_with_transform, priorᵥ, θ))
end

############################################################################################
#export
export sample_constraint, log_prior, log_prior_with_transform


text:
include("ModuleA.jl")

importall MyModule

mytype = MyType(4)

@show foo(mytype)
@show MyType
#@show bar

@show MyModule.foo(mytype)
@show MyModule.MyType
@show MyModule.bar(2)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function characterize_S2(	particle_type::String,
							R::Array{Float64, 2},
							Lx::Float64,
							Ly::Float64,
							Lz::Float64,
							X::Array{Float64, 1},
							Y::Array{Float64, 1},
							Z::Array{Float64, 1},
							Q0::Array{Float64, 1},
							Q1::Array{Float64, 1},
							Q2::Array{Float64, 1},
							Q3::Array{Float64, 1},
							A11::Array{Float64, 1},
							A12::Array{Float64, 1},
							A13::Array{Float64, 1},
							A21::Array{Float64, 1},
							A22::Array{Float64, 1},
							A23::Array{Float64, 1},
							A31::Array{Float64, 1},
							A32::Array{Float64, 1},
							A33::Array{Float64, 1},
							number_of_samples::Int64,
							d2::Array{Float64, 1},
							cell_lists::Array{Array{Int64, 1}, 3})

	# Number of cells.
	(number_of_cells_x::Int64, number_of_cells_y::Int64, number_of_cells_z::Int64) = size(cell_lists)

	current_particle_in_cell::Int64 = 0
	is_void::Bool = true

	number_of_points::Int64 = length(d2)
	x::Array{Float64, 1} = zeros(number_of_points)
	y::Array{Float64, 1} = zeros(number_of_points)
	z::Array{Float64, 1} = zeros(number_of_points)
	S2::Array{Float64, 1} = zeros(number_of_points)
	indicator::Array{Bool, 1} = falses(number_of_points)

	vx::Float64 = 0.0
	vy::Float64 = 0.0
	vz::Float64 = 0.0

	current_cell_x::Int64 = 0
	current_cell_y::Int64 = 0
	current_cell_z::Int64 = 0
	number_of_particles_current_cell::Int64 = 0

	for current_sample = 1:number_of_samples
		#if mod(current_sample, 1000) == 0
		#	println(current_sample)
		#end

		# Pick random intial sample position.
		x0 = Lx * rand()
		y0 = Ly * rand()
		z0 = Lz * rand()

		# Pick random direction vector and normalize.
		vx = randn()
		vy = randn()
		vz = randn()
		magnitude = sqrt(vx * vx + vy * vy + vz * vz)
		vx /= magnitude
		vy /= magnitude
		vz /= magnitude

		# Create vectors of sample coordinates in relative coordinates.
		x = x0 + vx * d2
		y = y0 + vy * d2
		z = z0 + vz * d2

		# Convert to absolute coordinates (w.r.t. to Lx, Ly, Lz)
		for current_point = 1:number_of_points
			x[current_point] = position_mod(x[current_point], Lx)
			y[current_point] = position_mod(y[current_point], Ly)
			z[current_point] = position_mod(z[current_point], Lz)
		end

		# Check whether void or solid for each sample point.
		for current_point = 1:number_of_points

			current_cell_x = convert(Int64, ceil(x[current_point] / Lx * convert(Float64, number_of_cells_x)))
			current_cell_y = convert(Int64, ceil(y[current_point] / Ly * convert(Float64, number_of_cells_y)))
			current_cell_z = convert(Int64, ceil(z[current_point] / Lz * convert(Float64, number_of_cells_z)))
			number_of_particles_current_cell = length(cell_lists[current_cell_x, current_cell_y, current_cell_z])
			current_cell_list = cell_lists[current_cell_x, current_cell_y, current_cell_z]

			is_void = true
			current_particle_in_cell = 0
			while current_particle_in_cell < number_of_particles_current_cell && is_void
				current_particle_in_cell += 1

				vx = signed_distance_mod(x[current_point], X[current_cell_list[current_particle_in_cell]], Lx)
				vy = signed_distance_mod(y[current_point], Y[current_cell_list[current_particle_in_cell]], Ly)
				vz = signed_distance_mod(z[current_point], Z[current_cell_list[current_particle_in_cell]], Lz)

				if particle_type == "sphere"
					if vx^2 + vy^2 + vz^2 <= R[current_cell_list[current_particle_in_cell], 1]^2
						is_void = false
					end
				elseif particle_type == "ellipse"
					# Not supported for S2.
				elseif particle_type == "ellipsoid"
					if vx * (A11[current_cell_list[current_particle_in_cell]] * vx + A12[current_cell_list[current_particle_in_cell]] * vy + A13[current_cell_list[current_particle_in_cell]] * vz) + vy * (A21[current_cell_list[current_particle_in_cell]] * vx + A22[current_cell_list[current_particle_in_cell]] * vy + A23[current_cell_list[current_particle_in_cell]] * vz) + vz * (A31[current_cell_list[current_particle_in_cell]] * vx + A32[current_cell_list[current_particle_in_cell]] * vy + A33[current_cell_list[current_particle_in_cell]] * vz) <= 1.0
						is_void = false
					end
				elseif particle_type == "cuboid"
					(vx, vy, vz) = (A11[current_cell_list[current_particle_in_cell]] * vx + A12[current_cell_list[current_particle_in_cell]] * vy + A13[current_cell_list[current_particle_in_cell]] * vz,
									A21[current_cell_list[current_particle_in_cell]] * vx + A22[current_cell_list[current_particle_in_cell]] * vy + A23[current_cell_list[current_particle_in_cell]] * vz,
									A31[current_cell_list[current_particle_in_cell]] * vx + A32[current_cell_list[current_particle_in_cell]] * vy + A33[current_cell_list[current_particle_in_cell]] * vz)

					if abs(vx) <= R[current_cell_list[current_particle_in_cell], 1] && abs(vy) <= R[current_cell_list[current_particle_in_cell], 2] && abs(vz) <= R[current_cell_list[current_particle_in_cell], 3]
						is_void = false
					end
				elseif particle_type == "superellipsoid"
					(vx, vy, vz) = (A11[current_cell_list[current_particle_in_cell]] * vx + A12[current_cell_list[current_particle_in_cell]] * vy + A13[current_cell_list[current_particle_in_cell]] * vz,
									A21[current_cell_list[current_particle_in_cell]] * vx + A22[current_cell_list[current_particle_in_cell]] * vy + A23[current_cell_list[current_particle_in_cell]] * vz,
									A31[current_cell_list[current_particle_in_cell]] * vx + A32[current_cell_list[current_particle_in_cell]] * vy + A33[current_cell_list[current_particle_in_cell]] * vz)

					if (abs(vx)/R[current_cell_list[current_particle_in_cell], 1])^R[current_cell_list[current_particle_in_cell], 4] +
						(abs(vy)/R[current_cell_list[current_particle_in_cell], 2])^R[current_cell_list[current_particle_in_cell], 4] +
						(abs(vz)/R[current_cell_list[current_particle_in_cell], 3])^R[current_cell_list[current_particle_in_cell], 4] <= 1.0
						is_void = false
					end
				end
			end

			if is_void
				indicator[current_point] = true
			else
				indicator[current_point] = false
			end
		end

		if indicator[1] # If the origin point is in void, do something, otherwise don't.
			for current_point = 1:number_of_points
				if indicator[current_point]
					S2[current_point] += 1.0
				end
			end
		end
	end

	return S2
end


text:
#
# Compares the performance of Dict based parameters versus type based parameters
# Dict based parameters is implemented in DictCurve type.
# Type based parameters is implemented in SvenssonCurve type.
# SvenssonCurveNullable is the same as SvenssonCurve but has nullable parameters.
# Array based parameters is implemented in ArrayCurve type.
#
# Type and Array based parameters are equivalent in performance.
# Dict based parameters doubles the time execution compared to the other methods.
#

using Base.Test

type DictCurve
	params::Dict{Symbol, Float64}
end

type SvenssonCurve
	beta1::Float64
	beta2::Float64
	beta3::Float64
	beta4::Float64
	lambda1::Float64
	lambda2::Float64
end

type SvenssonCurveNullable
	beta1::Nullable{Float64}
	beta2::Nullable{Float64}
	beta3::Nullable{Float64}
	beta4::Nullable{Float64}
	lambda1::Nullable{Float64}
	lambda2::Nullable{Float64}
end

type ArrayCurve
	params::Array{Float64}
end

function zero_rate(curve::ArrayCurve, t::Float64)

	# beta1 = param[1]
	# beta2 = param[2]
	# beta3 = param[3]
	# beta4 = param[4]
	# lambda1 = param[5]
	# lambda2 = param[6]

	_exp_lambda1_t_ = exp(-curve.params[5]*t)
	_exp_lambda2_t_ = exp(-curve.params[6]*t)
	F_beta2 = (1.0 - _exp_lambda1_t_) / (curve.params[5]*t)

	return curve.params[1] + curve.params[2]*F_beta2 + curve.params[3]*(F_beta2 - _exp_lambda1_t_) +
			curve.params[4]*( (1.0 - _exp_lambda2_t_)/(curve.params[6]*t) - _exp_lambda2_t_)
end

function zero_rate(curve::SvenssonCurve, t::Float64)

	# beta1 = param[1]
	# beta2 = param[2]
	# beta3 = param[3]
	# beta4 = param[4]
	# lambda1 = param[5]
	# lambda2 = param[6]

	_exp_lambda1_t_ = exp(-curve.lambda1*t)
	_exp_lambda2_t_ = exp(-curve.lambda2*t)
	F_beta2 = (1.0 - _exp_lambda1_t_) / (curve.lambda1*t)

	return curve.beta1 + curve.beta2*F_beta2 + curve.beta3*(F_beta2 - _exp_lambda1_t_) +
			curve.beta4*( (1.0 - _exp_lambda2_t_)/(curve.lambda2*t) - _exp_lambda2_t_)
end

function zero_rate(curve::SvenssonCurveNullable, t::Float64)

	# beta1 = param[1]
	# beta2 = param[2]
	# beta3 = param[3]
	# beta4 = param[4]
	# lambda1 = param[5]
	# lambda2 = param[6]

	_exp_lambda1_t_ = exp(-get(curve.lambda1)*t)
	_exp_lambda2_t_ = exp(-get(curve.lambda2)*t)
	F_beta2 = (1.0 - _exp_lambda1_t_) / (get(curve.lambda1)*t)

	return get(curve.beta1) + get(curve.beta2)*F_beta2 + get(curve.beta3)*(F_beta2 - _exp_lambda1_t_) +
			get(curve.beta4)*( (1.0 - _exp_lambda2_t_)/(get(curve.lambda2)*t) - _exp_lambda2_t_)
end

function zero_rate(curve::DictCurve, t::Float64)

	# beta1 = param[1]
	# beta2 = param[2]
	# beta3 = param[3]
	# beta4 = param[4]
	# lambda1 = param[5]
	# lambda2 = param[6]

	_exp_lambda1_t_ = exp(-curve.params[:lambda1]*t)
	_exp_lambda2_t_ = exp(-curve.params[:lambda2]*t)
	F_beta2 = (1.0 - _exp_lambda1_t_) / (curve.params[:lambda1]*t)

	return curve.params[:beta1] + curve.params[:beta2]*F_beta2 + curve.params[:beta3]*(F_beta2 - _exp_lambda1_t_) +
			curve.params[:beta4]*( (1.0 - _exp_lambda2_t_)/(curve.params[:lambda2]*t) - _exp_lambda2_t_)
end

params = [0.1, 0.2, 0.3, 0.4, 0.5, 0.8]

params_dict = Dict(:beta1=>params[1], :beta2=>params[2], :beta3=>params[3], :beta4=>params[4],
				:lambda1=>params[5], :lambda2=>params[6])

d_curve = DictCurve(params_dict)
s_curve = SvenssonCurve(params...)
sn_curve = SvenssonCurveNullable(Nullable(params[1]), Nullable(params[2]), Nullable(params[3]), Nullable(params[4]), Nullable(params[5]), Nullable(params[6]))
a_curve = ArrayCurve(params)

mat_vec = linspace(0.1, 10.0, 1000000)

 # warmup
zero_rate(d_curve, mat_vec[1])
zero_rate(s_curve, mat_vec[1])
zero_rate(sn_curve, mat_vec[1])
zero_rate(a_curve, mat_vec[1])

println("Dict based parameters")
@time for t in mat_vec
	zero_rate(d_curve, t)
end

println("Type based parameters")
@time for t in mat_vec
	zero_rate(s_curve, t)
end

println("Type based nullable parameters")
@time for t in mat_vec
	zero_rate(sn_curve, t)
end

println("Array based parameters")
@time for t in mat_vec
	zero_rate(a_curve, t)
end


text:
const RNEAREST = Int32(4)

function add2!(x::ArbFloat{P}, y::ArbFloat{P}) where P
    rounddir = ccall(@libarb(arf_add), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), x, x, y, P, RNEAREST)
    return x
end

function sub2!(x::ArbFloat{P}, y::ArbFloat{P}) where P
    rounddir = ccall(@libarb(arf_sub), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), x, x, y, P, RNEAREST)
    return x
end

function mul2!(x::ArbFloat{P}, y::ArbFloat{P}) where P
    rounddir = ccall(@libarb(arf_mul_rnd_any), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), x, x, y, P, RNEAREST)
    return x
end

function div2!(x::ArbFloat{P}, y::ArbFloat{P}) where P
    rounddir = ccall(@libarb(arf_div), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), x, x, y, P, RNEAREST)
    return x
end

function add2!(x::ArbReal{P}, y::ArbReal{P}) where P
    ccall(@libarb(arb_add), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), x, x, y, P)
    return x
end

function sub2!(x::ArbReal{P}, y::ArbReal{P}) where P
    ccall(@libarb(arb_sub), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), x, x, y, P)
    return x
end

function mul2!(x::ArbReal{P}, y::ArbReal{P}) where P
    ccall(@libarb(arb_mul), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), x, x, y, P)
    return x
end

function div2!(x::ArbReal{P}, y::ArbReal{P}) where P
    ccall(@libarb(arb_div), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), x, x, y, P)
    return x
end

function add2!(x::ArbComplex{P}, y::ArbComplex{P}) where P
    ccall(@libarb(acb_add), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong), x, x, y, P)
    return x
end

function sub2!(x::ArbComplex{P}, y::ArbComplex{P}) where P
    ccall(@libarb(acb_sub), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong), x, x, y, P)
    return x
end

function mul2!(x::ArbComplex{P}, y::ArbComplex{P}) where P
    ccall(@libarb(acb_mul), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong), x, x, y, P)
    return x
end

function div2!(x::ArbComplex{P}, y::ArbComplex{P}) where P
    ccall(@libarb(acb_div), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong), x, x, y, P)
    return x
end

"""
    add2!(x, y)

x gets the sum of x and y

returns updated x
""" add2!

"""
    sub2!(x, y)

x gets the difference of x and y

returns updated x
""" sub2!

"""
    mul2!(x, y)

x gets the product of x and y

returns updated x
""" mul2!

"""
    div2!(x, y)

x gets the quotient of x and y

returns updated x
""" div2!


text:
function MOI.add_variable(model::Optimizer)
    model.gradient_cache = nothing
    vi = MOI.add_variable(model.optimizer)
    return vi
end

function MOI.add_variables(model::Optimizer, N::Int)
    model.gradient_cache = nothing
    return VI[MOI.add_variable(model) for i in 1:N]
end

function MOI.add_constraint(model::Optimizer, f::SUPPORTED_SCALAR_FUNCTIONS, s::SUPPORTED_SCALAR_SETS)
    model.gradient_cache = nothing
    return MOI.add_constraint(model.optimizer, f, s)
end

function MOI.add_constraint(model::Optimizer, vf::SUPPORTED_VECTOR_FUNCTIONS, s::SUPPORTED_VECTOR_SETS)
    model.gradient_cache = nothing
    return MOI.add_constraint(model.optimizer, vf, s)
end

function MOI.add_constraints(model::Optimizer, f::AbstractVector{F}, s::AbstractVector{S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    return CI{F, S}[MOI.add_constraint(model, f[i], s[i]) for i in eachindex(f)]
end

function MOI.set(model::Optimizer, attr::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES}, f::SUPPORTED_OBJECTIVES)
    model.gradient_cache = nothing
    MOI.set(model.optimizer, attr, f)
end

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.set(model::Optimizer, attr::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, sense)
end

function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute)
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, s)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, f)
end

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, s)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, f)
end

# `MOI.supports` methods

function MOI.supports(::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}})
    return true
end

function MOI.supports(::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}})
    return true
end

function MOI.supports(::Optimizer, ::MOI.AbstractModelAttribute)
    return true
end

function MOI.supports(::Optimizer, ::MOI.ObjectiveFunction)
    return false
end

function MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES})
    return true
end

MOI.supports_constraint(::Optimizer, ::Type{<: SUPPORTED_SCALAR_FUNCTIONS}, ::Type{<: SUPPORTED_SCALAR_SETS}) = true
MOI.supports_constraint(::Optimizer, ::Type{<: SUPPORTED_VECTOR_FUNCTIONS}, ::Type{<: SUPPORTED_VECTOR_SETS}) = true
function MOI.supports(model::Optimizer, attr::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.supports(model.optimizer, attr, CI{F,S})
end
function MOI.supports(model::Optimizer, attr::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.supports(model.optimizer, attr, CI{F,S})
end

MOI.get(model::Optimizer, attr::MOI.SolveTimeSec) = MOI.get(model.optimizer, attr)

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.optimizer)
    model.gradient_cache = nothing
    return
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
           model.gradient_cache === nothing
end

# now supports name too
MOI.supports_incremental_interface(::Optimizer) = true

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike)
    model.gradient_cache = nothing
    return MOIU.default_copy_to(model.optimizer, src)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    return MOI.get(model.optimizer, MOI.TerminationStatus())
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::VI, value::Float64)
    model.gradient_cache = nothing
    MOI.set(model.optimizer, MOI.VariablePrimalStart(), vi, value)
end

function MOI.supports(model::Optimizer, attr::MOI.AbstractVariableAttribute,
                      ::Type{MOI.VariableIndex})
    return MOI.supports(model.optimizer, attr, MOI.VariableIndex)
end

function MOI.set(model::Optimizer, attr::MOI.AbstractVariableAttribute, v::MOI.VariableIndex, value)
    MOI.set(model.optimizer, attr, v, value)
end

function MOI.get(model::Optimizer, attr::MOI.AbstractVariableAttribute, vi::MOI.VariableIndex)
    return MOI.get(model.optimizer, attr, vi)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    MOI.delete(model.optimizer, ci)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    MOI.delete(model.optimizer, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.is_valid(model::Optimizer, v::VI)
    return MOI.is_valid(model.optimizer, v::VI)
end

MOI.is_valid(model::Optimizer, con::CI) = MOI.is_valid(model.optimizer, con)

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintBasisStatus, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintBasisStatus(), ci)
end

# helper methods to check if a constraint contains a Variable
function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VariableIndex, S}) where {S <: SUPPORTED_SCALAR_SETS}
    return v == MOI.get(model, MOI.ConstraintFunction(), ci)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.ScalarAffineFunction{Float64}, S}) where {S <: SUPPORTED_SCALAR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.variable, func.terms)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorOfVariables, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return v in func.variables
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorAffineFunction{Float64}, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.scalar_term.variable, func.terms)
end

function MOI.delete(model::Optimizer, v::VI)
    model.gradient_cache = nothing
    MOI.delete(model.optimizer, v)
end

# for array deletion
function MOI.delete(model::Optimizer, indices::Vector{VI})
    model.gradient_cache = nothing
    for i in indices
        MOI.delete(model, i)
    end
end

function MOI.modify(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    chg::MOI.AbstractFunctionModification
)
    model.gradient_cache = nothing
    MOI.modify(
        model.optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        chg
    )
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    MOI.modify(model.optimizer, ci, chg)
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:MOI.VectorAffineFunction{Float64}, S <: SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    MOI.modify(model.optimizer, ci, chg)
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    return MOI.get(model.optimizer, MOI.VariableIndex, name)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintName, con::CI, name::String)
    MOI.set(model.optimizer, MOI.ConstraintName(), con, name)
end

function MOI.get(model::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    return MOI.get(model.optimizer, MOI.ConstraintIndex, name)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, con::CI)
    return MOI.get(model.optimizer, MOI.ConstraintName(), con)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{F,S})
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{VF, VS}}) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{VF,VS})
end

function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    MOI.set(model.optimizer, MOI.Name(), name)
end

function MOI.get(model::Optimizer, ::Type{CI{F, S}}, name::String) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, CI{F,S}, name)
end

function MOI.get(model::Optimizer, ::Type{CI{VF, VS}}, name::String) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, CI{VF,VS}, name)
end

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Real, Nothing})
    MOI.set(model.optimizer, MOI.TimeLimitSec(), value)
end
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model.optimizer, MOI.TimeLimitSec())
end

MOI.supports(model::Optimizer, ::MOI.Silent) = MOI.supports(model.optimizer, MOI.Silent())

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    MOI.set(model.optimizer, MOI.Silent(), value)
end

function MOI.get(model::Optimizer, ::MOI.Silent)
    return MOI.get(model.optimizer, MOI.Silent())
end

function MOI.get(model::Optimizer, ::MOI.SolverName)
    return MOI.get(model.optimizer, MOI.SolverName())
end


text:
cd(@__DIR__)

tests = readdir()
filter!(x -> x != "runtests.jl", tests)

res = map(tests) do t
    include(t)
end


text:
#-------------Max-Plus------------------------#

import Base.show,Base.zero,Base.one,Base.promote_rule,Base.+,Base.*,Base.convert
immutable MP  v::Float64 end
# promote_rule(MP,Real)=MP
promote_rule(::Type{MP},::Type{Int64})=MP
promote_rule(::Type{MP},::Type{Float64})=MP
*(a::MP,b::MP)= MP(a.v+b.v)
*(a::MP,b::Real)=MP(a.v+b)
*(a::Real,b::MP)=MP(a+b.v)
+(a::MP,b::MP)= MP(max(a.v,b.v))
+(a::MP,b::Real)=MP(max(a.v,b))
+(a::Real,b::MP)=MP(max(a,b.v))
show(io::IO,k::MP)=print(io,k.v)
zero(::Type{MP})=MP(typemin(Float64))
zero(::MP)=MP(typemin(Float64))
one(::Type{MP})= MP(zero(Float64))
convert(::Type{MP},x::Float64)=MP(x)
convert(::Type{MP},x::Int64)=MP(Float64(x))
convert(::Type{Float64},x::MP)=x.v
mparray(A::Array)=map(MP,A)
array(A::Array{MP})=map(x->x.v, A)
mpsparse(S::SparseMatrixCSC)=SparseMatrixCSC(S.m,S.n,S.colptr,S.rowval,mparray(S.nzval))
mpsparse(M::Array{Float64})=mpsparse(sparse(M))
mpsparse(M::Array{MP})=(k=0; R::Array{Int64,1}=[]; C::Array{Int64,1}=[]; V::Array{MP,1}=[];
    for i=1:size(M,1) for j=1:size(M,2)
     M[i,j]!=zero(MP) && (R=[R;i]; C=[C;j]; V=[V;M[i,j]])
    end end ;
    sparse(R,C,V) )

#------------Test-Max-Pus----- -----------------#

promote_type(MP,Float64)
promote_type(MP,Int64)
one(MP)
zero(MP)
ones(MP,2,2)
zeros(MP,3,4)
spzeros(MP,2,2)
speye(MP,2,3)
eye(MP,3,4)
MP(2.0)+MP(3)
MP(2.0)*MP(3.0)
MP(2.0)+3.0
MP(4.0)+2
[MP(2) 2.0]
[MP(2) 2]
MP[2 3 5]
MP[2.0 3.0]
MP[2 3;4 5]
MP[MP(4) 2]
a=[MP(2.0) 3.0;  1 2.0 ]
a+a
a*a
b=rand(2,2)
c=mparray(b)
a*b
d=array(c)
n=400
a = rand(n,n);  gc(); @time  a*a
a = mparray(rand(n,n));  gc(); @time a*a
c=rand(n,n); gc(); @time a*c
a=sparse(rand(n,n)); gc(); @time a*a
a=mpsparse(rand(n,n)); gc(); @time a*a

#----------------Ring-Operator------------------#

function ringmatmul(+, *, A::Matrix, B::Matrix)
    m, n, p = size(A,1), size(B,2), size(A,2)
    C = [A[i,1]*B[1,j] for i=1:m, j=1:n]
    for i=1:m, j=1:n, k=2:p
        C[i,j] += A[i,k]*B[k,j]
    end
    return C
end
*((+)::Function, (*)::Function) = (A,B)->ringmatmul(+,*,A,B)

#----------------Test-Ring-Operator--------------#

a = rand(n,n); b = rand(n,n); gc(); @time  *(max,+)(a,b)

max_stars_repo_name: Lecrapouille/BacASable

max_issues_repo_name: Lecrapouille/BacASable

max_forks_repo_name: Lecrapouille/BacASable


text:
# Exponential map via Rodriguez formula
# =====================================
function expmap(O)
  O = vec(O);
  dp = dot(O,O);
  ln = sqrt(dp);
  W = skew(O);
  if ln < 1E-15
    R = Matrix{Float64}(I, 3, 3);
  else
    R = I + (sin(ln)/ln)*W + ((1-cos(ln))/dp)*W*W;
  end
  return R;
end

# Skew-symmetric matrix out of 3-vector
# =====================================
function skew(O)
  [0.0 -O[3] O[2];
   O[3] 0.0 -O[1];
  -O[2] O[1] 0.0];
end

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
const ThreeAgentPos = SVector{6,Int}

const syma_ego = [[:up,], [:down,], [:left,], [:right,]]
const syma_adv = [[:up, :up], [:down, :up], [:left, :up], [:right, :up],
                [:up, :down], [:down, :down], [:left, :down], [:right, :down],
                [:up, :left], [:down, :left], [:left, :left], [:right, :left],
                [:up, :right], [:down, :right], [:left, :right], [:right, :right]]
const aind_ego = Dict(syma_ego[i] => i for i=1:length(syma_ego))
const aind_adv = Dict(syma_adv[i] => i for i=1:length(syma_adv))
const S2 = ThreeAgentPos
const A2 = Array{Symbol}

# Gridworld with adversary
@with_kw mutable struct DoubleAdversarialGridworldMDP <: MDP{S2, A2}
    size::Tuple{Int, Int} = (10,10)
    rewards::Dict{GWPos, Float64} = Dict()
    walls::Vector{GWPos} = []
    tprob::Float64 = 0.7
    discount::Float64 = 0.95
    agent_gets_action = :ego # :ego or :adversary
    ego_policy = (s, rng::AbstractRNG = Random.GLOBAL_RNG) -> rand(rng, syma_ego)
    adversary_policy = (s, rng::AbstractRNG = Random.GLOBAL_RNG) -> rand(rng, syma_adv)
    failure_penalty = 1.
end

valid_pos(mdp::DoubleAdversarialGridworldMDP, pos::GWPos) = !(pos in mdp.walls || any((pos .> mdp.size) .| (pos .< GWPos(1,1))))

function POMDPs.states(mdp::DoubleAdversarialGridworldMDP)
    lengths = (mdp.size[1], mdp.size[2], mdp.size[1], mdp.size[2], mdp.size[1], mdp.size[2])
    ss = S2[]
    for ijk in CartesianIndices(lengths)
        s = S2(ijk.I...)
        valid_pos(mdp, ego_pos(s)) && valid_pos(mdp, adversary1_pos(s)) && valid_pos(mdp, adversary2_pos(s))  &&  push!(ss, s)
    end
    push!(ss, S2(-1,-1,-1,-1,-1,-1))
    ss
end

function random_valid_pos(mdp::DoubleAdversarialGridworldMDP, rng::AbstractRNG = Random.GLOBAL_RNG, exclude = [], max_trials = 1000)
    trial = 0
    while trial < max_trials
        pos = GWPos(rand(rng, 1:mdp.size[1]), rand(rng, 1:mdp.size[2]))
        if valid_pos(mdp, pos) && !(haskey(mdp.rewards, pos) || pos in exclude)
            return pos
        end
        trial += 1
    end
 end

function POMDPs.initialstate(mdp::DoubleAdversarialGridworldMDP, rng::AbstractRNG = Random.GLOBAL_RNG)
    ego = random_valid_pos(mdp, rng)
    adversary1 = random_valid_pos(mdp, rng, [ego])
    adversary2 = random_valid_pos(mdp, rng, [ego, adversary1])
    Deterministic(S2(ego..., adversary1..., adversary2...))
end

POMDPs.actions(mdp::DoubleAdversarialGridworldMDP) = mdp.agent_gets_action == :ego ? syma_ego : syma_adv
POMDPs.actionindex(mdp::DoubleAdversarialGridworldMDP, a::A2) = mdp.agent_gets_action == :ego ? aind_ego[a] : aind_adv[a]

ego_pos(s::S2) = GWPos(s[1], s[2])
adversary1_pos(s::S2) = GWPos(s[3], s[4])
adversary2_pos(s::S2) = GWPos(s[5], s[6])
agents_overlap(s::S2) = (ego_pos(s) == adversary1_pos(s)) && (ego_pos(s) == adversary2_pos(s))

POMDPs.isterminal(mdp::DoubleAdversarialGridworldMDP, s::S2) = any(s .< 0)
POMDPs.discount(mdp::DoubleAdversarialGridworldMDP) = mdp.discount

# Returns a sample next state and reward
function POMDPs.gen(mdp::DoubleAdversarialGridworldMDP, s::S2, a::A2, rng::AbstractRNG = Random.GLOBAL_RNG)
    if haskey(mdp.rewards, ego_pos(s)) || agents_overlap(s) || isterminal(mdp, s)
        return (sp = S2(-1,-1,-1,-1,-1,-1), r = reward(mdp, s))
    else
        # If this MDP controls the agent then use the adversary policy for the adversary
        # Do the opposite if the adversary is being controlled by the MDP action
        if mdp.agent_gets_action == :ego
            rdir = (rand(rng) < mdp.tprob) ? dir[a[1]] : dir[rand(rng, syma[a[1] .!= syma])]
            new_ego = ego_pos(s) + rdir
            a_adv = mdp.adversary_policy(s, rng)
            new_adv1 = adversary1_pos(s) + dir[a_adv[1]]
            new_adv2 = adversary2_pos(s) + dir[a_adv[2]]
        else
            rdir1 = (rand(rng) < mdp.tprob) ? dir[a[1]] : dir[rand(rng, syma[a[1] .!= syma])]
            rdir2 = (rand(rng) < mdp.tprob) ? dir[a[2]] : dir[rand(rng, syma[a[2] .!= syma])]
            new_ego = ego_pos(s) + dir[mdp.ego_policy(s, rng)[1]]
            new_adv1 = adversary1_pos(s) + rdir1
            new_adv2 = adversary2_pos(s) + rdir2
        end

        # Make sure the moves are in bound and not hitting a wall
        new_ego = valid_pos(mdp, new_ego) ? new_ego : ego_pos(s)
        new_adv1 = valid_pos(mdp, new_adv1) ? new_adv1 : adversary1_pos(s)
        new_adv2 = valid_pos(mdp, new_adv2) ? new_adv2 : adversary2_pos(s)

        return (sp = S2(new_ego..., new_adv1..., new_adv2...), r = reward(mdp, s))
    end
end

# Returns the reward for the provided state
function POMDPs.reward(mdp::DoubleAdversarialGridworldMDP, s::S2)
    isterminal(mdp, s) && return 0
    r = (get(mdp.rewards, ego_pos(s), 0.0) - mdp.failure_penalty*agents_overlap(s)) /  mdp.failure_penalty
    mdp.agent_gets_action == :ego ? r : -r
end

function tocolor(mdp::DoubleAdversarialGridworldMDP, r::Float64)
    maxr = maximum(values(mdp.rewards))
    minr = -maxr
    frac = (r-minr)/(maxr-minr)
    return get(ColorSchemes.redgreensplit, frac)
end

# Renders the mdp
function POMDPModelTools.render(mdp::DoubleAdversarialGridworldMDP, s::S2)
    nx, ny = mdp.size
    cells = []
    for x in 1:nx, y in 1:ny
        pos = GWPos(x,y)
        reward_index = findfirst([pos] .== keys(mdp.rewards))
        wall_index = findfirst([pos] .== mdp.walls)
        ctx = context((x-1)/nx, (ny-y)/ny, 1/nx, 1/ny)
        color = "white"
        if !isnothing(reward_index)
            color = tocolor(mdp, get(mdp.rewards, pos, 0))
        elseif !isnothing(wall_index)
            color = "black"
        end
        cell = compose(ctx, Compose.rectangle(), fill(color))
        push!(cells, cell)
    end
    grid = compose(context(), Compose.stroke("gray"), cells...)
    outline = compose(context(), Compose.rectangle())

    ego = nothing
    x,y = ego_pos(s)
    if all([x,y] .> 0)
        agent_ctx = context((x-1)/nx, (ny-y)/ny, 1/nx, 1/ny)
        ego = compose(agent_ctx, Compose.circle(0.5, 0.5, 0.4), Compose.stroke("black"), fill("blue"))
    end

    adversary1 = nothing
    x,y = adversary1_pos(s)
    if all([x,y] .> 0)
        agent_ctx = context((x-1)/nx, (ny-y)/ny, 1/nx, 1/ny)
        adversary1 = compose(agent_ctx, Compose.circle(0.5, 0.5, 0.4), Compose.stroke("black"), fill("orange"))
    end

    adversary2 = nothing
    x,y = adversary2_pos(s)
    if all([x,y] .> 0)
        agent_ctx = context((x-1)/nx, (ny-y)/ny, 1/nx, 1/ny)
        adversary2 = compose(agent_ctx, Compose.circle(0.5, 0.5, 0.4), Compose.stroke("black"), fill("red"))
    end
    agents_comp = compose(context(), ego, adversary1, adversary2)

    sz = min(w, h)
    return compose(context((w-sz)/2, (h-sz)/2, sz, sz), agents_comp, grid, outline)
end


text:
@enum BranchEnum forward_branch=1 backward_branch=2 imaginary_branch=3

struct BranchPoint
  val::ComplexF64
  ref::Float64
  domain::BranchEnum
end

struct Branch
  domain::BranchEnum
  min_val::ComplexF64
  max_val::ComplexF64
  function Branch(d::BranchEnum, len::Real)
    @assert len > 0.0 "branch must have non-zero length (d = $d, len = $len)"
    if d == forward_branch
      return new(d, 0.0, len)
    elseif d == backward_branch
      return new(d, len, 0.0)
    else
      return new(d, 0.0, -1.0im * len)
    end
  end
end

"""
maps from the unit interval [0, 1] to a point on a branch
"""
function get_point(b::Branch, ref::Real)
  @assert 0.0 <= ref <= 1.0
  return BranchPoint(ref * b.max_val + (1 - ref) * b.min_val, ref, b.domain)
end

(b::Branch)(ref::Real) = get_point(b, ref)

function length(b::Branch)
  return abs(b.max_val - b.min_val)
end


text:
export companion

function companion(Z)

    T = eltype(Z)
    K = size(Z,3)
    @assert K > 1

    M, N = size(Z)[1:2]
    C = similar(Z, M*(K-1), N*(K-1))
    fill!(C,0)

    @assert M == N
    Id = Matrix{T}(I, M, N)
    for m in 2:K-1
        n = m-1
        C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = Id
    end

    W = -inv(Z[:,:,1])
    for n in 1:K-1
        m = 1
        C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = W*Z[:,:,n+1]
    end

    return C
end


text:
using DistributedArrays, Krylov, LinearOperators, MatrixMarket

function residuals(A, b, shifts, x)
  nshifts = size(shifts, 1);
  r = [ (b - A * x[:,i] - shifts[i] * x[:,i]) for i = 1 : nshifts ];
  return r;
end

# Parallel reduce.
preduce(func, darray) = reduce(func,
                               map(fetch,
                                   [ (@spawnat p reduce(func, localpart(darray))) for p = procs(darray) ] ));

"""The Lanczos version of the conjugate gradient method to solve a family
of shifted systems

  (A + αI) x = b  (α = α₁, α₂, ...)

The method does _not_ abort if A + αI is not definite.

In this version, the shifted systems are dispatched on the available processors,
and operations specific to each shift is carried out on the processor hosting it.
"""
function cg_lanczos_shift_par{Tb <: Real, Ts <: Real}(A :: AbstractLinearOperator, b :: Array{Tb,1}, shifts :: Array{Ts,1};
                                                      atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, itmax :: Int=0,
                                                      check_curvature :: Bool=false, verbose :: Bool=false)

  n = size(b, 1);
  (size(A, 1) == n & size(A, 2) == n) || error("Inconsistent problem size");

  dshifts = distribute(shifts);
  nshifts = size(shifts, 1);
  nchunks = length(dshifts.indexes)
  verbose && @printf("CG Lanczos: system of %d equations in %d variables with %d shifts\n", n, n, nshifts);

  # Initial state.
  ## Distribute x similarly to shifts.
  dx = dzeros((n, nshifts), workers(), [1, nchunks]);
  β = norm(b);
  β == 0 && return convert(Array, dx),
            Krylov.LanczosStats(true, zeros(nshifts), false, 0.0, 0.0, "x = 0
                                is a zero-residual solution")
  v = b / β;
  v_prev = v;

  # Distribute p similarly to shifts.
  # Initialize each p to b.
  dp = dzeros((n, nshifts), workers(), [1, nchunks]);
  @sync for proc in procs(dp)
    @spawnat proc begin
      p_loc = localpart(dp);
      for i = 1 : size(p_loc, 2)
        p_loc[:,i] = b;
      end
    end
  end

  # Keep track of shifted systems that have converged.
  dsolved = dfill(false, (nshifts,), workers(), [nchunks]);
  dconverged = dfill(false, (nshifts,), workers(), [nchunks]);
  dindefinite= dfill(false, (nshifts,), workers(), [nchunks]);
  iter = 0;
  itmax == 0 && (itmax = 2 * n);

  # Initialize some constants used in recursions below.
  dσ = dfill(β, (nshifts,), workers(), [nchunks]);
  dδ = dzeros((nshifts,), workers(), [nchunks]);
  dω = dzeros((nshifts,), workers(), [nchunks]);
  dγ = dones((nshifts,), workers(), [nchunks]);

  # Define stopping tolerance.
  drNorm = dfill(β, (nshifts,), workers(), [nchunks]);
  ε = atol + rtol * β;

  # Build format strings for printing.
  if verbose
    fmt = "%5d" * repeat("  %8.1e", nshifts) * "\n";
    c_printf(fmt, iter, drNorm...);
  end

  solved = false;
  tired = iter >= itmax;
  status = "unknown";

  # Main loop.
  while ! (solved || tired)
    # Form next Lanczos vector.
    v_next = A * v;
    δ = dot(v, v_next);
    BLAS.axpy!(n, -δ, v, 1, v_next, 1);  # Faster than v_next = v_next - δ * v;
    if iter > 0
      BLAS.axpy!(n, -β, v_prev, 1, v_next, 1);  # Faster than v_next = v_next - β * v_prev;
      v_prev = v;
    end
    β = norm(v_next);
    v = v_next / β;

    # Compute next CG iterate for each shift.
    @sync for proc in procs(dp)
      @spawnat proc begin
        solved_loc = localpart(dsolved);
        converged_loc = localpart(dconverged);
        indefinite_loc = localpart(dindefinite);
        shifts_loc = localpart(dshifts);
        δ_loc = localpart(dδ);
        γ_loc = localpart(dγ);
        ω_loc = localpart(dω);

        rNorm_loc = localpart(drNorm);

        # Check curvature: v'(A + sᵢI)v = v'Av + sᵢ ‖v‖² = δ + sᵢ because ‖v‖ = 1.
        # It is possible to show that σⱼ² (δⱼ + sᵢ - ωⱼ₋₁ / γⱼ₋₁) = pⱼᵀ (A + sᵢ I) pⱼ.
        # Stop iterating on indefinite problems if requested.
        δ_loc[:] = δ + shifts_loc
        γ_loc[:] = 1 ./ (δ_loc - ω_loc ./ γ_loc)
        indefinite_loc[:] |= (γ_loc .<= 0.0);
        not_cv = check_curvature ? find(! (converged_loc | indefinite_loc)) : find(! converged_loc);

        if length(not_cv) > 0
          # Fetch parts of relevant arrays for which the residual has not yet converged.
          σ_loc = localpart(dσ);
          x_loc = localpart(dx);
          p_loc = localpart(dp);
          shifts_loc = shifts_loc[not_cv];

          ω_loc[not_cv] = β * γ_loc[not_cv];
          σ_loc[not_cv] .*= -ω_loc[not_cv];
          ω_loc[not_cv] .*= ω_loc[not_cv];
          rNorm_loc[not_cv] = abs(σ_loc[not_cv]);

          # Can't seem to use BLAS3 GEMM calls here.
          x_loc[:, not_cv] += (p_loc[:, not_cv] * diagm(γ_loc[not_cv]));  # diagm?!
          p_loc[:, not_cv] = (v * σ_loc[not_cv]' + p_loc[:, not_cv] * diagm(ω_loc[not_cv]));

          # Update list of systems that have converged.
          converged_loc[not_cv] = rNorm_loc[not_cv] .<= ε;
        end

        solved_loc[:] = converged_loc | indefinite_loc;
      end
    end

    iter = iter + 1;
    verbose && c_printf(fmt, iter, drNorm...);

    solved = preduce(&, dsolved);
    tired = iter >= itmax;
  end

  status = tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol"
  stats = Krylov.LanczosStats(solved, convert(Array, drNorm), convert(Array, dindefinite), 0.0, 0.0, status);
  return (dx, stats);
end

cg_lanczos_shift_par{TA <: Number, Tb <: Number, Ts <: Number}(A :: Array{TA,2}, b :: Array{Tb,1}, shifts :: Array{Ts,1};
                                                               atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, itmax :: Int=0,
                                                               check_curvature :: Bool=false, verbose :: Bool=false) =
  cg_lanczos_shift_par(LinearOperator(A), b, shifts, atol=atol, rtol=rtol, itmax=itmax, check_curvature=check_curvature, verbose=verbose);

cg_lanczos_shift_par{TA <: Number, Tb <: Number, Ts <: Number, IA <: Integer}(A :: SparseMatrixCSC{TA,IA}, b :: Array{Tb,1}, shifts :: Array{Ts,1};
                                                                             atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, itmax :: Int=0,
                                                                             check_curvature :: Bool=false, verbose :: Bool=false) =
  cg_lanczos_shift_par(LinearOperator(A), b, shifts, atol=atol, rtol=rtol, itmax=itmax, check_curvature=check_curvature, verbose=verbose);

# mtx = "data/1138bus.mtx";
mtx = "data/bcsstk09.mtx";

A = MatrixMarket.mmread(mtx);
n = size(A, 1);
b = ones(n); b_norm = norm(b);

# Define a linear operator with preallocation.
Ap = zeros(n);
op = LinearOperator(n, n, true, true, p -> A_mul_B!(1.0,  A, p, 0.0, Ap))

# Solve (A+αI)x = b in parallel.
shifts = [1, 2, 3, 4];
(x, stats) = cg_lanczos_shift_par(op, b, shifts, verbose=false);
# @profile (x, stats) = cg_lanczos_shift_par(op, b, shifts);
@time (x, stats) = cg_lanczos_shift_par(op, b, shifts, verbose=false);
r = residuals(A, b, shifts, convert(Array, x));
resids = map(norm, r) / b_norm;
@printf("Relative residuals with shifts:\n");
for resid in resids
  @printf(" %8.1e", resid);
end
@printf("\n");


text:
module FiniteDiff

using LinearAlgebra, SparseArrays, StaticArrays, ArrayInterfaceCore, Requires

import Base: resize!

_vec(x) = vec(x)
_vec(x::Number) = x

_mat(x::AbstractMatrix) = x
_mat(x::StaticVector)   = reshape(x, (axes(x, 1),     SOneTo(1)))
_mat(x::AbstractVector) = reshape(x, (axes(x, 1), Base.OneTo(1)))

include("iteration_utils.jl")
include("epsilons.jl")
include("derivatives.jl")
include("gradients.jl")
include("jacobians.jl")
include("hessians.jl")

end # module


text:
using Documenter, MolecularIntegrals

makedocs(sitename="MolecularIntegrals.jl")
deploydocs(
    repo = "github.com/rpmuller/MolecularIntegrals.jl.git",
    devbranch="main"
)


text:
#nb # %% slideshow [markdown] {"slideshow": {"slide_type": "slide"}}
#md # # [Lecture 4 - Scope and Fruitful Functions](@id lecture4)
#!md # # Lecture 4 - Scope and Fruitful Functions
#
# **Essential Skills For Computational Biology**
#
# *2020-06-16*
#
# Kevin Bonham, PhD

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## "Scope" is where references live
#
# - A reference is a human-readable name for data or action
#   - eg. a variable, argument, or function name
# - All references have scope
#   - sometimes the scope is "global" (available everywhere)
#   - sometimes the scope is narrow (eg. inside a loop)
# - In different scopes, the same name may reference different things

#nb %% A slide [code] {"slideshow": {"slide_type": "subslide"}}
i = 3

for i in 1:5
    print(i, " ")
end

#-
#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
i

#-
#nb %% A slide [code] {"slideshow": {"slide_type": "subslide"}}
myvar = "Green eggs"
function ham(myvar)
    myvar = "ham"
    return myvar
end

#-
#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
myvar

#-
#nb %% A slide [code] {"slideshow": {"slide_type": "subslide"}}
function decrement4()
    j = 4
    while j > 0
        print(j, " ")
        j = j - 1
    end
    return j
end
decrement4()

#-
#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}

## copy to REPL
j = 4

while j > 0
    print(j, " ")
    j = j - 1
end

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Fruitful functions and the perils of printing
#
# - returned values are printed in the REPL
# - printed values are not (necessarily) returned

#nb %% A slide [code] {"slideshow": {"slide_type": "subslide"}}

function verbose(thing)
    println(thing)
    return thing
end
#-
#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
verbose(true)
#-
#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
verbose(1)
verbose(2)
verbose(3)
#-
#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
verbose("woah");

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## A reminder: I make mistakes
#
# - Several people have encountered tests that failed even
#   though they had correct answers
# - Don't bang your head against the wall (for too long)
# - We will still make you work for it
# - Not all difficulty is desirable

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## One other reminder: Helping eachother is ok

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# ### (Just don't copy-paste code from one another)

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Final project ideas
#
# - COVID data analysis
# - SARS-CoV2 evolution / strain analysis
# - Police violence data analysis
# - Your idea here


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
    x1'(t) = k_prod - k_deg * x1(t) - k1 * x1(t) * u(t),
    x2'(t) = -k3 * x2(t) - k_deg * x2(t) - a2 * x2(t) * x10(t) + t1 * x4(t) - a3 * x2(t) * x13(t) + t2 * x5(t) + (k1 * x1(t) - k2 * x2(t) * x8(t)) * u(t),
    x3'(t) = k3 * x2(t) - k_deg * x3(t) + k2 * x2(t) * x8(t) * u(t),
    x4'(t) = a2 * x2(t) * x10(t) - t1 * x4(t),
    x5'(t) = a3 * x2(t) * x13(t) - t2 * x5(t),
    x6'(t) = c_6a * x13(t) - a1 * x6(t) * x10(t) + t2 * x5(t) - i1 * x6(t),
    x7'(t) = i1 * kv * x6(t) - a1 * x11(t) * x7(t),
    x8'(t) = c4 * x9(t) - c5 * x8(t),
    x9'(t) = c2 + c1 * x7(t) - c3 * x9(t),
    x10'(t) = -a2 * x2(t) * x10(t) - a1 * x10(t) * x6(t) + c_4a * x12(t) - c_5a * x10(t) - i_1a * x10(t) + e_1a * x11(t),
    x11'(t) = -a1 * x11(t) * x7(t) + i_1a * kv * x10(t) - e_1a * kv * x11(t),
    x12'(t) = c_2a + c_1a * x7(t) - c_3a * x12(t),
    x13'(t) = a1 * x10(t) * x6(t) - c_6a * x13(t) - a3 * x2(t) * x13(t) + e_2a * x14(t),
    x14'(t) = a1 * x11(t) * x7(t) - e_2a * kv * x14(t),
    x15'(t) = c_2c + c_1c * x7(t) - c_3c * x15(t),
    y1(t) = x2(t),
    y2(t) = x10(t) + x13(t),
    y3(t) = x9(t),
    y4(t) = x1(t) + x2(t) + x3(t),
    y5(t) = x12(t),
    y6(t) = x7(t)
)

ode = SIAN.set_parameter_values(ode, OrderedDict(
    a1 => Nemo.QQ(1, 2),
    a2 => Nemo.QQ(1, 5),
    a3 => Nemo.QQ(1),
    c_1a => Nemo.QQ(5, 10^(7)),
    c_2a => Nemo.QQ(0),
    c_5a => Nemo.QQ(1, 10^(4)),
    c_6a => Nemo.QQ(2, 10^(5)),
    c1 => Nemo.QQ(5, 10^(7)),
    c2 => Nemo.QQ(0),
    c3 => Nemo.QQ(4, 10^(4)),
    c4 => Nemo.QQ(1, 2),
    kv => Nemo.QQ(5),
    e_1a => Nemo.QQ(5, 10^(4)),
    c_1c => Nemo.QQ(5, 10^(7)),
    c_2c => Nemo.QQ(0),
    c_3c => Nemo.QQ(4, 10^(4))
))

res = identifiability_ode(ode,  get_parameters(ode); p=0.99, p_mod=2^29 - 3, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia