text:
using SparseArrays, LinearAlgebra
struct Maze
    dimx::Int
    dimy::Int
    walls::Symmetric{Bool,SparseMatrixCSC{Bool,Int}}
end
Base.length(m::Maze) = m.dimx * m.dimy
Base.size(m::Maze) = (m.dimx, m.dimy)
function Maze(; nx = 40, ny = 40, nwalls = div(nx*ny, 20), rng = ENV_RNG) #Random.GLOBAL_RNG)
    m = emptymaze(nx, ny)
    for _ in 1:nwalls
        addrandomwall!(m, rng = rng)
    end
    breaksomewalls!(m, rng = rng)
end
function emptymaze(dimx, dimy)
    Maze(dimx, dimy, Symmetric(sparse([], [], Bool[], dimx * dimy, dimx * dimy)))
end
iswall(maze, pos, dir) = hitborder(maze, pos, dir) || maze.walls[pos, nextpos(maze, pos, dir)]
function setwall!(maze, pos, dir)
    npos = nextpos(maze, pos, dir)
    maze.walls.data[sort([pos, npos])...] = true
    maze
end
function breakwall!(maze, pos, dir)
    npos = nextpos(maze, pos, dir)
    maze.walls.data[sort([pos, npos])...] = false
    dropzeros!(maze.walls.data)
    maze
end
function nextpos(maze, pos, dir)
    dir == :up && return pos - 1
    dir == :down && return pos + 1
    dir == :left && return pos - maze.dimy
    dir == :right && return pos + maze.dimy
end
function hitborder(maze, pos, dir)
    dir == :up && pos % maze.dimy == 1 && return true
    dir == :down && pos % maze.dimy == 0 && return true
    dir == :left && pos <= maze.dimy && return true
    dir == :right && pos > maze.dimy * (maze.dimx - 1) && return true
    return false
end

function orthogonal_directions(dir)
    dir ∈ (:up, :down) && return (:left, :right)
    return (:up, :down)
end

function is_wall_neighbour(maze, pos)
    for dir in (:up, :down, :left, :right)
        iswall(maze, pos, dir) && return true
    end
    for npos in (nextpos(maze, pos, :up), nextpos(maze, pos, :down))
        for dir in (:left, :right, :up, :down)
            iswall(maze, npos, dir) && return true
        end
    end
    return false
end
function is_wall_tangential(maze, pos, dir)
    for ortho_dir in orthogonal_directions(dir)
        iswall(maze, pos, ortho_dir) && return true
    end
    return false
end
is_wall_ahead(maze, pos, dir) = iswall(maze, pos, dir)

function addrandomwall!(maze; rng = ENV_RNG) # Random.GLOBAL_RNG)
    potential_startpos = filter(x -> !is_wall_neighbour(maze, x), 1:maze.dimx * maze.dimy)
    if potential_startpos == []
        @warn("Cannot add a random wall.")
        return maze
    end
    pos = rand(rng, potential_startpos)
    direction = rand(rng, (:up, :down, :left, :right))
    while true
        setwall!(maze, pos, orthogonal_directions(direction)[2])
        pos = nextpos(maze, pos, direction)
        is_wall_tangential(maze, pos, direction) && break
        if is_wall_ahead(maze, pos, direction)
            setwall!(maze, pos, orthogonal_directions(direction)[2])
            break
        end
    end
    return maze
end

function n_effective(n, f, list)
    N = n === nothing ? div(length(list), Int(1/f)) : n
    min(N, length(list))
end
function breaksomewalls!(m; f = 1/50, n = nothing, rng = ENV_RNG) #Random.GLOBAL_RNG)
    wall_idxs = findall(m.walls.data)
    pos = sample(rng, wall_idxs, n_effective(n, f, wall_idxs), replace = false)
    m.walls.data[pos] .= false
    dropzeros!(m.walls.data)
    m
end
# function addobstacles!(m; f = 1/100, n = nothing, rng = Random.GLOBAL_RNG)
#     nz = findall(x -> x == 1, reshape(m, :))
#     pos = sample(rng, nz, n_effective(n, f, nz), replace = false)
#     m[pos] .= 0
#     m
# end
setTandR!(d) = for s in 1:length(d.maze) setTandR!(d, s) end
function setTandR!(d, s)
    # @show s
    T = d.mdp.trans_probs
    R = d.mdp.reward
    goals = d.goals
    # @show goals
    ns = d.mdp.observationspace.n
    maze = d.maze
    idx_goals = findfirst(x -> x == s, goals)
    # @show idx_goals
    # @show d.goalrewards
    #### TODO #####
    # Here R should be reset to 0s everwhere!!!!
    # Otherwise what you had before stay same!!!!!
    if idx_goals !== nothing
        # @show d.goalrewards[idx_goals]
        R.value[s] = d.goalrewards[idx_goals]
    end
    # @show unique(R.value)
    if (!(in(:chosenactionweight, fieldnames(typeof(d))))
        || (d.chosenactionweight == 1.))
        for (aind, a) in enumerate((:up, :down, :left, :right))
            # @show (aind, a)
            npos = iswall(maze, s, a) ? s : nextpos(maze, s, a)
            # @show npos
            if d.neighbourstateweight > 0
                positions = [npos]
                # @show positions
                weights = [1.]
                for dir in (:up, :down, :left, :right)
                    if !iswall(maze, npos, dir)
                        push!(positions, nextpos(maze, npos, dir))
                        push!(weights, d.neighbourstateweight)
                    end
                end
                # @show positions
                # @show weights
                weights /= sum(weights)
                # @show weights
                T[aind, s] = sparsevec(positions, weights, ns)
            else
                T[aind, s] = sparsevec([npos], [1.], ns)
            end
        end
    else
        positions = []
        for (aind, a) in enumerate((:up, :down, :left, :right))
            # @show (aind, a)
            npos = iswall(maze, s, a) ? s : nextpos(maze, s, a)
            push!(positions, npos)
        end
        # @show positions
        for (aind, a) in enumerate((:up, :down, :left, :right))
            # @show (aind, a)
            weights = (1. - d.chosenactionweight)/ 3. * ones(length(positions))
            weights[aind] = d.chosenactionweight
            weights /= sum(weights)
            # @show weights
            T[aind, s] = sparsevec(positions, weights, ns)
            # @show T[aind, s]
        end
    end
end

"""
    struct DiscreteMaze
        mdp::MDP
        maze::Array{Int, 2}
        goals::Array{Int, 1}
        statefrommaze::Array{Int, 1}
        mazefromstate::Array{Int, 1}
"""
struct DiscreteMaze{T}
    mdp::T
    maze::Maze
    goals::Array{Int, 1}
    goalrewards::Array{Float64, 1}
    neighbourstateweight::Float64
    chosenactionweight::Float64
end
"""
    DiscreteMaze(; nx = 40, ny = 40, nwalls = div(nx*ny, 10), ngoals = 1,
                   goalrewards = 1, stepcost = 0, stochastic = false,
                   neighbourstateweight = .05)
Returns a `DiscreteMaze` of width `nx` and height `ny` with `nwalls` walls and
`ngoals` goal locations with reward `goalreward` (a list of different rewards
for the different goal states or constant reward for all goals), cost of moving
`stepcost` (reward = -`stepcost`); if `stochastic = true` the actions lead with
a certain probability to a neighbouring state, where `neighbourstateweight`
controls this probability.
"""
DiscreteMaze(maze; kwargs...) = DiscreteMaze(; maze = maze, kwargs...)
function DiscreteMaze(; nx = 40, ny = 40, nwalls = div(nx*ny, 20),
                      rng = ENV_RNG, #Random.GLOBAL_RNG,
                      maze = Maze(nx = nx, ny = ny, nwalls = nwalls, rng = rng),
                      ngoals = 1,
                      goalrewards = 1.,
                      stepcost = 0,
                      stochastic = false,
                      neighbourstateweight = stochastic ? .05 : 0.,
                      chosenactionweight = 1.)
    na = 4
    ns = length(maze)
    legalstates = 1:ns
    T = Array{SparseVector{Float64,Int}}(undef, na, ns)
    goals = sort(sample(rng, legalstates, ngoals, replace = false))
    R = DeterministicNextStateReward(fill(-stepcost, ns))
    isterminal = zeros(Int, ns); isterminal[goals] .= 1
    isinitial = setdiff(legalstates, goals)
    res = DiscreteMaze(MDP(DiscreteSpace(ns, 1),
                           DiscreteSpace(na, 1),
                           rand(rng, legalstates),
                           T, R,
                           isinitial,
                           isterminal),
                       maze,
                       goals,
                       typeof(goalrewards) <: Number ? fill(goalrewards, ngoals) :
                                                       goalrewards,
                       neighbourstateweight,
                       chosenactionweight)
    setTandR!(res)
    res
end
interact!(env::DiscreteMaze, a) = interact!(env.mdp, a)
reset!(env::DiscreteMaze) = reset!(env.mdp)
getstate(env::DiscreteMaze) = getstate(env.mdp)
actionspace(env::DiscreteMaze) = actionspace(env.mdp)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
mutable struct ChangeDiscreteMaze{DiscreteMaze}
    discretemaze::DiscreteMaze
    stepcounter::Int
    switchsteps::Array{Int, 1}
    switchpos::Array{Int, 1}
    switchdir::Array{Symbol, 1}
    switchflag::Array{Bool, 2} # Used for RecordSwitches callback
    chosenactionweight::Float64
end
function ChangeDiscreteMaze(; nx = 40, ny = 40, nwalls = div(nx*ny, 20),
                            rng = ENV_RNG, #Random.GLOBAL_RNG,
                            maze = Maze(nx = nx, ny = ny, nwalls = nwalls, rng = rng),
                            stochastic = false,
                            ngoals = 1,  goalrewards = 1.,
                            nswitches = 1,
                            switchpos = rand(ENV_RNG, 1:length(maze), nswitches),
                            switchdir = rand(ENV_RNG, [:up, :down, :left, :right], nswitches),
                            switchsteps = 10^2 * ones(Int, length(switchpos)),
                            neighbourstateweight = stochastic ? .05 : 0.,
                            chosenactionweight = 1.)
    dm = DiscreteMaze(maze = maze, ngoals = ngoals, goalrewards=goalrewards,
                        neighbourstateweight = neighbourstateweight,
                        chosenactionweight = chosenactionweight)
    switchflag = fill(false, 4, length(dm.maze))
    ChangeDiscreteMaze(dm, 0, switchsteps, switchpos, switchdir, switchflag, chosenactionweight)
end
ChangeDiscreteMaze(maze; kwargs...) = ChangeDiscreteMaze(; maze = maze, kwargs...)
ChangeDiscreteMaze(maze, switchpos, switchdir; kwargs...) =
    ChangeDiscreteMaze(; maze = maze, switchpos = switchpos, switchdir = switchdir, kwargs...)
# --- For standard RL learners all we need is:
function interact!(env::ChangeDiscreteMaze, action)
    env.stepcounter += 1
    env.switchflag .= false
    #@show env.stepcounter
    if any(env.stepcounter .== env.switchsteps)
        # println(" ------- Switch now!")
        # @show env.stepcounter
        # && (env.discretemaze.mdp.isterminal[env.discretemaze.mdp.state] == 1))# Switch or not!
        #println("####################################################################################################")
        nswitches = findall(env.stepcounter .== env.switchsteps)
        previousT = deepcopy(env.discretemaze.mdp.trans_probs)
        for i in nswitches
            setupswitch!(env, i)
        end
        updateswitchflag!(env, previousT)
    end
    interact!(env.discretemaze.mdp, action)
end

function setupswitch!(env::ChangeDiscreteMaze, iswitch)
    if iswall(env.discretemaze.maze, env.switchpos[iswitch], env.switchdir[iswitch])
        breakwall!(env.discretemaze.maze, env.switchpos[iswitch], env.switchdir[iswitch])
    else
        setwall!(env.discretemaze.maze, env.switchpos[iswitch], env.switchdir[iswitch])
    end
    # env.discretemaze.maze[env.switchpos[iswitch]] = 1 - env.discretemaze.maze[env.switchpos[iswitch]]
    setTandR!(env.discretemaze)
end
# # --- For MDP learner we need this:
# function interact!(env::ChangeDiscreteMaze, action)
#     env.stepcounter += 1
#     env.switchflag .= false
#     #@show env.stepcounter
#     if any(env.stepcounter .== env.switchsteps)
#         previousT = deepcopy(env.discretemaze.mdp.trans_probs)
#         # && (env.discretemaze.mdp.isterminal[env.discretemaze.mdp.state] == 1))# Switch or not!
#         #println("###########################################################")
#         nswitches = findall(env.stepcounter .== env.switchsteps)
#         # @show env.switchsteps
#         # @show env.switchpos
#         if any(env.switchpos .== env.discretemaze.mdp.state)
#             for i in nswitches
#                 push!(env.switchsteps, env.stepcounter+1)
#                 push!(env.switchpos, env.switchpos[i])
#             end
#             # @show env.switchsteps
#             # @show env.switchpos
#             # println("i saw it!")
#         else
#             for i in nswitches
#                 # @show i
#                 setupswitch!(env, i)
#             end
#         end
#         updateswitchflag!(env, previousT)
#     end
#     interact!(env.discretemaze.mdp, action)
# end
# function setupswitch!(env::ChangeDiscreteMaze, iswitch)
#     #for i in 1:length(env.switchpos)
#     env.discretemaze.maze[env.switchpos[iswitch]] = 1 - env.discretemaze.maze[env.switchpos[iswitch]]
#     # @show env.discretemaze.maze
#     #end
#     # # !!!!! Reset the whole transition matrix, so that new walls become "undef"!!!
#     # NOTE: Not important for normal RL learners. But important for MDPlearner!
#     ns = env.discretemaze.mdp.observationspace.n
#     na = env.discretemaze.mdp.actionspace.n
#     env.discretemaze.mdp.trans_probs[:] = Array{SparseVector{Float64,Int}}(undef, na, ns)
#
#     setTandR!(env.discretemaze)
# end
function updateswitchflag!(env, previousT)
    for i in 1:length(env.switchflag)
        # @show i
        if !isassigned(previousT, i)
            if isassigned(env.discretemaze.mdp.trans_probs, i)
                env.switchflag[i] = true
                # @show i
                # @show env.discretemaze.mdp.trans_probs[i]
                # println("Case 1")
            end
        else
            if !isassigned(env.discretemaze.mdp.trans_probs, i)
                env.switchflag[i] = true
                # println("Case 2")
                # @show i
                # @show previousT[i]
            else
                if !all(previousT[i] .== env.discretemaze.mdp.trans_probs[i])
                    env.switchflag[i] = true
                    # println("Case 3")
                    # @show previousT[i]
                    # @show env.discretemaze.mdp.trans_probs[i]
                end
            end
        end
    end
end
reset!(env::ChangeDiscreteMaze) = reset!(env.discretemaze.mdp)
getstate(env::ChangeDiscreteMaze) = getstate(env.discretemaze.mdp)
actionspace(env::ChangeDiscreteMaze) = actionspace(env.discretemaze.mdp)
plotenv(env::ChangeDiscreteMaze) = plotenv(env.discretemaze)

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

mutable struct ChangeDiscreteMazeProbabilistic{DiscreteMaze}
    discretemaze::DiscreteMaze
    changeprobability::Float64
    switchpos::Array{Int, 1}
    switchdir::Array{Symbol, 1}
    switchflag::Array{Bool, 2} # Used for RecordSwitches callback
    chosenactionweight::Float64
    seed::Any
    rng::MersenneTwister # Used only for switches!
end
function ChangeDiscreteMazeProbabilistic(; nx = 40, ny = 40, nwalls = div(nx*ny, 20),
                            seed = 3,
                            rng = MersenneTwister(seed),
                            maze = Maze(nx = nx, ny = ny, nwalls = nwalls, rng = rng),
                            ngoals = 1,  goalrewards = 1.,
                            changeprobability = .01,
                            nswitches = 1,
                            switchpos = rand(ENV_RNG, 1:length(maze), nswitches),
                            switchdir = rand(ENV_RNG, [:up, :down, :left, :right], nswitches),
                            stochastic = false,
                            neighbourstateweight = stochastic ? .05 : 0.,
                            chosenactionweight = 1.)
    dm = DiscreteMaze(maze = maze, ngoals = ngoals, goalrewards=goalrewards,
                        neighbourstateweight = neighbourstateweight,
                        chosenactionweight = chosenactionweight)
    switchflag = fill(false, 4, length(dm.maze))
    ChangeDiscreteMazeProbabilistic(dm, changeprobability, switchpos, switchdir,
                            switchflag, chosenactionweight, seed, rng)
end
ChangeDiscreteMazeProbabilistic(maze; kwargs...) =
    ChangeDiscreteMazeProbabilistic(; maze = maze, kwargs...)
ChangeDiscreteMazeProbabilistic(maze, switchpos, switchdir; kwargs...) =
    ChangeDiscreteMazeProbabilistic(; maze = maze, switchpos = switchpos, switchdir = switchdir, kwargs...)

# --- For standard RL learners all we need is:
function interact!(env::ChangeDiscreteMazeProbabilistic, action)
    env.switchflag .= false
    r = rand(env.rng)

    if r < env.changeprobability # Switch or not!
        println("Switch!")
        @show env.discretemaze.mdp.state
        # @show env.switchflag[:, 76]
        # @show env.discretemaze.mdp.trans_probs[1, 76]
        nswitches = length(env.switchpos)
        previousT = deepcopy(env.discretemaze.mdp.trans_probs)
        for i in nswitches
            setupswitch!(env, i)
        end
        updateswitchflag!(env, previousT)
        # @show env.switchflag[:, 76]
        # @show env.discretemaze.mdp.trans_probs[1, 76]
        nswitches = length(env.switchpos)
    end
    interact!(env.discretemaze.mdp, action)
end
# # --- For MDP learner we need this:
# function interact!(env::ChangeDiscreteMaze, action)
#     env.stepcounter += 1
#     env.switchflag .= false
#     #@show env.stepcounter
#     if any(env.stepcounter .== env.switchsteps)
#         previousT = deepcopy(env.discretemaze.mdp.trans_probs)
#         # && (env.discretemaze.mdp.isterminal[env.discretemaze.mdp.state] == 1))# Switch or not!
#         #println("###########################################################")
#         nswitches = findall(env.stepcounter .== env.switchsteps)
#         # @show env.switchsteps
#         # @show env.switchpos
#         if any(env.switchpos .== env.discretemaze.mdp.state)
#             for i in nswitches
#                 push!(env.switchsteps, env.stepcounter+1)
#                 push!(env.switchpos, env.switchpos[i])
#             end
#             # @show env.switchsteps
#             # @show env.switchpos
#             # println("i saw it!")
#         else
#             for i in nswitches
#                 # @show i
#                 setupswitch!(env, i)
#             end
#         end
#         updateswitchflag!(env, previousT)
#     end
#     interact!(env.discretemaze.mdp, action)
# end
# function setupswitch!(env::ChangeDiscreteMaze, iswitch)
#     #for i in 1:length(env.switchpos)
#     env.discretemaze.maze[env.switchpos[iswitch]] = 1 - env.discretemaze.maze[env.switchpos[iswitch]]
#     # @show env.discretemaze.maze
#     #end
#     # # !!!!! Reset the whole transition matrix, so that new walls become "undef"!!!
#     # NOTE: Not important for normal RL learners. But important for MDPlearner!
#     ns = env.discretemaze.mdp.observationspace.n
#     na = env.discretemaze.mdp.actionspace.n
#     env.discretemaze.mdp.trans_probs[:] = Array{SparseVector{Float64,Int}}(undef, na, ns)
#
#     setTandR!(env.discretemaze)
# end
function setupswitch!(env::ChangeDiscreteMazeProbabilistic, iswitch)
    if iswall(env.discretemaze.maze, env.switchpos[iswitch], env.switchdir[iswitch])
        breakwall!(env.discretemaze.maze, env.switchpos[iswitch], env.switchdir[iswitch])
    else
        setwall!(env.discretemaze.maze, env.switchpos[iswitch], env.switchdir[iswitch])
    end
    # env.discretemaze.maze[env.switchpos[iswitch]] = 1 - env.discretemaze.maze[env.switchpos[iswitch]]
    setTandR!(env.discretemaze)
end
reset!(env::ChangeDiscreteMazeProbabilistic) = reset!(env.discretemaze.mdp)
getstate(env::ChangeDiscreteMazeProbabilistic) = getstate(env.discretemaze.mdp)
actionspace(env::ChangeDiscreteMazeProbabilistic) = actionspace(env.discretemaze.mdp)
plotenv(env::ChangeDiscreteMazeProbabilistic) = plotenv(env.discretemaze)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# NOTE: RandomChangeDiscreteMaze is not tested
# ------------------------------------------------------------------------------
mutable struct RandomChangeDiscreteMaze{DiscreteMaze}
    discretemaze::DiscreteMaze
    nbreak::Int
    nadd::Int
    changeprobability::Float64
    switchflag::Array{Bool, 2} # Used for RecordSwitches callback.
    seed::Any
    rng::MersenneTwister # Used only for switches!
    chosenactionweight::Float64
end
function RandomChangeDiscreteMaze(; nx = 20, ny = 20, nwalls = 10,
                                    seed = 3, rng = MersenneTwister(seed),
                            maze = Maze(nx = nx, ny = ny, nwalls = nwalls, rng = rng),
                            ngoals = 4, goalrewards = 1.,
                    stochastic = false,
                    neighbourstateweight = stochastic ? .05 : 0.,
                    chosenactionweight = 2. /3.,
                    nbreak = 2, nadd = 4, changeprobability = 0.01)
    dm = DiscreteMaze(maze = maze, ngoals = ngoals, goalrewards=goalrewards,
                        stochastic = stochastic,
                        neighbourstateweight = neighbourstateweight,
                        chosenactionweight = chosenactionweight)
    switchflag = fill(false, 4, length(dm.maze))
    RandomChangeDiscreteMaze(dm, nbreak, nadd, changeprobability, switchflag,
                            seed, rng, chosenactionweight)
end
RandomChangeDiscreteMaze(maze; kwargs...) = RandomChangeDiscreteMaze(; maze = maze, kwargs...)
function interact!(env::RandomChangeDiscreteMaze, action)
    env.switchflag .= false
    r = rand(env.rng)
    if r < env.changeprobability # Switch or not!
        previousT = deepcopy(env.discretemaze.mdp.trans_probs)
        setupswitch!(env)
        updateswitchflag!(env, previousT)
    end
    interact!(env.discretemaze.mdp, action)
end
function setupswitch!(env::RandomChangeDiscreteMaze)
    breaksomewalls!(env.discretemaze.maze, n = env.nbreak, rng = env.rng)
    addobstacles!(env.discretemaze.maze, n = env.nadd, rng = env.rng)
    setTandR!(env.discretemaze)
end
reset!(env::RandomChangeDiscreteMaze) = reset!(env.discretemaze.mdp)
getstate(env::RandomChangeDiscreteMaze) = getstate(env.discretemaze.mdp)
actionspace(env::RandomChangeDiscreteMaze) = actionspace(env.discretemaze.mdp)
plotenv(env::RandomChangeDiscreteMaze) = plotenv(env.discretemaze)
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
function plotenv(env::DiscreteMaze)
    goals = env.goals
    # px = 6
    px = 16
    m = zeros(px*env.maze.dimx, px*env.maze.dimy)
    for i in 1:env.maze.dimx
        for j in 1:env.maze.dimy
            s = (i - 1) * env.maze.dimy + j
            iswall(env.maze, s, :down) && (m[(i-1)*px+1:i*px, j*px] .= -3)
            iswall(env.maze, s, :up) && (m[(i-1)*px+1:i*px, (j-1)*px+1] .= -3)
            iswall(env.maze, s, :right) && (m[i*px, (j-1)*px+1:j*px] .= -3)
            iswall(env.maze, s, :left) && (m[(i-1)*px+1, (j-1)*px+1:j*px] .= -3)
            if s == env.mdp.state
                val = -1
            elseif s in goals
                val = 1
            else
                continue
            end
            m[(i-1)*px+2:i*px-1, (j-1)*px+2:j*px-1] .= val
        end
    end
    imshow(m, clim = (-3, 3), colormap = 42)
end


text:
@test "this file is loaded" begin
  @assert 1 == 1
end

ext: jl


text:
module Peaks

include("minmax.jl")
include("peakprom.jl")

end # module Peaks


text:
using Test
import IterTools

using ChpSim

@testset "ChpSim.jl" begin

@testset "convert_result" begin
    @test convert(Bool, ChpSim.MeasureResult(true, determined=false)) == true
    @test convert(Complex, ChpSim.MeasureResult(true, determined=false)) ==
        Complex(true, false)
end

@testset "identity" begin
    s = ChpState(1)
    @test measure!(s, 1) == ChpSim.MeasureResult(false, determined=true)
end

@testset "bit_packed_table" begin
    s = ChpState(2)
    hadamard!(s, 1)
    phase!(s, 1)
    phase!(s, 1)
    hadamard!(s, 1)
    cnot!(s, 1, 2)
    @test measure!(s, 1) == ChpSim.MeasureResult(true, determined=true)
    @test measure!(s, 2) == ChpSim.MeasureResult(true, determined=true)
end

@testset "bit_flip" begin
    s = ChpState(1)
    hadamard!(s, 1)
    phase!(s, 1)
    phase!(s, 1)
    hadamard!(s, 1)
    @test measure!(s, 1) == ChpSim.MeasureResult(true, determined=true)
end

@testset "identity_2" begin
    s = ChpState(2)
    @test measure!(s, 1) == ChpSim.MeasureResult(false, determined=true)
    @test measure!(s, 2) == ChpSim.MeasureResult(false, determined=true)
end

@testset "bit_flip_2" begin
    s = ChpState(2)
    hadamard!(s, 1)
    phase!(s, 1)
    phase!(s, 1)
    hadamard!(s, 1)
    @test measure!(s, 1) == ChpSim.MeasureResult(true, determined=true)
    @test measure!(s, 2) == ChpSim.MeasureResult(false, determined=true)
end

@testset "epr" begin
    s = ChpState(2)
    hadamard!(s, 1)
    cnot!(s, 1, 2)
    v1 = measure!(s, 1)
    v2 = measure!(s, 2)
    @test !v1.determined
    @test v2.determined
    @test v1.value == v2.value
end

@testset "_pauli_product_phase" begin
    paulis = Bool[0 0  # I
                  1 0  # X
                  1 1  # Y
                  0 1] # Z
    expected = [0  0  0  0
                0  0  1 -1
                0 -1  0  1
                0  1 -1  0]
    @test ChpSim._pauli_product_phase(true, false, true, true) == 1
    for i in 1:4, j in 1:4
        @test (ChpSim._pauli_product_phase(paulis[i, :]..., paulis[j, :]...)
               == expected[i, j])
    end
end

@testset "phase_kickback_consume_s_state" begin
    s = ChpState(2)
    hadamard!(s, 2)
    phase!(s, 2)
    hadamard!(s, 1)
    cnot!(s, 1, 2)
    v1 = measure!(s, 2)
    @test !v1.determined
    if v1.value
        phase!(s, 1)
        phase!(s, 1)
    end
    phase!(s, 1)
    hadamard!(s, 1)
    @test measure!(s, 1) == ChpSim.MeasureResult(true, determined=true)
end

@testset "phase_kickback_preserve_s_state" begin
    s = ChpState(2)

    hadamard!(s, 2)
    phase!(s, 2)

    hadamard!(s, 1)

    cnot!(s, 1, 2)
    hadamard!(s, 2)
    cnot!(s, 1, 2)
    hadamard!(s, 2)

    phase!(s, 1)
    hadamard!(s, 1)
    @test measure!(s, 1) == ChpSim.MeasureResult(true, determined=true)
    phase!(s, 2)
    hadamard!(s, 2)
    @test measure!(s, 2) == ChpSim.MeasureResult(true, determined=true)
end

@testset "kickback_vs_stabilizer" begin
    s = ChpState(3)
    hadamard!(s, 3)
    cnot!(s, 3, 1)
    cnot!(s, 3, 2)
    phase!(s, 1)
    phase!(s, 2)
    hadamard!(s, 1)
    hadamard!(s, 2)
    hadamard!(s, 3)
    @test repr(MIME("text/plain"), s) == """
        -Y..
        -.Y.
        +..X
        ----
        +X.X
        +.XX
        +YYZ"""
    v1 = measure!(s, 1, bias=0)
    @test repr(MIME("text/plain"), s) == """
        +X.X
        -.Y.
        +..X
        ----
        +Z..
        +.XX
        +ZYY"""
    v2 = measure!(s, 2, bias=0)
    @test repr(MIME("text/plain"), s) == """
        +X.X
        +.XX
        +..X
        ----
        +Z..
        +.Z.
        -ZZZ"""
    v3 = measure!(s, 3, bias=0)
    @test repr(MIME("text/plain"), s) == """
        +X.X
        +.XX
        +..X
        ----
        +Z..
        +.Z.
        -ZZZ"""
    @test v1 == ChpSim.MeasureResult(false, determined=false)
    @test v2 == ChpSim.MeasureResult(false, determined=false)
    @test v3 == ChpSim.MeasureResult(true, determined=true)
end

@testset "s_state_distillation_low_depth" begin
    for _ in 1:100
        s = ChpState(9)

        stabilizers = [
            (1, 2, 3, 4),
            (1, 2, 5, 6),
            (1, 3, 5, 7),
            (2, 3, 5, 8),
        ]
        checks = [
            (s=[1], q=stabilizers[1]),
            (s=[2], q=stabilizers[2]),
            (s=[3], q=stabilizers[3]),
        ]

        stabilizer_measurements = Bool[]
        anc = 9
        for stab in stabilizers
            hadamard!(s, anc)
            for k in stab
                cnot!(s, anc, k)
            end
            hadamard!(s, anc)
            v = measure!(s, anc)
            @test !v.determined
            if v.value
                hadamard!(s, anc)
                phase!(s, anc)
                phase!(s, anc)
                hadamard!(s, anc)
            end
            push!(stabilizer_measurements, v.value)
        end

        qubit_measurements = Bool[]
        for k in 1:7
            phase!(s, k)
            hadamard!(s, k)
            push!(qubit_measurements, measure!(s, k).value)
        end

        if (sum(stabilizer_measurements) + sum(qubit_measurements)) & 1 == 1
            phase!(s, 8)
            phase!(s, 8)
        end

        phase!(s, 8)
        hadamard!(s, 8)
        r = measure!(s, 8)

        @test r == ChpSim.MeasureResult(false, determined=true)
        for c in checks
            rvs = [stabilizer_measurements[k] for k in c.s]
            rms = [qubit_measurements[k] for k in c.q]
            @test (sum(rvs) + sum(rms)) & 1 == 0
        end
    end
end

@testset "s_state_distillation_low_space" begin
    for _ in 1:100
        s = ChpState(5)

        phasors = [
            (1,),
            (2,),
            (3,),
            (1, 2, 3),
            (1, 2, 4),
            (1, 3, 4),
            (2, 3, 4),
        ]

        anc = 5
        for phasor in phasors
            hadamard!(s, anc)
            for k in phasor
                cnot!(s, anc, k)
            end
            hadamard!(s, anc)
            phase!(s, anc)
            hadamard!(s, anc)
            v = measure!(s, anc)
            @test !v.determined
            if v.value
                for k in (phasor..., anc)
                    hadamard!(s, k)
                    phase!(s, k)
                    phase!(s, k)
                    hadamard!(s, k)
                end
            end
        end

        for k in 1:3
            @test measure!(s, k) == ChpSim.MeasureResult(false, determined=true)
        end
        phase!(s, 4)
        hadamard!(s, 4)
        @test measure!(s, 4) == ChpSim.MeasureResult(true, determined=true)
    end
end

@testset "count_s_state_distillation_failure_cases" begin
    distill(errors)::String = distill(Set(errors))
    function distill(errors::Set{Int})::String
        s = ChpState(5)

        phasors = [
            (1,),
            (2,),
            (3,),
            (1, 2, 3),
            (1, 2, 4),
            (1, 3, 4),
            (2, 3, 4),
        ]

        anc = 5
        for (e, phasor) in enumerate(phasors)
            hadamard!(s, anc)
            for k in phasor
                cnot!(s, anc, k)
            end
            hadamard!(s, anc)
            phase!(s, anc)

            if e in errors
                phase!(s, anc)
                phase!(s, anc)
            end

            hadamard!(s, anc)
            v = measure!(s, anc)
            @test !v.determined
            if v.value
                hadamard!(s, anc)
                phase!(s, anc)
                phase!(s, anc)
                hadamard!(s, anc)
                for k in phasor
                    hadamard!(s, k)
                    phase!(s, k)
                    phase!(s, k)
                    hadamard!(s, k)
                end
            end
        end

        phase!(s, 4)
        phase!(s, 4)
        phase!(s, 4)
        hadamard!(s, 4)
        v = measure!(s, 4)
        @test v.determined
        hadamard!(s, 4)
        phase!(s, 4)
        checks = [measure!(s, k) for k in 1:3]
        @test all(e.determined for e in checks)
        good_result = !v.value
        checks_passed = !any(e.value for e in checks)
        return if checks_passed
            good_result ? "good" : "ERROR"
        else
            good_result ? "victim" : "caught"
        end
    end

    function classify(errs)::Dict{String, Int}
        result::Dict{String, Int} = Dict()
        for err in errs
            r = distill(err)
            result[r] = get(result, r, 0) + 1
        end
        result
    end

    nones = IterTools.subsets(1:7, 0)
    singles = IterTools.subsets(1:7, 1)
    doubles = IterTools.subsets(1:7, 2)
    triples = IterTools.subsets(1:7, 3)

    @test classify(nones) == Dict("good"=>1)
    @test classify(singles) == Dict("caught"=>3, "victim"=>4)
    @test classify(doubles) == Dict("caught"=>12, "victim"=>9)
    @test classify(triples) == Dict("caught"=>12, "victim"=>16, "ERROR"=>7)
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 64500da0-15b3-11ec-2965-5dbae2ce1ca3
using GMT
mozambique = gadm("MOZ");

# ╔═╡ c235d97e-1587-11ec-0438-81130f8f5239
md"### Plot countries administrative units from [GADM](https://gadm.org)

We start by download and save in cache the geopackage file with Mozambique data.
Note, a downloading message will be printed only once, and in this case it will say

*Downloading geographic data for country MOZ provided by the https://gadm.org project. It may take a while.
The file gadm36_MOZ.gpkg (after uncompressing) will be stored in c:/j/.gmt\cache*
"

# ╔═╡ f1be5140-1588-11ec-0705-918d55480e3a
imshow(mozambique, proj=:guess, title="Mozambique")

# ╔═╡ 12178d42-15b5-11ec-11a3-01a7cb72aff1
md"
Next let us add all the provincies of Mozambique. For that we'll use the country code as parent and the option `children=true` to indicate that we want all provincies boundaries."

# ╔═╡ 89b7aa40-15b2-11ec-01d0-37cce7e20f24
mozambique = gadm("MOZ", children=true);
imshow(mozambique, proj=:guess, title="Provinces of Mozambique")

# ╔═╡ 734b2090-15b5-11ec-3394-d99ce6470590
md"To know the provinces names such that we can use them individually for example, we use the option `names=true`"

# ╔═╡ 54c19ee0-15b2-11ec-3cf0-b7e5cdc75d3a
gadm("MOZ", names=true)

# ╔═╡ a658e3f0-15b5-11ec-0809-69ce7675d096
md"
Now we can plot only one of those provinces and its children"

# ╔═╡ 75208430-15b2-11ec-0ae1-3724dade5f8e
CD = gadm("MOZ", "Cabo Delgado", children=true);
imshow(CD, proj=:guess, title="Cabo Delgado")

# ╔═╡ Cell order:
# ╠═c235d97e-1587-11ec-0438-81130f8f5239
# ╠═64500da0-15b3-11ec-2965-5dbae2ce1ca3
# ╠═f1be5140-1588-11ec-0705-918d55480e3a
# ╠═12178d42-15b5-11ec-11a3-01a7cb72aff1
# ╠═89b7aa40-15b2-11ec-01d0-37cce7e20f24
# ╠═734b2090-15b5-11ec-3394-d99ce6470590
# ╠═54c19ee0-15b2-11ec-3cf0-b7e5cdc75d3a
# ╠═a658e3f0-15b5-11ec-0809-69ce7675d096
# ╠═75208430-15b2-11ec-0ae1-3724dade5f8e


text:
using EquationsOfState
using Test

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


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 78c5a32c-fef3-11eb-0560-6592672bbca5
begin
    using PlutoUI
end

# ╔═╡ 7feba529-29c5-4be5-aeba-00db5397c4da
md"""
# Data, Information, and Knowledge in this Short Course

Definitions are tricky things, and trying to define what is meant by _data_, _information_, and _knowledge_ is a slippery task.

>[T]he problems of definition are perhaps best exposed in a study in 2007 by the information scientist Chaim Zins that compiled definitions of _information_ (and _data_ and _knowledge_) offered by forty-five information scholars from sixteen countries. The compilation revealed 130 distinct notions, with different degrees of compatibility with one another, but no one capable of encompassing all thedifferent variations. (_Information: A Historical Companion_ p. x)

In this short course, we are going to largely be following a simplified view of the data-information-knowledge map provided by Luciano Floridi:

![Semantic concept of information](https://plato.stanford.edu/entries/information-semantic/figure1.jpg)

A useful set of definitions to start with is provided by Jefrey Pomerantz in _Metadata (MIT Press Essential Knowledge series)_.

- "Data, according to this view, is the raw stuff: what is collected by instrumentation or machinery."
- "Information...is data that has been processed into a form that may be consumed by a human being."
- "Knowledge...is what you know, information that you have internalized."

In what follows we will expand, clarify, hedge around these definitions.

## What is Data?

### Big Picture View of Data

> Data...is the raw stuff: what is collected by instrumentation or machinery. The stream of bits sent to Earth by a Mars Rover, for example, is data. The signal carried by radio frequencies between your phone and the local cell tower is data....
>
>
> Data is stuff. It is raw, unprocessed, possibly even untouched by human hands, unviewed by human eyes, un-thought-about by human minds. We’re not used to thinking of information objects in this way: we’re used to thinking of information objects as things like books, or files on our computer, things that have been deliberately created by humans, and where human understanding is an integral part of their creation. However, think about the stream of bits sent to Earth by the Mars Rover, or a book in, for example, [Lushootseed](https://en.wikipedia.org/wiki/Lushootseed) (or some other language that you neither speak nor read… and apologies if you actually do know Lushootseed): you may know that a stream of bits or a book in Lushootseed has some meaning embedded in it, but without some processing, that meaning is not accessible to you. **Data is potential information, analogous to potential energy: work is required to release it.** (Pomerantz, Jeffrey. *Metadata (MIT Press Essential Knowledge series)* (pp. 20-21). The MIT Press. Kindle Edition.)

### Data as Difference: _The Diaphoric Definition of Data (DDD)_

Ultimately this data "stuff" is some difference, some difference in some physical quantity whether it is a voltage level on a transitor in a computer representing "ON" or the pigment on a sheet of paper depicting an "a", depending on the context.

### Examples

#### Letters

The modern English alphabet consists of 26 lower-case and 26 upper-case letters. It is sometimes referred to as the Latin alphabet (for example in computer type standards), although it differs slightly from the classical Latin alphabet (Classical Latin did not have _J_, _U_, or _W_). In this context __a__ and __D__ would each be a _datum_ as they are different symbols.

#### Words

In a different context we may want to consider a word as a datum. So instead of viewing __BRIAN__ as the data set __{B__, __R__, __I__, __A__, __N}__, __BRIAN__ is a single datum that differs from __BRAIN__, for example, within the context of English

## Exercise

Can you think of different contexts where "__Brian Chapman__" would represent a data set with

- 13 elements?
- 12 elements?
- 2 elements?
- 1 element?

$(@bind myNameContext TextField((80,7)))
"""

# ╔═╡ 2acf7039-a934-43fa-b956-d581b5d49169
let myNameContextAns = md"""

- __13 elements:__ There are 13 characters (counting the space between "Brian" and "Chapman" as a character) in the string. I might not want to differentiate between letters and spaces if I was figuring out how to typeset the string for printing.
- __12 elements:__ There are 12 letters in the string. I may be interested in the letters and no tthe space if were counting letter frequency, for example trying to encrypt or decrypt a message.
- __2 elements:__ There are two words if a view the space as delimiting words.
- __1 element:__ In the context of names, there is one name in the string.
"""

HTML("My Answers$(html(myNameContextAns))")
    end

# ╔═╡ 2ccc8246-67cd-47ac-82a3-a06b23eb0a7d
md"""

## Example

Here are some actual data from my life:

```
141/73 6/3/2021
```

### What do you think these data and are what do they represent?

"""

# ╔═╡ ccd85c80-99d8-48c5-8a23-143a402d6ff0
md"""
#### At the most "atomistic" view we have the following data:

```
'1' '4' '1' '/' '7' '3' ' ' '6' '/' '3' '/' '2' '0' '2' '1'
```

#### Or we could view numbers as belonging together and spaces and '/' representing something else:

```
'141' '/' '73' ' ' '6' '/' '3' '/' '2021'
```

#### Or we could view spaces as separating pieces of data:

```
'141/73' '6/3/2021'

```
"""

# ╔═╡ 9064db1d-3031-4e43-943d-1e601821e44f
md"""### What are some of the ways we can interpret the `141/73` data?

#### Be creative!"""

# ╔═╡ 7a611c30-c745-4ade-881e-cec71e4f6973
md"""#### Type your answers here
$(@bind resp1 TextField((80,7)))
"""

# ╔═╡ 3b245fa2-1f83-45c3-9b3b-5fa62d002818
let answer1 = md"""

#### 141/73?

- $(parse(Int, "141", base=8) / parse(Int, "73", base=8 )) Assuming Base 8  Integer
- $(parse(Int, "141", base=8) / parse(Int, "73", base=16)) Assuming Base 16 Integer
- $(parse(Int, "141", base=8) / parse(Int, "73", base=10)) Assumint Base 10 Integer
- Systolic blood pressure of 141 and Diastolic blood pressure of 73
"""

HTML("My Answers$(html(answer1))")
    end

# ╔═╡ 9d4a1119-bbbe-4d26-af88-6789420ab030
md"""### What if I tell you it is health-related data?"""

# ╔═╡ 9a4bd43f-c2ae-4dc1-8b33-6a4422c4bcfb

# ╔═╡ 97c45761-4f20-4063-800f-ce9d79111466
md"""#### Type your answers here
$(@bind resp1b TextField((80,7)))
"""

# ╔═╡ 15193b6b-98b5-4cf6-bf65-3133ca734112
md"""#### What are some of the ways we can interpret the '6/3/2021' data?"""

# ╔═╡ bccd1a45-33aa-43a1-8f70-68cd9a40c41f
md"""#### Type your answers here
$(@bind resp1c TextField((80,7)))
"""

# ╔═╡ 4efcfeaa-eb86-4cf3-b4e9-60210278991a
let answer2 = md"""#### 6/3/2021?

- $(6/3/2021) Assuming base 10 and `/` indicates division
- 6 March 2021
- 3 June 2021
"""

HTML("My Answers$(html(answer2))")
    end

# ╔═╡ bd257e42-52b5-4fd6-9214-596617a5c10d
md"""## What are they actually?"""

# ╔═╡ d95f2476-0286-4aa3-ab25-c6d50b113c04
let answer4 = md"""
- 141/73
  - My maximum and minimum heart rate (beats per minute) on this bikeride

#### Anyone guess this as the physisological quantity?
"""

HTML("My Answers$(html(answer4))")
    end

# ╔═╡ 08349618-e3a3-4d99-a56d-e145544266e2
let answer3 = md"""
- June 3, 2021
  - The date of a bike ride I recorded on Strava

#### Anyone guess this date?

"""
HTML("My Answers$(html(answer3))")
    end

# ╔═╡ bc23a5b8-b489-4d5a-9fad-572e2c9aea4d
md"""
# Creating Information from Data

If data is the "stuff" that information is made of, and if "information" is data that has been processed for human (and I would add computing) how is this done?

Our answer to this question is we create information from data by __adding__ __metadata__ to the data. The data we are adding metadata to is known as __primary data__.

__Warning:__ Terminology can quickly become a mess.

## What is metadata?

>[W]hen metadata is doing its job well, it just fades into the background, unnoticed and nearly invisible. (Pomerantz, Jeffrey. _Metadata (MIT Press Essential Knowledge series)_ p. 1)

---------

>The most common—---and perhaps least useful---definition of metadata is that it is “data about data.” As catchy as this definition is, however, it is entirely ambiguous. First of all, what is data? And second, what does “about” mean? (Pomerantz, Jeffrey. _Metadata (MIT Press Essential Knowledge series_, p. 19)

Despite Dr. Pomerantz's concerns, metadata as data about data is probably an adequate definition for our purposes. And, as he says, it is so common that you are likely to hear it many, many times, and perhaps say it many times yourself. Nonetheless, for completness, I give you Dr. Pomerantz's preferred definition:

>Metadata Is a Statement about a Potentially Informative Object. (_Metadata (MIT Press Essential Knowledge series)_ p. 26).

### Examples

Here are some clips from the hospital bill from my birth. For each clip, we'll identify the primary data and associated meta as well as mention any metadata that seems to be missing.

#### Clip 1
![metadata 1](https://github.com/chapmanbe/raw-stuff/blob/main/MetaData/metadata1.png?raw=true)

- Primary Data: `F`
  - Metadata: `SEX`
  - Missing metadata: The meaning of `F`. We infer that `F` stands for `FEMALE`. We will later see examples of how sex is encoded in a more complex manner.
- Primary Data: `30`
    - Metadata: `PATIENT AGE`
    - Missing metadata: The units for age. We infer from the context (a maternity bill) that the age is in years. In a different context age could easily be in days or weeks. By common usage we assume that `30` is written in base 10 (per our example earlier).
----------

#### Clip 2

![metadata 2](https://github.com/chapmanbe/raw-stuff/blob/main/MetaData/metadata3.png?raw=true)

- Primary Data: `5.75`
    - Metadata:
        - `Pharmacy`
        - `Amount`
        -  `106`. This is interesting because this datum, which is metadata for `5.75` has its own metadata: `Pharmacy Code` to which if we look at the back of the bill we can see that `Pharmacy Code 106` means "Central nervous system drugs".
    - Missing Metadata: The units for amount.

![metadata 3](https://github.com/chapmanbe/raw-stuff/blob/main/MetaData/metadata4.png?raw=true)

"""

# ╔═╡ 7f2e9d8b-9085-4f65-b5e0-d45eb80de42d
md"""
# Knowledge

Pomerantz provided a starting definition of knowledge as "what you know, information that you have internalized." This definition has a number of relevant challenges.

The most well known definition of knowledge comes from the ancient Greek philosopher Plato who defined knowledge as __justified, true belief.__

How does Plato's definition compare to what Pomerantz provided and how do both it affect us in our Data To Knowledge work?

- __Knowledge must be justified__, that is you can explain why/how you know something. In the context of our use case, to "know" what we learn from the data entails being able to explain the model(s) we will develop in week 4.
- __Knowledge must be true,__ so we must have internalized true information and not misinformation or disinformation (see Floridi's figure at the top of this notebook. Perhaps more importantly, we should have an element of humility in what we say we know as medical/scientific knowledge is continually being revised or even reversed. As [Karl Popper](https://plato.stanford.edu/entries/popper/) famously pointed out, all scientific “knowledge” is tentative, something like “approximate temporary truth.”  /Popper
- __Knowledge is belief/internalized information__: This subjective view of knowledge, is something we will eventually be departing from as we will think of AI applications as having "knowledge."
"""

# ╔═╡ 4ea542ea-af63-4bf2-a242-ac35f1756bbe
md"""
## Detailed Blood Glucose Data Example

To wrap up this section, I'm sharing example data for a blood glucose measurement, as represented using HL7-FHIR standards.

What do you feel the primary data is here?

```Yaml
code:
  coding:
  - code: 15074-8
    display: Glucose [Moles/volume] in Blood
    system: http://loinc.org
effectivePeriod:
  start: '2013-04-02T09:30:10+01:00'
id: f001
identifier:
- system: http://www.bmc.nl/zorgportal/identifiers/observations
  use: official
  value: '6323'
interpretation:
- coding:
  - code: H
    display: High
    system: http://terminology.hl7.org/CodeSystem/v3-ObservationInterpretation
issued: '2013-04-03T15:30:10+01:00'
performer:
- display: A. Langeveld
  reference: Practitioner/f005
referenceRange:
- high:
    code: mmol/L
    system: http://unitsofmeasure.org
    unit: mmol/l
    value: 6.2
  low:
    code: mmol/L
    system: http://unitsofmeasure.org
    unit: mmol/l
    value: 3.1
resourceType: Observation
status: final
subject:
  display: P. van de Heuvel
  reference: Patient/f001
valueQuantity:
  code: mmol/L
  system: http://unitsofmeasure.org
  unit: mmol/l
  value: 6.3
```

In my view the primary data here is `6.3` and everything else is metadata to help make this `6.3` useable.

Notice the difference in metadata detail for this one data versus the metadata for the primary data in my birth bill. This is a nice lead-in to our next topic which is data collection as abstraction.

"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─7feba529-29c5-4be5-aeba-00db5397c4da
# ╟─2acf7039-a934-43fa-b956-d581b5d49169
# ╟─2ccc8246-67cd-47ac-82a3-a06b23eb0a7d
# ╟─ccd85c80-99d8-48c5-8a23-143a402d6ff0
# ╟─9064db1d-3031-4e43-943d-1e601821e44f
# ╟─7a611c30-c745-4ade-881e-cec71e4f6973
# ╟─3b245fa2-1f83-45c3-9b3b-5fa62d002818
# ╟─9d4a1119-bbbe-4d26-af88-6789420ab030
# ╟─9a4bd43f-c2ae-4dc1-8b33-6a4422c4bcfb
# ╟─97c45761-4f20-4063-800f-ce9d79111466
# ╟─15193b6b-98b5-4cf6-bf65-3133ca734112
# ╟─bccd1a45-33aa-43a1-8f70-68cd9a40c41f
# ╟─4efcfeaa-eb86-4cf3-b4e9-60210278991a
# ╟─bd257e42-52b5-4fd6-9214-596617a5c10d
# ╟─d95f2476-0286-4aa3-ab25-c6d50b113c04
# ╟─08349618-e3a3-4d99-a56d-e145544266e2
# ╟─bc23a5b8-b489-4d5a-9fad-572e2c9aea4d
# ╟─7f2e9d8b-9085-4f65-b5e0-d45eb80de42d
# ╟─4ea542ea-af63-4bf2-a242-ac35f1756bbe
# ╟─78c5a32c-fef3-11eb-0560-6592672bbca5
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


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

JLLWrappers.@generate_main_file_header("LAMMPS")
JLLWrappers.@generate_main_file("LAMMPS", UUID("5b3ab26d-9607-527c-88ea-8fe5ba57cafe"))
end  # module LAMMPS_jll


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using ACE, Random
using Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing
using JuLIP: evaluate, evaluate_d, evaluate!, evaluate_d!,
             alloc_temp, alloc_temp_d
using ACE: alloc_B, alloc_dB
using BenchmarkTools
using Juno

# (:Si, [:C, :O])
# N = 2:2:6
# maxdeg in [7, 12, 17]
degrees = Dict(2 => [7, 12, 17],
               4 => [7, 10, 13],
               6 => [7, 9, 11])
for species in (:Si, ), N = 2:2:6, maxdeg in degrees[N]
   r0 = 2.3
   rcut = 5.0
   trans = PolyTransform(1, r0)
   Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
   P1 = ACE.BasicPSH1pBasis(Pr; species = species)
   D = ACE.SparsePSHDegree()

   basis = ACE.PIBasis(P1, N, D, maxdeg, evaluator = :classic)
   basisdag = ACE.PIBasis(P1, N, D, maxdeg, evaluator = :dag)

   @info("species = $species; N = $N; length = $(length(basis))")

   for Nat in [5, 50]
      Rs, Zs, z0 = ACE.Random.rand_nhd(Nat, basis.basis1p.J,
                                         species)
      tmp = alloc_temp(basis, Rs, Zs, z0)
      tmpd = alloc_temp_d(basis, Rs, Zs, z0)
      tmpdag = alloc_temp(basisdag, Rs, Zs, z0)
      tmpdagd = alloc_temp_d(basisdag, Rs, Zs, z0)
      B = alloc_B(basis, Rs)
      dB = alloc_dB(basis, Rs)
      # ------------------------------
      println("   Nat = $Nat")
      print("   classic eval : ");
      @btime evaluate!($B, $tmp, $basis, $Rs, $Zs, $z0)
      print("     graph eval : ");
      @btime evaluate!($B, $tmpdag, $basisdag, $Rs, $Zs, $z0)
      print("   classic grad : ");
      @btime evaluate_d!($B, $dB, $tmpd, $basis, $Rs, $Zs, $z0)
      print("     graph grad : ")
      @btime evaluate_d!($B, $dB, $tmpdagd, $basisdag, $Rs, $Zs, $z0)
      println()
   end
end

# species = :Si
# N  = 4
# maxdeg = 7
# r0 = 2.3
# rcut = 5.0
# trans = PolyTransform(1, r0)
# Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
# P1 = ACE.BasicPSH1pBasis(Pr; species = species)
# D = ACE.SparsePSHDegree()
#
# Rs, Zs, z0 = ACE.Random.rand_nhd(30, Pr, species)
# basis = ACE.PIBasis(P1, N, D, maxdeg, evaluator = :classic)
# basisdag = ACE.PIBasis(P1, N, D, maxdeg, evaluator = :dag)
# tmp = alloc_temp(basis, Rs, Zs, z0)
# tmpd = alloc_temp_d(basis, Rs, Zs, z0)
# tmpdag = alloc_temp(basisdag, Rs, Zs, z0)
# tmpdagd = alloc_temp_d(basisdag, Rs, Zs, z0)
# B = alloc_B(basis, Rs)
# dB = alloc_dB(basis, Rs)

# @btime evaluate_d!($B, $dB, $tmpdagd, $basisdag, $Rs, $Zs, $z0)
# @btime evaluate_d!($B, $dB, $tmpdagd, $basisdag, $Rs, $Zs, $z0)

#
# function runN(N, f, args...)
#    for n = 1:N
#       f(args...)
#    end
#    return nothing
# end
#
# runN(10, evaluate_d!, B, dB, tmpd, basis, Rs, Zs, z0)
#
# Juno.@profiler runN(10_000, evaluate_d!, B, dB, tmpd, basis, Rs, Zs, z0)


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

@testset "GeoMean" begin
    bridged_mock = MOIB.Constraint.GeoMean{Float64}(mock)

    MOIT.basic_constraint_tests(
        bridged_mock, config,
        include = [(F, S)
                   for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64},
                             MOI.VectorQuadraticFunction{Float64}]
                   for S in [MOI.GeometricMeanCone]])

    @testset "geomean1test" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [ones(4); 2; √2; √2])
        MOIT.geomean1vtest(bridged_mock, config)
        MOIT.geomean1ftest(bridged_mock, config)

        @testset "Test mock model" begin
            var_names = ["t", "x", "y", "z", "x21", "x11", "x12"]
            MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
            lessthan = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
            @test length(lessthan) == 2
            MOI.set.(mock, MOI.ConstraintName(), lessthan, ["lessthan1", "lessthan2"])
            rsoc = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone}())
            @test length(rsoc) == 3
            MOI.set.(mock, MOI.ConstraintName(), rsoc, ["rsoc21", "rsoc11", "rsoc12"])

            s = """
            variables: t, x, y, z, x11, x12, x21
            lessthan1: t + -0.5 * x21 in MathOptInterface.LessThan(0.0)
            lessthan2: x + y + z in MathOptInterface.LessThan(3.0)
            rsoc11: [1.0x, y, x11] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc12: [z, 0.5 * x21, x12] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc21: [1.0x11, x12, x21] in MathOptInterface.RotatedSecondOrderCone(3)
            maxobjective: t
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(mock, model, var_names, ["lessthan1", "lessthan2", "rsoc11", "rsoc12", "rsoc21"])
        end

        @testset "Test bridged model" begin
            var_names = ["t", "x", "y", "z"]
            MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
            lessthan = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
            @test length(lessthan) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), lessthan[1], "lessthan")
            geomean = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}())
            @test length(geomean) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")

            s = """
            variables: t, x, y, z
            lessthan: x + y + z in MathOptInterface.LessThan(3.0)
            geomean: [1.0t, x, y, z] in MathOptInterface.GeometricMeanCone(4)
            maxobjective: t
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(bridged_mock, model, var_names, ["lessthan", "geomean"])
        end

        # Dual is not yet implemented for GeoMean bridge
        ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}()))
        test_delete_bridge(bridged_mock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),
                                                (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0)))
    end

    @testset "geomean2test" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, vcat(ones(10), 4.0, fill(2 * √2, 2), fill(2, 4), fill(√2, 8)))
        MOIT.geomean2vtest(bridged_mock, config)
        MOIT.geomean2ftest(bridged_mock, config)

        @testset "Test mock model" begin
            x_names = ["x$(i)" for i in 1:9]
            l1_names = ["x1$(i)" for i in 1:8]
            l2_names = ["x2$(i)" for i in 1:4]
            l3_names = ["x3$(i)" for i in 1:2]
            var_names = vcat("t", x_names, "x41", l3_names, l2_names, l1_names)
            MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
            lessthan = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
            @test length(lessthan) == 1
            MOI.set(mock, MOI.ConstraintName(), lessthan[1], "lessthan")
            equalto = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}())
            @test length(equalto) == 9
            equalto_names = ["equalto$(i)" for i in 1:9]
            MOI.set.(mock, MOI.ConstraintName(), equalto, equalto_names)
            rsoc = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone}())
            @test length(rsoc) == 15
            rsoc_names = vcat("rsoc41", ["rsoc3$(i)" for i in 1:2], ["rsoc2$(i)" for i in 1:4], ["rsoc1$(i)" for i in 1:8])
            MOI.set.(mock, MOI.ConstraintName(), rsoc, rsoc_names)

            varnames_str = vcat("t", [", " * v for v in var_names[2:end]])
            s = """
            variables: $(varnames_str...)
            lessthan: t + -0.25 * x41 in MathOptInterface.LessThan(0.0)
            equalto1: 1.0x1 in MathOptInterface.EqualTo(1.0)
            equalto2: 1.0x2 in MathOptInterface.EqualTo(1.0)
            equalto3: 1.0x3 in MathOptInterface.EqualTo(1.0)
            equalto4: 1.0x4 in MathOptInterface.EqualTo(1.0)
            equalto5: 1.0x5 in MathOptInterface.EqualTo(1.0)
            equalto6: 1.0x6 in MathOptInterface.EqualTo(1.0)
            equalto7: 1.0x7 in MathOptInterface.EqualTo(1.0)
            equalto8: 1.0x8 in MathOptInterface.EqualTo(1.0)
            equalto9: 1.0x9 in MathOptInterface.EqualTo(1.0)
            rsoc41: [1.0x31, x32, x41] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc31: [1.0x21, x22, x31] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc32: [1.0x23, x24, x32] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc21: [1.0x11, x12, x21] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc22: [1.0x13, x14, x22] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc23: [1.0x15, x16, x23] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc24: [1.0x17, x18, x24] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc11: [1.0x1, x2, x11] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc12: [1.0x3, x4, x12] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc13: [1.0x5, x6, x13] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc14: [1.0x7, x8, x14] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc15: [1.0x9, 0.25 * x41, x15] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc16: [0.25 * x41, 0.25 * x41, x16] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc17: [0.25 * x41, 0.25 * x41, x17] in MathOptInterface.RotatedSecondOrderCone(3)
            rsoc18: [0.25 * x41, 0.25 * x41, x18] in MathOptInterface.RotatedSecondOrderCone(3)
            maxobjective: t
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(mock, model, var_names, vcat("lessthan", equalto_names, rsoc_names))
        end

        @testset "Test bridged model" begin
            var_names = vcat("t", ["x$(i)" for i in 1:9])
            MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
            equalto = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}())
            @test length(equalto) == 9
            equalto_names = ["equalto$(i)" for i in 1:9]
            MOI.set.(bridged_mock, MOI.ConstraintName(), equalto, equalto_names)
            geomean = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}())
            @test length(geomean) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")

            s = """
            variables: t, x1, x2, x3, x4, x5, x6, x7, x8, x9
            equalto1: 1.0x1 in MathOptInterface.EqualTo(1.0)
            equalto2: 1.0x2 in MathOptInterface.EqualTo(1.0)
            equalto3: 1.0x3 in MathOptInterface.EqualTo(1.0)
            equalto4: 1.0x4 in MathOptInterface.EqualTo(1.0)
            equalto5: 1.0x5 in MathOptInterface.EqualTo(1.0)
            equalto6: 1.0x6 in MathOptInterface.EqualTo(1.0)
            equalto7: 1.0x7 in MathOptInterface.EqualTo(1.0)
            equalto8: 1.0x8 in MathOptInterface.EqualTo(1.0)
            equalto9: 1.0x9 in MathOptInterface.EqualTo(1.0)
            geomean: [1.0t, x1, x2, x3, x4, x5, x6, x7, x8, x9] in MathOptInterface.GeometricMeanCone(10)
            maxobjective: t
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(bridged_mock, model, var_names, vcat(equalto_names, "geomean"))
        end

        # Dual is not yet implemented for GeoMean bridge
        ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}()))
        test_delete_bridge(bridged_mock, ci, 10, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
                                                (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0)))
    end

    @testset "geomean3test" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, 2 * ones(2))
        MOIT.geomean3vtest(bridged_mock, config)
        MOIT.geomean3ftest(bridged_mock, config)

        @testset "Test mock model" begin
            MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), ["t", "x"])
            lessthan = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
            @test length(lessthan) == 2
            MOI.set.(mock, MOI.ConstraintName(), lessthan, ["lessthan1", "lessthan2"])
            rsoc = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone}())
            @test isempty(rsoc)
            nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
            @test length(nonneg) == 1
            MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")

            s = """
            variables: t, x
            lessthan1: t + -1.0x in MathOptInterface.LessThan(0.0)
            lessthan2: 1.0x in MathOptInterface.LessThan(2.0)
            nonneg: [1.0x] in MathOptInterface.Nonnegatives(1)
            maxobjective: 2.0 * t
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(mock, model, ["t", "x"], ["lessthan1", "lessthan2", "nonneg"])
        end

        @testset "Test bridged model" begin
            MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), ["t", "x"])
            lessthan = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
            @test length(lessthan) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), lessthan[1], "lessthan")
            geomean = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}())
            @test length(geomean) == 1
            MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")

            s = """
            variables: t, x
            lessthan: 1.0x in MathOptInterface.LessThan(2.0)
            geomean: [1.0t, x] in MathOptInterface.GeometricMeanCone(2)
            maxobjective: 2.0 * t
            """
            model = MOIU.Model{Float64}()
            MOIU.loadfromstring!(model, s)
            MOIU.test_models_equal(bridged_mock, model, ["t", "x"], ["lessthan", "geomean"])
        end

        # Dual is not yet implemented for GeoMean bridge
        ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}()))
        test_delete_bridge(bridged_mock, ci, 2, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),
                                                (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0)))
    end

end


text:
import CartesianGrids: convective_derivative!, divergence!, grad!, curl!

for f in [:divergence!, :grad!, :curl!]
  @eval $f(a,b,sys::ILMSystem) = $f(a,b,sys.base_cache)
end

_scale_derivative!(w,cache::BasicILMCache{N,IndexScaling}) where {N} = w
_scale_derivative!(w,cache::BasicILMCache{N,GridScaling}) where {N} = w ./= cellsize(cache)

_scale_laplacian!(w,cache::BasicILMCache{N,IndexScaling}) where {N} = w
_scale_laplacian!(w,cache::BasicILMCache{N,GridScaling}) where {N} = w # ./= cellsize(cache)^2

_scale_inverse_laplacian!(w,cache::BasicILMCache{N,IndexScaling}) where {N} = w
_scale_inverse_laplacian!(w,cache::BasicILMCache{N,GridScaling}) where {N} = w  #.*= cellsize(cache)^2

"""
    divergence!(p::Nodes{Primal},v::Edges{Primal},cache::BasicILMCache)
    divergence!(p::Nodes{Primal},v::Edges{Primal},sys::ILMSystem)

Compute the discrete divergence of `v` and return it in `p`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function divergence!(p::Nodes{Primal,NX,NY},q::Edges{Primal,NX,NY},cache::BasicILMCache) where {NX,NY}
    _unscaled_divergence!(p,q,cache)
    _scale_derivative!(p,cache)
end

_unscaled_divergence!(p,q,cache::BasicILMCache) = (fill!(p,0.0); divergence!(p,q))

"""
    grad!(v::Edges{Primal},p::Nodes{Primal},cache::BasicILMCache)
    grad!(v::Edges{Primal},p::Nodes{Primal},sys::ILMSystem)

Compute the discrete gradient of `p` and return it in `v`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function grad!(q::Edges{Primal,NX,NY},p::Nodes{Primal,NX,NY},cache::BasicILMCache) where {NX,NY}
    _unscaled_grad!(q,p,cache)
    _scale_derivative!(q,cache)
end

_unscaled_grad!(q,p,cache::BasicILMCache) = (fill!(q,0.0); grad!(q,p))

"""
    curl!(v::Edges{Primal},s::Nodes{Dual},cache::BasicILMCache)
    curl!(v::Edges{Primal},s::Nodes{Dual},sys::ILMSystem)

Compute the discrete curl of `s` and return it in `v`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function curl!(q::Edges{Primal,NX,NY},s::Nodes{Dual,NX,NY},cache::BasicILMCache) where {NX,NY}
    _unscaled_curl!(q,s,cache)
    _scale_derivative!(q,cache)
end

_unscaled_curl!(q::Edges,s::Nodes,cache::BasicILMCache) = (fill!(q,0.0); curl!(q,s))

"""
    curl!(w::Nodes{Dual},v::Edges{Primal},cache::BasicILMCache)
    curl!(w::Nodes{Dual},v::Edges{Primal},sys::ILMSystem)

Compute the discrete curl of `v` and return it in `w`, scaling it
by the grid spacing if `cache` (or `sys`) is of `GridScaling` type, or leaving it
as a simple differencing if `cache` (or `sys`) is of `IndexScaling` type.
"""
function curl!(w::Nodes{Dual,NX,NY},q::Edges{Primal,NX,NY},cache::BasicILMCache) where {NX,NY}
    _unscaled_curl!(w,q,cache)
    _scale_derivative!(w,cache)
end

_unscaled_curl!(w::Nodes,q::Edges,cache::BasicILMCache) = (fill!(w,0.0); curl!(w,q))

"""
    laplacian!(w::GridData,s::GridData,sys::ILMSystem)

Compute the Laplacian of grid data `s`, and divide the result
by unity or by the grid cell size, depending on whether `sys` has `IndexScaling` or `GridScaling`,
respectively, and return the result as `w`.
"""
laplacian!(w,s,sys::ILMSystem) = laplacian!(w,s,sys.base_cache)

"""
    laplacian!(w::GridData,s::GridData,cache::BasicILMCache)

Compute the Laplacian of grid data `s`, and divide the result
by unity or by the grid cell size, depending on whether `cache` has `IndexScaling` or `GridScaling`,
respectively, and return the result as `w`.
"""
function laplacian!(w::GridData,s::GridData,cache::BasicILMCache)
    @unpack L = cache
    _unscaled_laplacian!(w,s,L)
    _scale_laplacian!(w,cache)
end

_unscaled_laplacian!(w::GridData,s::GridData,L::CartesianGrids.Laplacian) = w .= L*s

"""
    inverse_laplacian!(w::GridData,sys::ILMSystem)

Compute the in-place inverse Laplacian of grid data `w`, and multiply the result
by unity or by the grid cell size, depending on whether `sys` has `IndexScaling` or `GridScaling`,
respectively.
"""
inverse_laplacian!(w,sys::ILMSystem) = inverse_laplacian!(w,sys.base_cache)

"""
    inverse_laplacian!(w::GridData,cache::BasicILMCache)

Compute the in-place inverse Laplacian of grid data `w`, and multiply the result
by unity or by the grid cell size, depending on whether `cache` has `IndexScaling` or `GridScaling`,
respectively.
"""
function inverse_laplacian!(w::GridData,cache::BasicILMCache)
    @unpack L = cache
    _unscaled_inverse_laplacian!(w,L)
    _scale_inverse_laplacian!(w,cache)
end

_unscaled_inverse_laplacian!(w::GridData,L::CartesianGrids.Laplacian) = w .= L\w

###  Convective derivative of velocity-like data by itself ###

"""
    ConvectiveDerivativeCache(dv::EdgeGradient)

Create a cache (a subtype of [`AbstractExtraILMCache`](@ref)) for computing
the convective derivative, using `dv` to define the cache data.
"""
struct ConvectiveDerivativeCache{VTT} <: AbstractExtraILMCache
   vt1_cache :: VTT
   vt2_cache :: VTT
   vt3_cache :: VTT
   ConvectiveDerivativeCache(vt::EdgeGradient) = new{typeof(vt)}(similar(vt),similar(vt),similar(vt))
end

"""
    ConvectiveDerivativeCache(cache::BasicILMCache)

Create a cache for computing the convective derivative, based on
the basic ILM cache `cache`.
"""
ConvectiveDerivativeCache(cache::BasicILMCache) = ConvectiveDerivativeCache(EdgeGradient(Primal,cache.gdata_cache))

"""
    convective_derivative(v::Edges,base_cache::BasicILMCache)

Compute the convective derivative of `v`, i.e., ``v\\cdot \\nabla v``. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`.
"""
function convective_derivative(u::Edges{Primal},base_cache::BasicILMCache)
    extra_cache = ConvectiveDerivativeCache(EdgeGradient(Primal,u))
    udu = similar(u)
    convective_derivative!(udu,u,base_cache,extra_cache)
end

"""
    convective_derivative!(vdv::Edges,v::Edges,base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)

Compute the convective derivative of `v`, i.e., ``v\\cdot \\nabla v``, and
return the result in `vdv`. The result is either divided by unity or
the grid cell size depending on whether `base_cache` is of type `IndexScaling`
or `GridScaling`. This version of the method uses `extra_cache` of type
[`ConvectiveDerivativeCache`](@ref).
"""
function convective_derivative!(udu::Edges{Primal},u::Edges{Primal},base_cache::BasicILMCache,extra_cache::ConvectiveDerivativeCache)
    _unscaled_convective_derivative!(udu,u,extra_cache)
    _scale_derivative!(udu,base_cache)
end

function _unscaled_convective_derivative!(udu::Edges{Primal},u::Edges{Primal},extra_cache::ConvectiveDerivativeCache)
    @unpack vt1_cache, vt2_cache, vt3_cache = extra_cache

    fill!(vt1_cache,0.0)
    grid_interpolate!(vt1_cache,u)
    CartesianGrids.transpose!(vt2_cache,vt1_cache)
    fill!(vt1_cache,0.0)
    grad!(vt1_cache,u)
    product!(vt3_cache,vt2_cache,vt1_cache)
    fill!(udu,0.0)
    grid_interpolate!(udu,vt3_cache)
    udu
end


text:
"""
    Spin{T1, T2}

Spin parameters of an asteroid

# Fields

## Spin pole @ equatorial coordinate
- `α`  # Right ascension (RA)
- `δ`  # Declination (Dec)

## Spin pole @ ecliptic coordinate
- `λ`  # Ecliptic longitude
- `β`  # Ecliptic latitude

## Other parameters
- `T`  # Spin period [sec]
- `ω`  # Angular velocity [rad/sec]
- `ŝ`  # Spin pole direction (normalized)
- `ε`  # Obliquity
- `γ`  # Vernal equinox lon. from the direction of perihelion
"""
struct Spin{T1, T2}
    α::T1
    δ::T1

    λ::T1
    β::T1

    T::T1
    ω::T1
    ŝ::T2
    ε::T1

    γ::T1
end

function Base.show(io::IO, spin::Spin)
    println(io, "Spin parameters")
    println(io, "---------------")

    println("Right ascension (RA) : α = ", rad2deg(spin.α), " [deg]")
    println("Declination (Dec)    : δ = ", rad2deg(spin.δ), " [deg]")
    println("Ecliptic longitude   : λ = ", rad2deg(spin.λ), " [deg]")
    println("Ecliptic latitude    : β = ", rad2deg(spin.β), " [deg]")
    println("Obliquity            : ε = ", rad2deg(spin.ε), " [deg]")
    println("Spin period          : P = ", spin.T / 3600,   " [hours]")
    println("Spin rate            : ω = ", spin.ω,          " [rad/sec]")
    println("Vernal equinox lon.  : γ = ", rad2deg(spin.γ), " [deg]")
    println("                           (longitude from the periheion direction)")
end

function setSpinParams(params, orbit)

    if haskey(params, :α) && haskey(params, :δ)
        α = deg2rad(params[:α])
        δ = deg2rad(params[:δ])
        λ, β = eq2ec(α, δ)
    elseif haskey(params, :λ) && haskey(params, :β)
        λ = deg2rad(params[:λ])
        β = deg2rad(params[:β])
        α, δ = ec2eq(λ, β)
    else
        println("Give spin pole direction.")
    end

    ŝ = getSpinNormal(λ, β, orbit)  # Orbital plane frame
    ε = acos(ŝ[3])
    γ = getVernalEquinox(ŝ)

    T = params[:T] * 3600
    ω = 2π / T

    return Spin(α, δ, λ, β, T, ω, ŝ, ε, γ)
end

function setSpinParams(T::Float64, ŝ::AbstractVector, orbit)

    T *= 3600
    ω = 2π / T
    ε = acos(ŝ[3])
    γ = getVernalEquinox(ŝ)

    λ = - atan(ŝ[2], ŝ[1])
    β = asin(ŝ[3])
    α, δ = ec2eq(λ, β)

    return Spin(α, δ, λ, β, T, ω, ŝ, ε, γ)
end

##############################################################################
##############################################################################

"""
Get a spin pole direction in the orbital plane frame

# Paramters
- λ : ecliptic longitude
- β : ecliptic latitude
- orbit : orbital elements
# Return
- ŝ : spin pole direction (normalized)
"""
function getSpinNormal(λ, β, orbit)
    ŝ = SA_F64[cos(β) * cos(λ), - cos(β) * sin(λ), sin(β)]  # inertial frame
    ŝ = inertia_to_orbit(ŝ, orbit)                          # orbital plane frame
end

"""
Get a unit vector ê⟂1 and ê⟂2
"""
function getSpinUnits(spin::Spin)
    N̂ = SA_F64[0, 0, 1]
    ê1 = (spin.ŝ * cos(spin.ε) - N̂) / sin(spin.ε)
    ê2 = (spin.ŝ × N̂) / sin(spin.ε)

    ê1, ê2
end

"""
Get a unit vector ê⟂1 only
"""
function getSpinUnit1(spin::Spin)
    N̂ = SA_F64[0, 0, 1]
    ê1 = (spin.ŝ * cos(spin.ε) - N̂) / sin(spin.ε)
end

"""
Get a unit vector ê⟂2 only
"""
function getSpinUnit2(spin::Spin)
    N̂ = SA_F64[0, 0, 1]
    ê2 = (spin.ŝ × N̂) / sin(spin.ε)
end

"""
Get a longitude of the vernal equinox with respect to the perihelion
"""
function getVernalEquinox(spin::Spin)
    ê2 = getSpinUnit2(spin)
    vernal_equinox = atan(ê2[2], ê2[1]) + π
end

function getVernalEquinox(ŝ)
    N̂ = SA_F64[0, 0, 1]
    ê2 = ŝ × N̂
    vernal_equinox = atan(ê2[2], ê2[1]) + π
end

"""
Get a longitude of the autumnal equinox with respect to the perihelion
"""
function getAutumnalEquinox(spin)
    ê2 = getSpinUnit2(spin)
    autumnal_equinox = atan(ê2[2], ê2[1])
    autumnal_equinox < 0 && (autumnal_equinox += 2π)
    autumnal_equinox
end

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
# As is discussed in the paper, you can see that the FD derivatives here are so
# inaccurate that Ipopt can't even converge. The FD Hessians are worse than
# garbage and completely unusable without very high-order and adaptive methods,
# which are so brutal to performance that they aren't even worth considering.

using Ipopt, Serialization

include("shared.jl")
include("gradient.jl")
include("ipopt_helpers.jl")

function fitter(objects_initval)
  ((case, gradfun, hesfun, maxiter), ini) = objects_initval
  is_bfgs = in(case, (:FD_BFGS, :AD_BFGS))
  println("Case: $case")
  cache = Vector{Float64}[]
  box_l = is_bfgs ? [1e-2, 1e-2, 0.25] : [0.0, 0.0, 0.25]
  prob  = createProblem(3, box_l, fill(1e22, (3,)), 0,
                        Float64[], Float64[], 0, div(3*4, 2),
                        _p->caching_nll(_p, cache),
                        (args...)->nothing,
                        gradfun,
                        (args...)->nothing,
                        (x,m,r,c,o,l,v)->ipopt_hessian(x,m,r,c,o,l,v,
                                                       hesfun,Function[],0))
  addOption(prob, "tol", 1e-5)
  addOption(prob, "max_iter", maxiter)
  if is_bfgs
    addOption(prob, "hessian_approximation", "limited-memory")
    addOption(prob, "nlp_scaling_method", "none")
  end
  prob.x = deepcopy(ini) # for safety to avoid weird persistent pointer games.
  try
    @time status = solveProblem(prob)
    _h = is_bfgs ? fil(NaN, 3, 3) : hesfun(prob.x)
    return (prob, status, case, deepcopy(prob.x), _nll(prob.x), _h, cache, ini)
  catch er
    println("\n\nOptimization failed with error $er\n\n")
    return (prob, :FAIL_OPT_ERR, case, deepcopy(prob.x), NaN, fill(NaN, 3, 3), cache, ini)
  end
end

const cases = ((:FD_FISH, grad_fd!, fishfd, 100),
               (:AD_FISH, grad_ad!, fishad, 100),
               (:FD_HESS, grad_fd!, hessfd, 100),
               (:AD_HESS, grad_ad!, _nllh,  100),
               (:FD_BFGS, grad_fd!, no_hessian, 100),
               (:AD_BFGS, grad_ad!, no_hessian, 100))

const inits = (ones(3), [1.0, 0.1, 2.0])
const test_settings = vec(collect(Iterators.product(cases, inits)))

if !isinteractive()
  const res   = map(fitter, test_settings)
  serialize("fit_results.serialized", res)
end


text:
export SoftScore

"""
    SoftScore(vs::Vector{I}, ss::Vector{Float64})

Return a `LogScore` of the log values in `ss` vector for
the associated keys in `vs`.
"""
function SoftScore(vs::Vector{I}, ss::Vector{Float64}) where I
    return LogScore(vs, [log(s) for s in ss])
end

"""
    SoftScore(scores::Dict{I,Float64})

Return a `LogScore` of the keys and log values in `score`.
"""
function SoftScore(scores::Dict{I,Float64}) where I
    d = Dict([(k,log(v)) for (k,v) in scores])
    return LogScore(d)
end


text:
__precompile__()

module SimpleDiffEq

using Reexport
@reexport using DiffEqBase
using StaticArrays

include("functionmap.jl")
include("euler_maruyama.jl")

end # module


text:
module NeuralArithmetic

using ChainRulesCore
using Flux
using Requires
using Zygote

using Flux: glorot_uniform, destructure

include("nau.jl")
include("nmu.jl")
include("npu.jl")

function __init__()
    @require DiffEqFlux="aae7a2af-3d4f-5e19-a356-7da93b79d9d0" include("fast_nau.jl")
    @require DiffEqFlux="aae7a2af-3d4f-5e19-a356-7da93b79d9d0" include("fast_npu.jl")
end

include("nac.jl")
include("nalu.jl")
include("inalu.jl")

end # module


text:
const INIT_CONSTS_CODE = []

const INIT_MODULES = Dict(
    :pybuiltinsmodule => "builtins",
    :pysysmodule => "sys",
    :pyosmodule => "os",
    :pynumbersmodule => "numbers",
    :pydatetimemodule => "datetime",
    :pycollectionsabcmodule => "collections.abc",
)

for (j,m) in INIT_MODULES
    @eval const $j = pynew()
    push!(INIT_CONSTS_CODE, :(pycopy!($j, pyimport($m))))
end

const INIT_ATTRS = Dict(
    :pyfractiontype => "fractions" => "Fraction",
    :pydatetype => "datetime" => "date",
    :pytimetype => "datetime" => "time",
    :pydatetimetype => "datetime" => "datetime",
    :pytimedeltatype => "datetime" => "timedelta",
)

for (j,k) in INIT_ATTRS
    @eval const $j = pynew()
    push!(INIT_CONSTS_CODE, :(pycopy!($j, pyimport($k))))
end

const BUILTINS = Set([
    # consts
    :True,
    :False,
    :NotImplemented,
    :None,
    :Ellipsis,
    # classes/functions
    :abs,
    :all,
    :any,
    :ascii,
    :bin,
    :bool,
    :bytes,
    :bytearray,
    :callable,
    :chr,
    :classmethod,
    :compile,
    :complex,
    :delattr,
    :dict,
    :dir,
    :divmod,
    :enumerate,
    :eval,
    :exec,
    :filter,
    :float,
    :format,
    :frozenset,
    :getattr,
    :globals,
    :hasattr,
    :hash,
    :help,
    :hex,
    :id,
    :input,
    :int,
    :isinstance,
    :issubclass,
    :iter,
    :len,
    :list,
    :locals,
    :map,
    :max,
    :memoryview,
    :min,
    :next,
    :object,
    :oct,
    :open,
    :ord,
    :pow,
    :print,
    :property,
    :range,
    :repr,
    :reversed,
    :round,
    :set,
    :setattr,
    :slice,
    :sorted,
    :staticmethod,
    :str,
    :sum,
    :super,
    :tuple,
    :type,
    :vars,
    :zip,
    # exceptions
    :BaseException,
    :Exception,
    :StopIteration,
    :GeneratorExit,
    :ArithmeticError,
    :LookupError,
    :AssertionError,
    :AttributeError,
    :BufferError,
    :EOFError,
    :FloatingPointError,
    :OSError,
    :ImportError,
    :IndexError,
    :KeyError,
    :KeyboardInterrupt,
    :MemoryError,
    :NameError,
    :OverflowError,
    :RuntimeError,
    :RecursionError,
    :NotImplementedError,
    :SyntaxError,
    :IndentationError,
    :TabError,
    :ReferenceError,
    :SystemError,
    :SystemExit,
    :TypeError,
    :UnboundLocalError,
    :UnicodeError,
    :UnicodeEncodeError,
    :UnicodeDecodeError,
    :UnicodeTranslateError,
    :ValueError,
    :ZeroDivisionError,
    :BlockingIOError,
    :BrokenPipeError,
    :ChildProcessError,
    :ConnectionError,
    :ConnectionAbortedError,
    :ConnectionRefusedError,
    :FileExistsError,
    :FileNotFoundError,
    :InterruptedError,
    :IsADirectoryError,
    :NotADirectoryError,
    :PermissionError,
    :ProcessLookupError,
    :TimeoutError,
    :EnvironmentError,
    :IOError,
    :Warning,
    :UserWarning,
    :DeprecationWarning,
    :PendingDeprecationWarning,
    :SyntaxWarning,
    :RuntimeWarning,
    :FutureWarning,
    :ImportWarning,
    :UnicodeWarning,
    :BytesWarning,
    :ResourceWarning,
])

@eval baremodule pybuiltins
    $([:(const $k = $pynew()) for k in BUILTINS]...)
end
export pybuiltins

for k in BUILTINS
    push!(INIT_CONSTS_CODE, :(pycopy!(pybuiltins.$k, pybuiltinsmodule.$k)))
end

@eval function init_consts()
    $(INIT_CONSTS_CODE...)
    return
end


text:
struct StandardNMDT <: AbstractRelaxation
    num_layers::Int
end

function reformulate_unary_quadratic_term!(
    model::MOI.ModelLike,
    reformulation::Reformulation{StandardNMDT},
    _x::VI,
    _y::VI,
    _x_l::Float64,
    _x_u::Float64,
)
    x, y = _map_to_interval_if_needed(model, _x, _y, _x_l, _x_u)
    x_l, x_u = 0.0, 1.0

    L = reformulation.relaxation.num_layers
    α_v, α_c =
        MOI.add_constrained_variables(model, [MOI.ZeroOne() for i in 1:L])
    Δx_v, Δx_c =
        MOI.add_constrained_variable(model, MOI.Interval(0.0, 2.0^(-L)))
    u_v = MOI.add_variables(model, L)
    Δz_v = MOI.add_variable(model)

    α = [SV(α_v[i]) for i in 1:L]
    u = [SV(u_v[i]) for i in 1:L]
    x = 1.0 * SV(x)
    y = 1.0 * SV(y)
    Δx = SV(Δx_v)
    Δz = SV(Δz_v)

    MOI.add_constraint(
        model,
        x - Δx - sum(2.0^(-i) * α[i] for i in 1:L),
        MOI.EqualTo(0.0),
    )
    MOI.add_constraint(
        model,
        y - Δz - sum(2.0^(-i) * u[i] for i in 1:L),
        MOI.EqualTo(0.0),
    )
    for i in 1:L
        a, b, c = x, α[i], u[i]
        a_l, a_u = x_l, x_u
        b_l, b_u = 0.0, 1.0
        MOIU.normalize_and_add_constraint(
            model,
            c - (a_u * b + a * b_u - a_u * b_u),
            MOI.GreaterThan(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            c - (a_l * b + a * b_l - a_l * b_l),
            MOI.GreaterThan(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            c - (a_u * b + a * b_l - a_u * b_l),
            MOI.LessThan(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            c - (a_l * b + a * b_u - a_l * b_u),
            MOI.LessThan(0.0),
        )
    end
    a, b, c = x, Δx, Δz
    a_l, a_u = x_l, x_u
    b_l, b_u = 0.0, 2.0^(-L)
    MOIU.normalize_and_add_constraint(
        model,
        c - (a_u * b + a * b_u - a_u * b_u),
        MOI.GreaterThan(0.0),
    )
    MOIU.normalize_and_add_constraint(
        model,
        c - (a_l * b + a * b_l - a_l * b_l),
        MOI.GreaterThan(0.0),
    )
    MOIU.normalize_and_add_constraint(
        model,
        c - (a_u * b + a * b_l - a_u * b_l),
        MOI.LessThan(0.0),
    )
    MOIU.normalize_and_add_constraint(
        model,
        c - (a_l * b + a * b_u - a_l * b_u),
        MOI.LessThan(0.0),
    )
    return α_v
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

High-level pnml labels are expected to have <text> and <structure> elements.
This concrete type is for "unclaimed" labels in a high-level petri net.
Some "claimed" `HLAnnotation` labels are [`Condition`](@ref),
[`Declaration`](@ref), [`HLMarking`](@ref), [`HLInscription`](@ref).
"""
struct HLLabel <: HLAnnotation
    text::Maybe{String}
    structure::Maybe{Structure}
    com::ObjectCommon    #TODO labels, toolinfos, graphics
    xml::XMLNode
    #TODO validate in constructor: must have text or structure
end

#------------------------------------------------------------------------
"""
$(TYPEDEF)
$(TYPEDFIELDS)

Name is for display, possibly in a tool specific way.
"""
struct Name <: AbstractLabel
    text::String
    graphics::Maybe{Graphics}
    tools::Maybe{Vector{ToolInfo}}
end

Name(; text::AbstractString = "", graphics=nothing, tools=nothing) =
    Name(text, graphics, tools)


text:
#! /usr/bin/julia

# Rosetta Code, Knuth's algorithm S

function makesofn(n::Int)
    buf = Any[]
    i = 0
    function sofn(item)
        i += 1
        if i <= n
            push!(buf, item)
        else
            j = rand(1:i)
            if j <= n
                buf[j] = item
            end
        end
        return buf
    end
    return sofn
end

nhist = zeros(Int, 10)

for i in 1:10^5
    kas = makesofn(3)
    for j in 0:8
        kas(j)
    end
    for k in kas(9)
        nhist[k+1] += 1
    end
end

println("Simulating sof3(0:9) 100000 times:")
for (i, c) in enumerate(nhist)
    println(@sprintf "   %2d => %5d" i-1 c)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using LShapedSolvers
using LinearAlgebra
using Test
using Logging
using JuMP
using StochasticPrograms
using GLPKMathProgInterface
using Ipopt

τ = 1e-5
reference_solver = GLPKSolverLP()
qp_solver = IpoptSolver(print_level=0)
lsolvers = [(LShapedSolver(reference_solver, log=false),"L-Shaped"),
            (LShapedSolver(qp_solver, subsolver=reference_solver, crash=Crash.EVP(), regularization=:rd, autotune=true, log=false),"RD L-Shaped"),
            (LShapedSolver(reference_solver, crash=Crash.EVP(), regularization=:rd, autotune=true, log=false, linearize=true),"Linearized RD L-Shaped"),
            (LShapedSolver(reference_solver, crash=Crash.EVP(), regularization=:tr, autotune=true, log=false),"TR L-Shaped"),
            (LShapedSolver(reference_solver, projectionsolver=qp_solver, regularization=:lv, log=false),"Leveled L-Shaped"),
            (LShapedSolver(reference_solver, regularization=:lv, log=false, linearize=true),"Linearized Leveled L-Shaped")]
problems = Vector{Tuple{<:StochasticProgram,String}}()
@info "Loading test problems..."
@info "Loading simple..."
include("simple.jl")
@info "Loading farmer..."
include("farmer.jl")
@info "Loading infeasible..."
include("infeasible.jl")
@info "Test problems loaded. Starting test sequence."

@testset "Sequential solvers" begin
    @testset "Simple problems" begin
        @testset "$lsname Solver: $name" for (lsolver,lsname) in lsolvers, (sp,name) in problems
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            optimize!(sp, solver=lsolver)
            @test abs(optimal_value(sp) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Bundled $lsname Solver: $name" for (lsolver,lsname) in lsolvers, (sp,name) in problems
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            add_params!(lsolver, bundle=2)
            optimize!(sp, solver=lsolver)
            @test abs(optimal_value(sp) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
    end
    @testset "Infeasible problem" begin
        @testset "$lsname Solver: Feasibility cuts" for (lsolver,lsname) in lsolvers
            optimize!(infeasible, solver=reference_solver)
            x̄ = optimal_decision(infeasible)
            Q̄ = optimal_value(infeasible)
            with_logger(NullLogger()) do
                @test optimize!(infeasible, solver=lsolver) == :Infeasible
            end
            add_params!(lsolver, checkfeas=true)
            optimize!(infeasible, solver=lsolver)
            @test abs(optimal_value(infeasible) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(infeasible) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Bundled $lsname Solver: Feasibility cuts" for (lsolver,lsname) in lsolvers
            optimize!(infeasible, solver=reference_solver)
            x̄ = optimal_decision(infeasible)
            Q̄ = optimal_value(infeasible)
            add_params!(lsolver, checkfeas=false, bundle=2)
            with_logger(NullLogger()) do
                @test optimize!(infeasible, solver=lsolver) == :Infeasible
            end
            add_params!(lsolver, checkfeas=true, bundle=2)
            optimize!(infeasible, solver=lsolver)
            @test abs(optimal_value(infeasible) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(infeasible) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
    end
end

@info "Starting distributed tests..."

include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))
disttestfile = joinpath(@__DIR__, "run_dtests.jl")
push!(test_exeflags.exec,"--color=yes")
cmd = `$test_exename $test_exeflags $disttestfile`

if !success(pipeline(cmd; stdout=stdout, stderr=stderr)) && ccall(:jl_running_on_valgrind,Cint,()) == 0
    error("Distributed test failed, cmd : $cmd")
end


text:
export edistance

#######################################################################################
# %% (Mostly) space agnostic helper functions
#######################################################################################

"""
    edistance(a, b, model::ABM)

Return the euclidean distance between `a` and `b` (either agents or agent positions),
respecting periodic boundary conditions (if in use). Works with any space where it makes
sense: currently `GridSpace` and `ContinuousSpace`.
"""
edistance(
    a::A,
    b::B,
    model::ABM{C,<:Union{ContinuousSpace,GridSpace}},
) where {A<:AbstractAgent,B<:AbstractAgent,C} = edistance(a.pos, b.pos, model)

function edistance(
    a::ValidPos,
    b::ValidPos,
    model::ABM{A,<:Union{ContinuousSpace{D,false},GridSpace{D,false}}},
) where {A,D}
    sqrt(sum(abs2.(a .- b)))
end

function edistance(
    p1::ValidPos,
    p2::ValidPos,
    model::ABM{A,<:ContinuousSpace{D,true}},
) where {A,D}
    total = 0.0
    for (a, b, d) in zip(p1, p2, model.space.extent)
        delta = abs(b - a)
        if delta > d - delta
            delta = d - delta
        end
        total += delta^2
    end
    sqrt(total)
end

function edistance(p1::ValidPos, p2::ValidPos, model::ABM{A,<:GridSpace{D,true}}) where {A,D}
    total = 0.0
    for (a, b, d) in zip(p1, p2, size(model.space))
        delta = abs(b - a)
        if delta > d - delta
            delta = d - delta
        end
        total += delta^2
    end
    sqrt(total)
end