text:
function tocategoricaldf(df::AbstractDataFrame; categoricalexcept)
    df = typeof(df) <: DataFrame ? copy(df) : DataFrame(df)
    categorical = setdiff(names(df), [categoricalexcept])
    categorical!(df, categorical)
end

copysubarray(x::SubArray) = copy(x)
copysubarray(x) = x
function splitplotanova(df, locals; data, kwargs...)
    aoverror(locals, copysubarray(getproperty(df, data)))
end
function splitplotanova(df, locals::Base.RefValue;
                        betweenfactors, withinfactors, id, data,
                        categoricalexcept = string.(data))
    formula = term(data) ~ predictor_product([betweenfactors; withinfactors])
    df = tocategoricaldf(df; categoricalexcept)
    res, locals[] = aoverror(formula, id, withinfactors, df; return_locals = true)
    res
end
anova(df, formula, locals) = aov(locals, copysubarray(getproperty(df, formula.lhs.sym)))
function anova(df, formula, locals::Base.RefValue)
    df = tocategoricaldf(df; categoricalexcept = string(formula.lhs.sym))
    res, locals[] = aov(formula, df; return_locals = true)
    res
end

abstract type AbstractSimpleTest end
struct UTest <: AbstractSimpleTest
    u::Float64
    p::Float64
    n::Int
end
function utest(x, y)
    res = MannWhitneyUTest(x, y)
    UTest(res.U, HypothesisTests.pvalue(res), length(x) + length(y))
end
function Base.show(io::IO, t::UTest)
    println(io, "Mann-Whitney U test")
    println(io, "u = $(t.u) p = $(t.p) n = $(t.n)")
end
struct TTest <: AbstractSimpleTest
    t::Float64
    p::Float64
    df::Int
end
struct NoTest <: AbstractSimpleTest
    comment::String
    p::Float64
end
# this is a hack for Correia07_exp2 where the content can be a string
NoTest(comment) = NoTest(comment, 1.)
function ttest(x::AbstractVector, y::AbstractVector, μ0::Real=0;
               paired = false, twosided = true)
    res = paired ? OneSampleTTest(x, y, μ0) : EqualVarianceTTest(x, y, μ0)
    TTest(res.t, HypothesisTests.pvalue(res, tail = twosided ? :both : :right), res.df)
end
function Base.show(io::IO, t::TTest)
    println(io, "Student's t-test")
    println(io, "t = $(t.t) p = $(t.p) DF = $(t.df)")
end
function is_significantly_bigger(df;
                                 conditionkey::Symbol,
                                 condbigger, condsmaller,
                                 valuekey::Symbol,
                                 paired = true)
    x = DataFramesMeta.where(df, d -> getproperty(d, conditionkey) .== condbigger)
    y = DataFramesMeta.where(df, d -> getproperty(d, conditionkey) .== condsmaller)
    ttest(getproperty(x, valuekey), getproperty(y, valuekey), paired = paired, twosided = false)
end
function isnot_significantly_different(df; conditionkey::Symbol,
                                       condition1, condition2,
                                       valuekey::Symbol,
                                       paired = true)
    x = DataFramesMeta.where(df, d -> getproperty(d, conditionkey) .== condition2)
    y = DataFramesMeta.where(df, d -> getproperty(d, conditionkey) .== condition1)
    ttest(getproperty(x, valuekey), getproperty(y, valuekey), paired = paired, twosided = true)
end
struct ChisqTest <: AbstractSimpleTest
    χ²::Float64
    p::Float64
    df::Int
end
function mergeanddropzerorows(x, y)
    tmp = Int[]
    for (xi, yi) in zip(x, y)
        xi == yi == 0 && continue
        push!(tmp, xi); push!(tmp, yi)
    end
    reshape(tmp, 2, :)
end
chisqtest(x::Vector{Int}, y::Vector{Int}) = chisqtest(mergeanddropzerorows(x, y))
function chisqtest(x::Array{Int, 2})
    χ² = HypothesisTests.ChisqTest(x)
    ChisqTest(χ².stat, HypothesisTests.pvalue(χ²), χ².df)
end
function Base.show(io::IO, c::ChisqTest)
    println(io, "χ²-test")
    println(io, "χ²($(c.df)) = $(c.χ²) p = $(c.p)")
end
struct FisherExactTest <: AbstractSimpleTest
    table::NTuple{4, Int}
    p::Float64
end
function Base.show(io::IO, f::FisherExactTest)
    println(io, "Fisher exact test")
    println(io, "$(f.table[1])|$(f.table[2]) vs $(f.table[3])|$(f.table[4])  p = $(f.p)")
end
function fisherexacttest(df; cond1, cond2, data)
    table = combine(groupby(df, [cond1, cond2]), d -> sum(getproperty(d, data)))
    tabletuple = tuple(table.x1...)
    if tabletuple[1] == tabletuple[3] == 0 || tabletuple[2] == tabletuple[4] == 0
        FisherExactTest(tabletuple, 1.) # this is a hack to avoid errors
    else
        p = HypothesisTests.pvalue(HypothesisTests.FisherExactTest(tabletuple...))
        FisherExactTest(tabletuple, p)
    end
end
struct BinomialTest <: AbstractSimpleTest
    x::Int
    n::Int
    p::Float64
end
function Base.show(io::IO, f::BinomialTest)
    println(io, "Binomial test")
    println(io, "x/n = $(f.x)/$(f.n), p = $(f.p)")
end
function binomialtest(v::Vector)
    x = sum(v)
    n = length(v)
    p = HypothesisTests.pvalue(HypothesisTests.BinomialTest(x, n))
    BinomialTest(x, n, p)
end

struct Test{T,Ta}
    id::String
    content::T
    target::Ta
    pvalues::Vector{Float64}
end
function Test(id, content, target)
    Test(id, content, target, pvalues(target.accessor, content, target))
end
function Test(id, content, target::Dict)
    Test(id, content, target[id])
end
function Base.show(io::IO, t::Test)
    printstyled(io, "Test ID:  $(t.id)\n", bold = true)
    pvalues = t.target.pvalues
    for (k, p) in zip(t.target.keys, pvalues)
        color = true ? :green : :red
        printstyled(io, "$(rpad(k, 25)) p = $(round(p, digits = 4))\n", color = color)
    end
    println(io, t.content)
end

pvalues(::Any, c::AbstractSimpleTest, ::Any) = [c.p]
df(::Any, c::AbstractSimpleTest, ::Any) = ""
pvalues(a, c::DataFrame, ::Any) = c.p[a]
df(a, c::DataFrame, ::Any) = map(x -> "$(c.DF[x]),$(c.DF[end])", a)
pvalues(a, c::Vector{DataFrame}, ::Any) = vcat([c[i].p[idx] for (i, idx) in a]...)
df(a, c::Vector{DataFrame}, x::Any) = vcat([df([idx], c[i], x) for (i, idx) in a]...)
function pvalues(a::Base.RefValue, c::DataFrame, t)
    a[] = map(name -> findfirst(x -> x == name, c.names), t.keys)
    pvalues(a[], c, t)
end
function accessor_tuple(k, c)
    for (i, d) in enumerate(c)
       idx = findfirst(x -> x == k, d.names)
       idx === nothing || return (i, idx)
    end
end
function pvalues(a::Base.RefValue, c::Vector{DataFrame}, t)
    a[] = [accessor_tuple(k, c) for k in t.keys]
    pvalues(a[], c, t)
end

struct TestCollection{Tt}
    tests::Tt
    pvalues::Vector{Float64}
end
function Base.show(io::IO, ts::TestCollection)
    for test in ts.tests
        Base.show(io, test)
        println()
    end
end
TestCollection(test::Test) = TestCollection([test])
function TestCollection(tests::Vector{<:Test})
    TestCollection(tests, vcat([t.pvalues for t in tests]...))
end

function significancecode(p; levels = SIGLEVELS)
    i = one(Float64)
    for v in levels
        p < v && return i
        i += 1.
    end
    i
end
const SIGLEVELS = [.001, .01, .05, .1]
function setrsiglevels!()
    empty!(SIGLEVELS)
    append!(SIGLEVELS, [.001, .01, .05, .1])
end
function setbinarysiglevels!()
    empty!(SIGLEVELS)
    push!(SIGLEVELS, .05)
end


text:
import Adapt: adapt
import .Tracker: IdSet

children(x) = ()
mapchildren(f, x) = x

children(x::Tuple) = x
mapchildren(f, x::Tuple) = map(f, x)

function treelike(T, fs = fieldnames(T))
  @eval current_module() begin
    Flux.children(x::$T) = ($([:(x.$f) for f in fs]...),)
    Flux.mapchildren(f, x::$T) = $T(f.($children(x))...)
  end
end

isleaf(x) = isempty(children(x))

function mapleaves(f, x; cache = ObjectIdDict())
  haskey(cache, x) && return cache[x]
  cache[x] = isleaf(x) ? f(x) : mapchildren(x -> mapleaves(f, x, cache = cache), x)
end

function prefor(f, x; seen = IdSet())
  x ∈ seen && return
  f(x)
  foreach(x -> prefor(f, x, seen = seen), children(x))
  return
end

function params(m)
  ps = []
  prefor(p ->
    Tracker.istracked(p) && Tracker.isleaf(p) &&
      !any(p′ -> p′ === p, ps) && push!(ps, p),
    m)
  return ps
end

params(m...) = params(m)

function loadparams!(m, xs)
  for (p, x) in zip(params(m), xs)
    size(p) == size(x) ||
      error("Expected param size $(size(p)), got $(size(x))")
    copy!(data(p), data(x))
  end
end

# CPU/GPU movement conveniences

cpu(m) = mapleaves(x -> adapt(Array, x), m)

gpu_adaptor = identity

@require CuArrays begin
  global gpu_adaptor = CuArrays.cu
end

gpu(x) = mapleaves(gpu_adaptor, x)

ext: jl


text:
# This file extends / overrides some of the Julia basic capability for linear algebra to add useful things for FE
# Most of the time these are "in-place" operations that Julia doesn't normally have an interface for, but we need.
# Also: these functions help us have a common interface with PETSc (which has more in-place capability than Julia)

"""
    Set all entries to zero
"""
function zero!(v)
    fill!(v, 0.)
end

"""
    v[i] += u
"""
function plusEquals!(v::Array, u::Array, i)
    v[i] += u
end

"""
    mat[i,j] += u
"""
function plusEquals!(mat::Matrix, u::Matrix, i, j)
    mat[i,j] += u
end

"""
    Zero rows in the matrix
"""
function zeroRows!(mat::Matrix, rows::Array)
    mat[rows,:] = 0
end

"""
    No-op (no need to assemble Julia objects)
"""
function assemble!(::Any)
end

"""
    Pass through
"""
function serializeToZero(array::Array)
    return array
end


text:
module ModelCheck

import
    GLM,
    Gadfly,
    Distributions,
    Random,
    StatsModels,
    StatsBase,
    Statistics,
    LinearAlgebra,
    LinRegOutliers,
    Colors

# Separate .jl files for each function (have to manually add file, otherwise error)

# Function #1: Calculate deviance residuals
include("calc_deviance_resids.jl")
export calc_deviance_resids

# Function #2: Plot quantile-quantile plot
include("plot_qqplot.jl")
export plot_qqplot

# Function #3: Plot residual vs fitted plot
include("plot_residual_fitted.jl")
export plot_residual_fitted

# Function #4: Plot scale-location plot
include("plot_scale_location.jl")
export plot_scale_location

# Function #5: Plot all diagnostics plot
include("plot_model_check.jl")
export plot_model_check

end


text:
import Base: *, -

"""
    *(f, g::AbstractGP)

Produce an `AbstractGP` `h` satisfying to `h(x) = f(x) * g(x)`, for some deterministic
function `f`.

If `f isa Real`, then `h(x) = f * g(x)`.
"""
*(f, g::AbstractGP) = CompositeGP((*, f, g), g.gpc)
*(f::AbstractGP, g) = CompositeGP((*, g, f), f.gpc)
*(f::AbstractGP, g::AbstractGP) = throw(ArgumentError("Cannot multiply two GPs together."))

const prod_args{Tf} = Tuple{typeof(*), Tf, <:AbstractGP}

#
# Scale by a function
#

mean((_, σ, g)::prod_args, x::AV) = σ.(x) .* mean(g, x)

function cov((_, σ, g)::prod_args, x::AV)
    σx = σ.(x)
    return σx .* cov(g, x) .* σx'
end
var((_, σ, g)::prod_args, x::AV) = σ.(x).^2 .* var(g, x)

function cov((_, σ, g)::prod_args, x::AV, x′::AV)
    return σ.(x) .* cov(g, x, x′) .* σ.(x′)'
end
function var((_, σ, g)::prod_args, x::AV, x′::AV)
    return σ.(x) .* var(g, x, x′) .* σ.(x′)
end

cov((_, σ, f)::prod_args, f′::AbstractGP, x::AV, x′::AV) = σ.(x) .* cov(f, f′, x, x′)
cov(f::AbstractGP, (_, σ, f′)::prod_args, x::AV, x′::AV) = cov(f, f′, x, x′) .* (σ.(x′))'

function var((_, σ, f)::prod_args, f′::AbstractGP, x::AV, x′::AV)
    return σ.(x) .* var(f, f′, x, x′)
end
function var(f::AbstractGP, (_, σ, f′)::prod_args, x::AV, x′::AV)
    return var(f, f′, x, x′) .* σ.(x′)
end

#
# Scale by a constant
#

mean((_, σ, g)::prod_args{<:Real}, x::AV) = σ .* mean(g, x)

cov((_, σ, g)::prod_args{<:Real}, x::AV) = (σ^2) .* cov(g, x)
var((_, σ, g)::prod_args{<:Real}, x::AV) = (σ^2) .* var(g, x)

cov((_, σ, g)::prod_args{<:Real}, x::AV, x′::AV) = (σ^2) .* cov(g, x, x′)
var((_, σ, g)::prod_args{<:Real}, x::AV, x′::AV) = (σ^2) .* var(g, x, x′)

cov((_, σ, f)::prod_args{<:Real}, f′::AbstractGP, x::AV, x′::AV) = σ .* cov(f, f′, x, x′)
cov(f::AbstractGP, (_, σ, f′)::prod_args{<:Real}, x::AV, x′::AV) = cov(f, f′, x, x′) .* σ

function var((_, σ, f)::prod_args{<:Real}, f′::AbstractGP, x::AV, x′::AV)
    return σ .* var(f, f′, x, x′)
end
function var(f::AbstractGP, (_, σ, f′)::prod_args{<:Real}, x::AV, x′::AV)
    return var(f, f′, x, x′) .* σ
end

# Use multiplication to define the negation of a GP
-(f::AbstractGP) = (-1) * f


text:
@testset "Property interface" begin
    nt = (a=1,b=2.0,c=false)
    @test @inferred(getproperty(:b)(nt)) === 2.0

    @test @inferred(TypedTables.getproperties(:b)(nt)) === (b = 2.0,)

    c1 = @Compute($b)
    @test c1 isa TypedTables.GetProperty
    @test @inferred(c1(nt)) === 2.0

    c2 = @Compute(2*$b)
    @test @inferred(c2(nt)) === 4.0

    c3 = @Compute($a + 2*$b)
    @test @inferred(c3(nt)) === 5.0

    s1 = @Select(a)
    @test s1 isa TypedTables.GetProperties
    @test @inferred(s1(nt)) === (a = 1,)

    s2 = @Select(c,b,a)
    @test s1 isa TypedTables.GetProperties
    @test @inferred(s2(nt)) === (c = false, b = 2.0, a = 1)

    s3 = @Select(a, c = $c, d = $a + 2*$b)
    @test @inferred(s3(nt)) === (a = 1, c = false, d = 5.0)
end


text:
using Catalyst, Unitful, Test
const MT = ModelingToolkit

@parameters α [unit=u"μM/s"] β [unit=u"s"^(-1)] γ [unit=u"μM*s"^(-1)]
@variables t [unit=u"s"] A(t) [unit=u"μM"] B(t) [unit=u"μM"] C(t) [unit=u"μM"]
rxs = [Reaction(α, nothing, [A]),
       Reaction(β, [A], [B]),
       Reaction(γ, [A,B], [B], [1,1], [2])]
rs = ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester"))
@test_nowarn ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester"))

odeunit = u"μM/s"
#jumpunit = u"s^(-1)"
for rx in reactions(rs)
    @test MT.get_unit(oderatelaw(rx)) == odeunit

    # we don't currently convert units, so they will be the same as for ODEs
    @test MT.get_unit(jumpratelaw(rx)) == odeunit
end

@test_nowarn convert(ODESystem,rs)
@test_nowarn convert(SDESystem,rs)
@test_nowarn convert(JumpSystem,rs)

@parameters β [unit=u"M"]
rxs = [Reaction(α, nothing, [A]),
       Reaction(β, [A], [B]),
       Reaction(γ, [A,B], [B], [1,1], [2])]
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester")))) == false

@parameters β [unit=u"s"^(-1)]
@variables B(t) [unit=u"M"]
rxs = [Reaction(α, nothing, [A]),
       Reaction(β, [A], [B]),
       Reaction(γ, [A,B], [B], [1,1], [2])]
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(ReactionSystem(rxs, t, [A,B,C], [α,β,γ], name=Symbol("unittester")))) == false

@variables B(t) [unit=u"μM"] D(t) [unit=u"M"]
badrx1 = Reaction(α, [A], [D])
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx1)) == false
badrx2 = Reaction(α, [A], [B,D])
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx2)) == false
badrx3 = Reaction(α, [A,D], [B])
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx3)) == false
badrx4 = Reaction(α + β, [A], [B])
@test (@test_logs (:warn, ) match_mode=:any ModelingToolkit.validate(badrx4)) == false


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

"""
    UpROOT.testdatadir()::String

Path to some test data installed by UpROOT.jl for testing and experimentation
purposes.

Use `UpROOT.testfiles` to get `Dict` of all test files.
"""
testdatadir() = joinpath(dirname(@__DIR__), "deps", "testdata")

const _testfile_keys = [
    "leaflist",
    "HZZ",
    "hepdata-example",
    "nesteddirs"
]

_testfilename(k::AbstractString) = joinpath(testdatadir(), "$k.root")

const testfiles = Dict{String,String}()


text:
module History

using Gtk, Gtk.ShortNames
using Maple
using ..Ahorn

include("history.jl")

struct RoomSnapshot <: Snapshot
    description::String
    room::Maple.Room
    layers::Array{String, 1}

    RoomSnapshot(description::String, room::Maple.Room, layers::Array{String, 1}=String[]) = new(description, deepcopy(room), layers)
end

struct SelectionSnapshot <: Snapshot
    description::String
    room::Maple.Room
    selections::Set{Ahorn.SelectedObject}

    SelectionSnapshot(description::String, room::Maple.Room, selections::Set{Ahorn.SelectedObject}) = new(description, room, deepcopy(selections))
end

struct MultiSnapshot <: Snapshot
    description::String
    snapshots::Array{Snapshot, 1}
end

function MapSnapshot(description::String, map::Maple.Map)
    snapshots = Snapshot[]

    for room in map.rooms
        push!(snapshots, RoomSnapshot("Map backup", room))
    end

    selections = getToolSelections()
    if Ahorn.loadedState.room !== nothing && isa(selections, Set)
        push!(snapshots, SelectionSnapshot("Map backup", Ahorn.loadedState.room, selections))
    end

    return MultiSnapshot(description, snapshots)
end

function restoreSnapshot!(map::Maple.Map, snapshot::SelectionSnapshot)
    return Ahorn.eventToModule(Ahorn.currentTool, "setSelections", map, snapshot.room, snapshot.selections)
end

function restoreSnapshot!(map::Maple.Map, snapshot::RoomSnapshot)
    room = Maple.getRoomByName(map, snapshot.room.name)

    if isa(room, Maple.Room)
        for n in fieldnames(typeof(room))
            setfield!(room, n, deepcopy(getfield(snapshot.room, n)))
        end

    else
        room = deepcopy(snapshot.room)
        push!(map.rooms, room)
    end

    dr = Ahorn.getDrawableRoom(map, room)
    for layer in dr.layers
        if layer.name in snapshot.layers || isempty(snapshot.layers)
            layer.redraw = true
        end
    end

    Ahorn.draw(Ahorn.canvas)
end

function restoreSnapshot!(map::Maple.Map, snapshot::MultiSnapshot)
    for s in snapshot.snapshots
        restoreSnapshot!(map, s)
    end
end

function getToolSelections()
    selectionRes = Ahorn.eventToModule(Ahorn.currentTool, "getSelections")
    if isa(selectionRes, Set)
        return deepcopy(selectionRes)
    end

    return Set{Ahorn.SelectedObject}()
end

const historyTimelines = Dict{Maple.Map, HistoryTimeline}()

currentMap() = Ahorn.loadedState.map

function getHistory(map::Maple.Map)
    if !haskey(historyTimelines, map)
        historyTimelines[map] = HistoryTimeline()
    end

    return historyTimelines[map]
end

function undo!(map::Maple.Map=currentMap())
    res = false
    Ahorn.eventToModule(Ahorn.currentTool, "beforeUndo", map)

    history = getHistory(map)
    if history.index > 0
        if history.index == length(history.snapshots)
            push!(history.snapshots, MapSnapshot("Map backup", map))
        end

        history.skip = true
        snapshot = pop!(history)
        @Ahorn.catchall restoreSnapshot!(map, snapshot)
        res = true
    end

    Ahorn.eventToModule(Ahorn.currentTool, "afterUndo", map)

    return res
end

function redo!(map::Maple.Map=currentMap())
    res = false
    Ahorn.eventToModule(Ahorn.currentTool, "beforeRedo", map)

    history = getHistory(map)

    if history.skip
        history.skip = false
        history.index += 1
    end

    if history.index < length(history.snapshots)
        history.index += 1
        snapshot = history.snapshots[history.index]

        if history.index == length(history.snapshots)
            pop!(history.snapshots)
            history.index -= 1
        end

        @Ahorn.catchall restoreSnapshot!(map, snapshot)
        res = true
    end

    Ahorn.eventToModule(Ahorn.currentTool, "afterRedo", map)

    return res
end

function undo!(widget::Ahorn.MenuItemsTypes)
    if !undo!()
        info_dialog("Cannot undo.\nStart of history.", Ahorn.window)
    end
end

function redo!(widget::Ahorn.MenuItemsTypes)
    if !redo!()
        info_dialog("Cannot redo.\nEnd of history.", Ahorn.window)
    end
end

function addSnapshot!(snapshot::Snapshot, map::Maple.Map=currentMap())
    history = getHistory(map)
    push!(history, snapshot)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################
"""
Example considering a double single feeder section taken from:

> Industrial Energy System Availability Management by T. Van Acker, pg. 92 (2020)

Power is transmitted from sources S₁/S₂ to users U₁/U₂ using two independent
single feeders (Fig. 1)

      n₁     fdr 1    n₃
    S₁ ⋅-□----///----□-⋅ U₁

    S₂ ⋅-□----///----□-⋅ U₂
      n₂     fdr 2    n₄
Fig. 1: Double single feeder section
"""

# load pkgs
using Unitful
using MultiStateSystems

# initialize the network corresponding with the double single feeder section
ntwᵈˢᶠ = Network()

# add the user, source and components to the network ntwᵈˢᶠ
add_users!(ntwᵈˢᶠ, node = [3,4])
add_sources!(ntwᵈˢᶠ, node = [1,2])
add_components!(ntwᵈˢᶠ, edge = [(1,3),(2,4)],
                        name = ["fdr 1", "fdr 2"],
                        std  = [solvedSTD(power = [0.0u"MW", 1.0u"MW"],
                                          prob  = [0.1,0.9]),
                                solvedSTD(power = [0.0u"MW", 1.0u"MW"],
                                          prob  = [0.1,0.9])])

# return the network
return ntwᵈˢᶠ


text:
@testset "promotion" begin
    for (mr,br,sr) in ((OneToMRange(10),         OneTo(10),              OneToSRange(10)),
                       (UnitMRange(1, 10),       UnitRange(1, 10),       UnitSRange(1, 10)),
                       (StepMRange(1, 1, 10),    StepRange(1, 1, 10),    StepSRange(1, 1, 10)),
                       (LinMRange(1, 10, 10),    LinRange(1, 10, 10),    LinSRange(1, 10, 10)),
                       (StepMRangeLen(1, 1, 10), StepRangeLen(1, 1, 10), StepSRangeLen(1, 1, 10))
                      )
        @testset "$(typeof(br)) lower range" begin
            @test last(promote(br, sr)) == br
            @test last(promote(sr, br)) == br
            @test last(promote(br, mr)) == mr
            @test last(promote(mr, br)) == mr
        end
        for (mr2,br2,sr2) in ((OneToMRange(10),         OneTo(10),              OneToSRange(10)),
                              (UnitMRange(1, 10),       UnitRange(1, 10),       UnitSRange(1, 10)),
                              (StepMRange(1, 1, 10),    StepRange(1, 1, 10),    StepSRange(1, 1, 10)),
                              (LinMRange(1, 10, 10),    LinRange(1, 10, 10),    LinSRange(1, 10, 10)),
                              (StepMRangeLen(1, 1, 10), StepRangeLen(1, 1, 10), StepSRangeLen(1, 1, 10))
                             )
            if br != br2
                @testset "$(typeof(br)) & $(typeof(br2))" begin
                    br3 = last(promote(br, br2))
                    @test eltype(last(promote(mr, mr2))) == eltype(br3)
                    @test eltype(last(promote(mr2, mr))) == eltype(br3)

                    @test eltype(last(promote(sr, sr2))) == eltype(br3)
                    @test eltype(last(promote(sr2, sr))) == eltype(br3)
                end
            end
        end
    end
end


text:
function convexHullSIMP!( points::Array{Tuple{T,T,T},1}; typ=UInt16 ) where {T<:Real}
	println("debugging")

    @assert length(points) >= 4

    # By now, the total length of faces is 2*length(points), this should provide
	# more than enough space for the all convex hull faces
    faces = Array{Tuple{typ,typ,typ},1}(undef,length(points)*2);
    norms = Array{Tuple{Float32,Float32,Float32},1}(undef,length(points)*2);
	visibleFaces = Array{typ,1}(undef,length(points)*2);

	# Making sure the 4 first points points are not the same, colinear or coplanar
	# These will form the initial tetrahedron that nucleates the convex hull
	fixSameColinearCoplanar!( points )

    faces[1] = ( 1, 2, 3 );
    faces[2] = ( 2, 3, 4 );
    faces[3] = ( 3, 4, 1 );
    faces[4] = ( 4, 1, 2 );
    lastFace = 4;

	center = ( points[1] .+ points[2] .+ points[3] .+ points[4] ) ./ 4;
    for fidx in 1:4
		norms[fidx] = faceNormalCenterSIMP( points, faces[fidx], center )
    end

	@inbounds for pidx in 5:length(points)

		new_point = points[pidx];

		numVisibleFaces = typ(0);

		for fidx in 1:lastFace;
			# vector from a point in the current face to new_point
			vp = new_point .- points[faces[fidx][1]];

			# new_point is outside/visible (positive dot product) from the current face
			if sign( sum( vp .* norms[ fidx ] ) ) == 1
				numVisibleFaces += typ(1);
				visibleFaces[numVisibleFaces] = fidx;
			end
		end

		# After adding the index of visible faces to "visibleFaces", we can proceed to
		# construct new faces over these visible faces to the new_point
		lastFace = addNewFacesSIMP!( points, faces, norms, center, lastFace,
						             pidx, visibleFaces, numVisibleFaces );

	end # @inbounds for pidx in 5:length(points)

    return faces[1:lastFace], norms[1:lastFace]
end

"""
	needs to be tested --> tested in TEST 1
"""
@inline function faceNormalCenterSIMP( points, face, ctr )
    v = points[face[2]] .- points[face[1]];
    w = points[face[3]] .- points[face[1]];
	n = ( v[2]*w[3] - v[3]*w[2], v[3]*w[1] - v[1]*w[3], v[1]*w[2] - v[2]*w[1] )
	mag = sqrt( sum( n .* n ) );
	n = n ./ mag
	p = ctr .- points[face[1]];
	return ( sign( sum( p .* n ) ) == -1 ) ? n : -1 .* n;
end

"""
	needs to be tested
"""
@inline function addNewFacesSIMP!( points,
                               faces::Array{Tuple{T,T,T},1},
							   norms::Array{Tuple{Float32,Float32,Float32},1},
							   center,
							   lastFace,
							   newpointindex,
							   visibleFaces::Array{T,1},
                               numVisibleFaces
                             ) where {T<:Integer}

	"""
	   1-. Isolating non-repeated edges among visibleFaces. These form the boundary
	       between visible faces and non-visible faces.
	"""
	boundary   = Array{Tuple{T,T},1}(undef, numVisibleFaces*3);
	isRepeated = zeros( Bool, numVisibleFaces*3 );
	lidx, nrep = 1, 0;

	@inbounds for face in faces[ visibleFaces[1:numVisibleFaces] ]
		edges = [(face[1],face[2]), (face[2],face[3]), (face[3],face[1])];
		for (p1, p2) in edges
			boundary[lidx] = ( p1, p2 )
			for bidx in 1:lidx-1
				boundary_edge = boundary[bidx]
				if (p1,p2) == boundary_edge || (p2,p1) == boundary_edge
					isRepeated[lidx] = true;
					isRepeated[bidx] = true;
					nrep += 2;
				end
			end # each previous boundary --> check that the current edge has not already been added
		lidx += 1
		end # each of the 3 edges --> add it to "boundary" and check if it is not repeated
	end # each visible face --> go through its 3 edges

	"""
		2-. From each unique edge in the boundary -> build a new face with the correct neighbouring data
	"""
	numNewFaces = lidx - 1 - nrep;
	change      = numNewFaces - numVisibleFaces;

	#( numNewFaces < numVisibleFaces ) && ( println( "Delete faces" ) );

	lidx = 0;
	for bidx in 1:length(boundary)
		if !(isRepeated[bidx])
			lidx += 1
			newfaceidx = ( lidx <= numVisibleFaces ) ? visibleFaces[lidx] : lastFace + lidx - numVisibleFaces;
			faces[newfaceidx] = ( boundary[bidx][1], boundary[bidx][2], newpointindex );
			norms[newfaceidx] = faceNormalCenterSIMP( points, ( boundary[bidx][1], boundary[bidx][2], newpointindex ), center )
		end
	end

	if ( change < 0 )
		cont = 0
		for start in visibleFaces[numVisibleFaces+change+1:numVisibleFaces]
			for fidx in start-cont:lastFace-cont
				faces[fidx] = faces[fidx+1]
				norms[fidx] = norms[fidx+1]
			end
			cont += 1;
		end
	end

	return lastFace + change
end

cross( a, b ) = a[2]*b[3] - a[3]*b[2], a[3]*b[1] - a[1]*b[3], a[1]*b[2] - a[2]*b[1]

function areSame( x, y )

	#println( "same test: ", (x, y) )
	return x == y
end

function areColinear( x, y, z )

	v1 = x .- z
	v2 = y .- z
	n  = cross( v1, v2 )

	#println( "colinear test: ", n )

	return n[1] == 0 && n[2] == 0 && n[3] == 0
end

function areCoplanar( x, y, z, p )

	v1 = ( x[1] - p[1], x[2] - p[2], x[3] - p[3] )
	v2 = ( y[1] - p[1], y[2] - p[2], y[3] - p[3] )
	v3 = ( z[1] - p[1], z[2] - p[2], z[3] - p[3] )

 	test = sum( v1 .* cross( v2, v3 ) )

	#println( "coplanar test: ", test )

    return test == 0
end

function fixSameColinearCoplanar!( points )

	idx = 2
	while areSame( points[1], points[idx] ) && idx < length(points)
		idx += 1;
	end
	if ( idx !== 2 )
		temp = points[idx]
		points[idx] = points[2]
		points[2] = temp
	end

	idx = 3
	while areColinear( points[1], points[2], points[idx] ) && idx < length(points)
		idx += 1;
	end
	if ( idx !== 3 )
		temp = points[idx]
		points[idx] = points[3]
		points[3] = temp
	end

	idx = 4
	while areCoplanar( points[1], points[2], points[3], points[idx] ) && idx < length(points)
		idx += 1;
	end
	if ( idx !== 4 )
		temp = points[idx]
		points[idx] = points[4]
		points[4] = temp
	end
end

function convexHullSIMP2( points::Array{Tuple{T,T,T},1}; typ=UInt16 ) where {T<:Real}

    @assert length(points) >= 4
    # length(points) for now, to have a static size

    faces = Array{Tuple{typ,typ,typ},1}(undef,0);
    norms = Array{Tuple{Float32,Float32,Float32},1}(undef,0);

    push!( faces, (1,2,3) )
	push!( faces, (2,3,4) )
	push!( faces, (3,4,1) )
	push!( faces, (4,1,2) )

	center   = ( points[1] .+ points[2] .+ points[3] .+ points[4] ) ./ 4;
    for fidx in 1:length(faces)
		push!( norms, faceNormalCenterSIMP( points, faces[fidx], center ) )
    end

	@inbounds for pidx in 5:length(points);

		new_point = points[pidx];

		visiblefaces = Array{typ,1}(undef,0)

		for fidx in 1:length(faces);

			vp = new_point .- points[faces[fidx][1]];
			if sign( sum( vp .* norms[ fidx ] ) ) == 1
				push!( visibleFaces, fidx );
			end
		end

		numVisibleFaces = length(visibleFaces);

		lastFace = addNewFacesSIMP2!( points, faces, norms, center, lastFace,
						             pidx, visibleFaces, numVisibleFaces );

	end # each input point

    return faces[1:lastFace], norms[1:lastFace]
end

"""
	needs to be tested
"""
# can this be made more compact?
@inline function addNewFacesSIMP2!( points,
                               faces::Array{Tuple{T,T,T},1},
							   norms::Array{Tuple{Float32,Float32,Float32},1},
							   center,
							   lastFace,
							   newpointindex,
							   visibleFaces::Array{T,1},
                               numVisibleFaces
                             ) where {T<:Integer}

	boundary   = Array{Tuple{T,T},1}(undef, numVisibleFaces*3);
	isRepeated = zeros( Bool, numVisibleFaces*3 );

	lidx, nrep = 1, 0;

	@inbounds for face in faces[ visibleFaces[1:numVisibleFaces] ]
		edges = [(face[1],face[2]), (face[2],face[3]), (face[3],face[1])];
		for (p1, p2) in edges
			boundary[lidx] = ( p1, p2 )
			for bidx in 1:lidx-1
				boundary_edge = boundary[bidx]
				if (p1,p2) == boundary_edge || (p2,p1) == boundary_edge
					isRepeated[lidx] = true;
					isRepeated[bidx] = true;
					nrep += 2;
				end
			end # each previous boundary --> check that the current edge has not already been added
		lidx += 1
		end # each of the 3 edges --> add it to "boundary" and check if it is not repeated
	end # each visible face --> go through its 3 edges

	"""
		2-. From each unique edge in the boundary -> build a new face with the correct neighbouring data
	"""
	numNewFaces = lidx - 1 - nrep;

	change = numNewFaces - numVisibleFaces;

	#( numNewFaces < numVisibleFaces ) && ( println( "Delete faces" ) );

	lidx = 0;
	for bidx in 1:length(boundary)
		if !(isRepeated[bidx])
			lidx += 1
			newfaceidx = ( lidx <= numVisibleFaces ) ? visibleFaces[lidx] : lastFace + lidx - numVisibleFaces;
			faces[newfaceidx] = ( boundary[bidx][1], boundary[bidx][2], newpointindex );
			norms[newfaceidx] = faceNormalCenterSIMP( points, ( boundary[bidx][1], boundary[bidx][2], newpointindex ), center )
		end
	end

	if ( change < 0 )
		cont = 0
		for start in visibleFaces[numVisibleFaces+change+1:numVisibleFaces]
			for fidx in start-cont:lastFace-cont
				faces[fidx] = faces[fidx+1]
				norms[fidx] = norms[fidx+1]
			end
			cont += 1;
		end
	end

	return lastFace + change
end

function areCoplanarSIMP( x, y, z, p )

    a1 = x[2] - x[1] ;
    b1 = y[2] - y[1] ;
    c1 = z[2] - z[1] ;
    a2 = x[3] - x[1] ;
    b2 = y[3] - y[1] ;
    c2 = z[3] - z[1] ;
    a = b1 * c2 - b2 * c1 ;
    b = a2 * c1 - a1 * c2 ;
    c = a1 * b2 - b1 * a2 ;
    d = - a * x[1] - b * y[1] - c * z[1];

    return a * p[1] + b * p[2] + c * p[3] + d == 0
end


text:
using Makie
using MultipleViewGeometry, Test, Random
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using MultipleViewGeometry.ModuleConstraints
using MultipleViewGeometry.ModuleConstruct
using MultipleViewGeometry.ModuleDraw
using MultipleViewGeometry.ModuleMove
using LinearAlgebra
using StaticArrays
using GeometryTypes

# Fix random seed.
Random.seed!(1234)
# Construct two camera matrices and parametrise two planar surfaces.
f = 50
image_width = 640 / 10
image_height = 480 / 10
𝐊₁ = @SMatrix [f 0 0 ;
               0 f 0 ;
               0 0 1 ]
𝐑₁ = SMatrix{3,3,Float64,9}(rotxyz(0, 25*(pi/180), 0))
𝐭₁ = [-30.0, 0.0, -5.0]

𝐊₂ = @SMatrix [f 0 0 ;
               0 f 0 ;
               0 0 1 ]

𝐑₂ = SMatrix{3,3,Float64,9}(rotxyz(0, -25*(pi/180), 0))
𝐭₂ = [30.0, 0.0, 5.0]

# Normals and distance from origin
𝐧₁ = [0.0, 0.0, 1.0]
d₁ = 55.0

# Normals and distance from origin
𝐧₁ = [0.0, 0.0, 1.0]
d₁ = 55.0

𝐧₂ = [0.5, -0.2, 2.0]
d₂ = 145.0

# Generates a random point on the plane centered around a point on the plane
# that is closest to the origin.
function random_points_on_plane(𝐧::AbstractArray, d::Real, extent::Real, N::Int)
    # Generate vector 𝐰 on a plane through the origin with normal vector 𝐧.
    first(𝐧) == 0 ? 𝐰 = cross(𝐧,SVector(1.0,0.0,0.0)) : 𝐰 = cross(𝐧,SVector(0.0,0.0,1.0))
    points = Array{SVector{3,Float64},1}(undef,N)
    for n = 1:N
        # Rotate 𝐰 randomly around the axis 𝐧.
        θ = rand() * 2*pi
        𝐤 = 𝐧 / norm(𝐧)
        𝐯 = 𝐰*cos(θ) + cross(𝐤,𝐰)*sin(θ) + 𝐤*dot(𝐤,𝐰)*(1-cos(θ))
        # Scale the vector so that it lies in the interval [0, extent)
        𝐯 = (rand() * extent) * 𝐯
        # Translate the vector so that it lies on the plane parametrised by 𝐧 and d.
        𝐯 = 𝐯 + d*(𝐧/norm(𝐧)^2)
        points[n] = 𝐯
    end
    points
end

𝒳₁ = random_points_on_plane(𝐧₁,d₁, 20, 50)
𝒳₂ = random_points_on_plane(𝐧₂,d₂, 20, 50)
X₁ = reshape(reinterpret(Float64,𝒳₁),(3,length(𝒳₁)))
X₂ = reshape(reinterpret(Float64,𝒳₂),(3,length(𝒳₂)))

scale = 20.0f0
x = Vec3f0(0); baselen = 0.2f0 * scale ; dirlen = 1f0 * scale
# create an array of differently colored boxes in the direction of the 3 axes
rectangles = [
    (HyperRectangle(Vec3f0(x), Vec3f0(dirlen, baselen, baselen)), RGBAf0(1,0,0,1)),
    (HyperRectangle(Vec3f0(x), Vec3f0(baselen, dirlen, baselen)), RGBAf0(0,1,0,1)),
    (HyperRectangle(Vec3f0(x), Vec3f0(baselen, baselen, dirlen)), RGBAf0(0,0,1,1))
]
meshes = map(GLNormalMesh, rectangles)

scene = mesh(merge(meshes))
scatter!(scene, X₁[1,:],X₁[2,:], X₁[3,:], markersize = 3, color = :red)
scatter!(scene, X₂[1,:],X₂[2,:], X₂[3,:], markersize = 3, color = :blue)

world_basis = (Vec(1.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(0.0, 0.0, 1.0))
camera_basis = (Point(0.0, 0.0, 0.0), Vec(-1.0, 0.0, 0.0), Vec(0.0, -1.0, 0.0), Vec(0.0, 0.0, 1.0))
picture_basis = (Point(0.0, 0.0), Vec(-1.0, 0.0), Vec(0.0, -1.0))

camera₁ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)
camera₂ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)
relocate!(camera₁, 𝐑₁, 𝐭₁)
relocate!(camera₂, 𝐑₂, 𝐭₂)

draw!(camera₁, scene)
draw!(camera₂, scene)

𝐑₁′, 𝐭₁′ = ascertain_pose(camera₁, world_basis... )
𝐊₁′ = obtain_intrinsics(camera₁, CartesianSystem())
𝐑₂′, 𝐭₂′ = ascertain_pose(camera₂, world_basis... )
𝐊₂′ = obtain_intrinsics(camera₂, CartesianSystem())

𝐏₁ = construct(ProjectionMatrix(),𝐊₁′,𝐑₁′,𝐭₁′)
𝐏₂ = construct(ProjectionMatrix(),𝐊₂′,𝐑₂′,𝐭₂′)

# Set of corresponding points.
ℳ₁ = project(camera₁,𝐏₁,𝒳₁)
ℳ₁ʹ= project(camera₂,𝐏₂,𝒳₁)
ℳ₂ = project(camera₁,𝐏₁,𝒳₂)
ℳ₂ʹ= project(camera₂,𝐏₂,𝒳₂)

M₁ = reshape(reinterpret(Float64,ℳ₁),(2,length(ℳ₁)))
M₂ = reshape(reinterpret(Float64,ℳ₂),(2,length(ℳ₂)))
M₁ʹ = reshape(reinterpret(Float64,ℳ₁ʹ),(2,length(ℳ₁ʹ)))
M₂ʹ = reshape(reinterpret(Float64,ℳ₂ʹ),(2,length(ℳ₂ʹ)))

# Visualize the set of corresponding points
scene = Scene()
scatter!(scene,M₁[1,:], M₁[2,:], markersize = 1, color = :red, limits = FRect(0, 0, 64, 48.0))
scatter!(scene,M₂[1,:], M₂[2,:], markersize = 1, color = :blue, limits = FRect(0, 0, 64, 48.0))

scene = Scene()
scatter!(scene,M₁ʹ[1,:], M₁ʹ[2,:], markersize = 1, color = :red, limits = FRect(0, 0, 64, 48.0))
scatter!(scene,M₂ʹ[1,:], M₂ʹ[2,:], markersize = 1, color = :blue, limits = FRect(0, 0, 64, 48.0))

# We will construct a pair of homography matrices and then construct a pair of
# projection matrices from the homographies.

𝐇₁ = construct(HomographyMatrix(),𝐊₁′,𝐑₁′,-𝐭₁′,𝐊₂′,𝐑₂′,-𝐭₂′,𝐧₁,d₁)
𝐇₂ = construct(HomographyMatrix(),𝐊₁′,𝐑₁′,-𝐭₁′,𝐊₂′,𝐑₂′,-𝐭₂′,𝐧₂,d₂)

context = ProjectionMatrices(HomographyMatrices(), Chojnacki(), TwoViews())
𝐎₁, 𝐎₂ = construct(context, (𝐇₁,𝐇₂))
𝐅₁ = construct(FundamentalMatrix(), 𝐏₁, 𝐏₂)
𝐅₂ = construct(FundamentalMatrix(), 𝐎₁, 𝐎₂)

𝐅₁ = 𝐅₁ / 𝐅₁[3,3]
𝐅₂ = 𝐅₂ / 𝐅₂[3,3]

# Use homographies to determine corresponding coordinates in the second image.
𝒪₁ = similar(ℳ₁ʹ)
𝒪₂ = similar(ℳ₂ʹ)
for n = 1:length(ℳ₁)
    𝐦 = hom(ℳ₁[n])
    𝐦ʹ = 𝐇₁*𝐦
    𝒪₁[n] = hom⁻¹(𝐦ʹ)
end
for n = 1:length(ℳ₂)
    𝐦 = hom(ℳ₂[n])
    𝐦ʹ = 𝐇₂*𝐦
    𝒪₂[n] = hom⁻¹(𝐦ʹ)
end

#TODO Test that 𝒪₁ .== ℳ₁ʹ and 𝒪₂ .== M₂ʹ


text:
using StringDistances, Random
Random.seed!(2)
x = map(Random.randstring, rand(5:25,500_000))
y = map(Random.randstring, rand(5:25,500_000))

function f(t, x, y; min_score = 0.0)
    [compare(x[i], y[i], t; min_score = min_score) for i in 1:length(x)]
end

@time f(Hamming(), x, y)
#0.05s
@time f(Jaro(), x, y)
#0.3s
@time f(Levenshtein(), x, y)
# 0.35s. A bit faster than StringDist
@time f(Levenshtein(), x, y, min_score = 0.8)
# 0.11
@time f(DamerauLevenshtein(), x, y)
# 0.45s.  Much faster than StringDist
@time f(DamerauLevenshtein(), x, y, min_score = 0.8)
# 0.08

@time find_best(x[1], y, Levenshtein())
# 0.41
@time find_best(x[1], y, DamerauLevenshtein())
# 0.41

@time find_all(x[1], y, Levenshtein())
# 0.14
@time find_all(x[1], y, DamerauLevenshtein())
# 0.07
@time find_all(x[1], y, Partial(DamerauLevenshtein()))
# 2.9

@time find_all(x[1], y, TokenSort(DamerauLevenshtein()))
# 0.7
@time find_all(x[1], y, TokenSet(DamerauLevenshtein()))
# 1.5
@time find_all(x[1], y, TokenMax(DamerauLevenshtein()))
# 5.4

# 1.6s slower compared to StringDist

# check
function h(t, x, y; min_score = 1/3)
	out = fill(false, length(x))
	for i in eachindex(x)
		if compare(x[i], y[i], t) <  min_score
			out[i] = compare(x[i], y[i], t ; min_score = min_score) ≈ 0.0
			else
			out[i] = compare(x[i], y[i], t ; min_score = min_score) ≈ compare(x[i], y[i], t)
		end
	end
	all(out)
end
h(Levenshtein(), x, y)
h(DamerauLevenshtein(), x, y)

#= Rcode
library(stringdist)
x <- sapply(sample(5:25,5 * 1e5,replace=TRUE), function(n) paste(sample(letters,n,replace=TRUE),collapse=""))
y <- sapply(sample(5:25,5 * 1e5,replace=TRUE), function(n) paste(sample(letters,n,replace=TRUE),collapse=""))
system.time(stringdist(x,y,method='lv', nthread = 1))
system.time(stringdist(x,y,method='dl', nthread = 1))

#  0.472
system.time(stringdist(x,y,method='jaccard', nthread = 1))
# 0.739
system.time(stringdist(x,y,method='cosine', nthread = 1))
system.time(stringdist(x,y,method='qgram', nthread = 1))

=#


text:
using TaskGraphs

## set up the environment
using GraphUtils

indicator_grid = IndicatorGrid(zeros(Bool,7,5))                           # initialize a 10x10 environment with no obstacles
indicator_grid[4,4] = true
indicator_grid[1,3] = true
indicator_grid[2,3] = true                                                             # obstacle at position 4,4
indicator_grid[1:6,2] .= true
indicator_grid[1:6,4] .= true                                                          # obstacle covering positions 6:7,6:7
env = construct_factory_env_from_indicator_grid(indicator_grid) # initialize the factory environment
@show env.vtx_map                                                                     # visualize the mapping from grid position to vertex indices

#   1   0   0   0   2
#   3   0   0   0   4
#   5   0   6   0   7
#   8   0   9   0  10
#  11   0  12   0  13
#  14   0  15   0  16
#  17  18  19  20  21

## set up the environment
## Define the initial conditions of the robots
robot_ics = [
    ROBOT_AT(2,12), # robot 1 starts at vertex 2
    ROBOT_AT(1,15), # robot 2 starts at vertex 9
]

## Define the manufacturing project
spec = ProjectSpec()
# set initial conditions of "raw materials"
set_initial_condition!(spec,OBJECT_AT(1,15)) # object 1 starts at vertex 4
set_initial_condition!(spec,OBJECT_AT(2,12))  # object 2 starts at vertex

# define the operations that need to take place

op1 = Operation(Δt=0)
op2 = Operation(Δt=0)

# inputs
set_precondition!(op1,OBJECT_AT(1,))
set_precondition!(op2,OBJECT_AT(2,12))

# outputs
# set_postcondition!(op1,OBJECT_AT(3,14))
# set_postcondition!(op2,OBJECT_AT(4,8))
add_operation!(spec,op1)
add_operation!(spec,op2)

# define assignment_dict
assignment_dict = Dict(
  1=>[1],2=>[2],
)
# export
#     apply_assignment_dict!

## define solver
solver = NBSSolver()
# finalize problem construction (the solver is passed as an argument here
# because it determines what cost model is used for the problem)
prob = pctapf_problem(solver,spec,env,robot_ics)
TaskGraphs.apply_assignment_dict!(get_schedule(get_env(prob)),assignment_dict,get_problem_spec(get_env(prob)))
prob = PC_MAPF(get_env(prob))
route_planner = solver.path_planner
# solve the problem
solution, cost = solve!(route_planner,prob)
println(env.vtx_map)
println(solution)
println(cost)
# check if the problem was solved to optimality
@show feasible_status(solver)
@show optimal_status(solver)


text:
"Simulation Controls"
global NUMBER_OF_WALKS = 1000
global MAX_STEPS_PER_WALK = 1485604
global PPD_ON = true  # no PPD == false
const FLOW_BIAS = true
const CATALASE_ON_WALLS = true

"Solvent Diffusion Step Sizes (microns per 2.7 microseconds)"
function set_PPD_STEP_SIZE(ppd_step::AbstractFloat=0.001)
    if PPD_ON
        global PPD_STEP_SIZE = ppd_step
    else
        global PPD_STEP_SIZE = WATER_STEP_SIZE
    end
    global STEP_SIZE_DICT = Dict("water" => WATER_STEP_SIZE,
                                "enz" => 0,
                                "ppd" => PPD_STEP_SIZE)
    nothing
end

function set_PPD_ON(trueorfalse::Bool=true)
    global PPD_ON = trueorfalse
    set_PPD_STEP_SIZE()
    nothing
end

const SECONDS_PER_STEP = 0.0000027077
const WATER_STEP_SIZE = 0.1
const ENZ_STEP_SIZE = 0
set_PPD_STEP_SIZE()

"Sensor/Walls Coordinates"
const SENSOR_HALF_WIDTH = 0.5 * 25
const WALL_Y = 1.5

"PPD Layer Dimensions"
const PPD_MAX_Y = 0.15
const PPD_MIN_Y = 0

"No Enzyme Layer (Variable Stubs)"
const ENZYME_LEFT_X = 0
const ENZYME_RIGHT_X = 0
const ENZYME_MAX_Y_FROM_WALL = 0
const ENZYME_MAX_Y = 0
const BORDER_CORRECTION = 0
const SPAWN_RIGHT_X = 0
const SPAWN_ENZYME_MAX_Y = 0
const CORNER_CUT_IN_ENZ = 0
const CORNER_CUT_IN_WATER = 0
const WATER_TO_ENZ_NORTHEAST = [0, 0]
const WATER_TO_ENZ_NORTHWEST = [0, 0]
const ENZ_TO_WATER_NORTHEAST = [0, 0]
const ENZ_TO_WATER_NORTHWEST = [0, 0]

"Rectangular Spawn Region from Water"
const MARGIN_OF_SPAWN = 5 * WATER_STEP_SIZE
const SPAWN_LEFT_X = MARGIN_OF_SPAWN - SENSOR_HALF_WIDTH
const SPAWN_TOP_LEFT_XY = [SPAWN_LEFT_X, WALL_Y - MARGIN_OF_SPAWN]
const SPAWN_X_RANGE = 2 * SENSOR_HALF_WIDTH - 2 * MARGIN_OF_SPAWN
const SPAWN_Y_RANGE = WALL_Y - PPD_MAX_Y - 2 * MARGIN_OF_SPAWN

"Escape Bound Limits"
const ESCAPE_X = 5 * SENSOR_HALF_WIDTH
const ESCAPE_Y = ESCAPE_X

"Safe Bound Limits"
const MARGIN_OF_COLLISION = 2 * WATER_STEP_SIZE
const SAFE_MIN_Y = WALL_Y + MARGIN_OF_COLLISION
const SAFE_MAX_Y = ESCAPE_Y - MARGIN_OF_COLLISION
const SAFE_MAX_X = ESCAPE_X - MARGIN_OF_COLLISION

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
function solver(Nt::Int , Δt::Float64, L::Float64, N::Int, ɸ::Array, x::Array{Float64,1}, ẋ::Array{Float64,1},
    parameter_ɸ, para_x, para_V, ɸ_fct::Function, x_fct::Function, V_fct::Function, times_save=reverse([i for i = 1:Nt]))
    V = ([real.(zero(ɸ)) for _ = 1:length(size(ɸ))+1])
    xarray = [i for i in range(-L / 2.0, L / 2.0, length = N)]
    update_V = ((V, x, t) -> V_fct(V, x, t, xarray, para_V...))
    ɸ_step = ((ɸ, V, t) -> ɸ_fct(ɸ, V, t, Δt, parameter_ɸ...))
    a_step = ((ɸ, V) -> x_fct(ɸ, V, para_x...))
    ɸ_r = [zero(ɸ) for _ = 1:length(times_save)]
    x_r = [zero(x) for _ = 1:length(times_save)]
    ẋ_r = [zero(ẋ) for _ = 1:length(times_save)]
    ẍ_r = [zero(ẋ) for _ = 1:length(times_save)]
    j = 1
    V = update_V(V, x, 0.0)
    temp = real.(a_step(ɸ, V))
    for i = 1:Nt
        a = temp
        x = (x + Δt* ẋ + 0.5 * Δt^2 * a )
        ɸ = ɸ_step(ɸ, V[1], i * Δt)
        V = update_V(V, x, i * Δt)
        temp = a_step(ɸ, V)
        ẋ = (ẋ + 0.5 * Δt * (a+temp))

        if i == times_save[length(times_save)] # O(1)
            ɸ_r[j] = deepcopy(ɸ)
            x_r[j] = deepcopy(x)
            ẋ_r[j] = deepcopy(ẋ)
            ẍ_r[j] = deepcopy(a)
            j += 1
            pop!(times_save) # O(1)
        end
    end
    return ɸ_r, x_r, ẋ_r, ẍ_r
end

max_stars_repo_name: phyjonas/ColloidFluctuatingField


text:
module ContentTypeDetection

using URIs
import ..sniff
import ..Form
using ..Messages
import ..IOExtras
import ..@debug, ..DEBUG_LEVEL

export contenttypedetectionlayer

function contenttypedetectionlayer(handler)
    return function(req; detect_content_type::Bool=false, kw...)
        if detect_content_type && (!hasheader(req.headers, "Content-Type")
            && !isa(req.body, Form)
            && isbytes(req.body))

            sn = sniff(bytes(req.body))
            setheader(req.headers, "Content-Type" => sn)
            @debug 1 "setting Content-Type header to: $sn"
        end
        return handler(req; kw...)
    end
end

end # module


text:
"""
    _rategradient(∂X, ∂Y)

Returns the rate of change in units of the values stored in the grid, and the
angle of the change in wind direction, *i.e.* an angle of 180 means that the
value is increasing *from* the south. When both ∂X and ∂Y are equal to 0, the
angle is assumed to be 0.
"""
function _rategradient(∂X::T, ∂Y::T) where {T <: Number}
    if ∂X == ∂Y == 0.0
        return (0.0, 0.0)
    end
    m = sqrt(∂X^2 + ∂Y^2)
    Δ = ∂X >= 0.0 ? 0.0 : 180.0
    θ = rad2deg(atan(∂X, ∂Y)) + Δ
    θ = ∂X > 0.0 ? θ + 180.0 : θ
    return (m, θ)
end

"""
    _rateofchange(x::Vector{T}, y::Vector{T}, z::Vector{T})

Rate of change for a series of three points, defined as a series of `x` and `y`
coordinates and a value `z`. Returns a rate of change (in units of `z`) and a
gradient in degrees.
"""
function _rateofchange(x::Vector{T}, y::Vector{T}, z::Vector{T}) where {T <: Number}

    # Check that all three vectors have the same length
    length(x) == length(y) || throw(DimensionMismatch("x and y must have the same length"))
    length(x) == length(z) || throw(DimensionMismatch("x and z must have the same length"))

    # Get the matrix of coefficients
    C = cat(y, x, fill(one(T), length(x)); dims=(2, 2))
    coeff = Base.inv(C) * z

    X = sum(C[:, 1]) / 3.0
    Y = sum(C[:, 2]) / 3.0

    ∂X = coeff[2] * Y + coeff[3]
    ∂Y = coeff[1] * X + coeff[3]

    # Rate of change and direction
    return _rategradient(∂X, ∂Y)
end

"""
    _rateofchange(A::Matrix{T}) where {T <: Number}

Rate of change for a series of a lattice points. Returns the rate of change and
the gradient for a 2x2 grid of numbers.
"""
function _rateofchange(x::Vector{T}, y::Vector{T}, z::Matrix{T};) where {T <: Number}
    size(z) == (2, 2) || throw(DimensionMismatch("the matrix z must have size (2,2)"))

    # The values in x and y should be sorted, but we'll still enforce it
    Δx = abs(first(diff(x)))
    Δy = abs(first(diff(y)))

    # We can get the values directly from the matrix
    Z₄, Z₁, Z₃, Z₂ = z

    ∂X = Z₂ - Z₁ + Δy * (Z₁ - Z₂ + Z₃ - Z₄)
    ∂Y = Z₄ - Z₁ + Δx * (Z₁ - Z₂ + Z₃ - Z₄)

    # Rate of change and direction
    return _rategradient(∂X, ∂Y)
end

_rateofchange(z::Matrix{T}) where {T <: Number} = _rateofchange(one(T) / 2, one(T) / 2, z)


text:
"""
   Implements basic ideas in

      S. Olver and A. Townsend,
      SIAM review, 2013,
      arXiv:1202.1347

   For more reviews, see Sec. 3 of

      A. Townsend and S. Olver,
      The automatic solution of partial differential
      equations using a global spectral method,
      Journal of Computational Physics, 2015
      arXiv:1409.2789

   and/or Sec. 6 of

      A. Townsend,
      Computing with functions in two dimensions,
      PhD Thesis, University of Oxford
      (can be accessed at https://pi.math.cornell.edu/~ajt/)

"""
module Gegenbauer

export compute_S, compute_D, compute_MX1, compute_MXp

include("CustomTypes.jl")
include("Chebyshev.jl")
using .CustomTypes
using .Chebyshev
using LinearAlgebra, SparseArrays

"""
   compute_S(n::myI, lam::myI)

   Computes n×n S_{λ} matrix.
"""
function compute_S(n::myI, lam::myI)
   @assert lam>=0
   @assert n>=3
   X = Vector{myI}(undef,0)
   Y = Vector{myI}(undef,0)
   V = Vector{myF}(undef,0)
   push!(X,1)
   push!(Y,1)
   push!(V,tomyF(1.0))
   if lam==0
      push!(X,1)
      push!(Y,3)
      push!(V,tomyF(-0.5))
      for i=2:(n-2)
         push!(X,i)
         push!(Y,i)
         push!(V,tomyF(0.5))
         push!(X,i)
         push!(Y,i+2)
         push!(V,tomyF(-0.5))
      end
      push!(X,n-1)
      push!(Y,n-1)
      push!(V,tomyF(0.5))
      push!(X,n)
      push!(Y,n)
      push!(V,tomyF(0.5))
   else
      push!(X,1)
      push!(Y,3)
      push!(V,tomyF(-lam/(lam+2.0)))
      for i=2:(n-2)
         push!(X,i)
         push!(Y,i)
         push!(V,tomyF(lam/(lam+i-1.0)))
         push!(X,i)
         push!(Y,i+2)
         push!(V,tomyF(-lam/(lam+i+1.0)))
      end
      push!(X,n-1)
      push!(Y,n-1)
      push!(V,tomyF(lam/(lam+n-2.0)))
      push!(X,n)
      push!(Y,n)
      push!(V,tomyF(lam/(lam+n-1.0)))
   end
   return sparse(X,Y,V)
end

"""
   compute_D(n::myI, lam::myI)

   Computes n×n D_{λ} matrix.
"""
function compute_D(n::myI, lam::myI)
   @assert lam>=1
   @assert n>=lam
   X = Vector{myI}(undef,0)
   Y = Vector{myI}(undef,0)
   V = Vector{myF}(undef,0)
   for i=1:(n-lam)
      push!(X,i)
      push!(Y,i+lam)
      push!(V,tomyF(lam+i-1))
   end
   for i=1:lam
      push!(X,n-lam+i)
      push!(Y,n)
      push!(V,tomyF(0))
   end
   return dropzeros(tomyF((2^(lam-1))*factorial(lam-1)).*sparse(X,Y,V))
end

"""
   compute_MX1(n::myI, lam::myI)

   Computes n×n M_{λ} matrix for polynomial x.
   For a simple derivation see Sec. 6.3.1 of
      Computing with functions in two dimensions,
      Townsend, Alex
      PhD Thesis, University of Oxford
   This can be accessed at
   https://pi.math.cornell.edu/~ajt/
"""
function compute_MX1(n::myI, lam::myI)
   @assert lam>=1
   @assert n>=lam
   X = Vector{myI}(undef,0)
   Y = Vector{myI}(undef,0)
   V = Vector{myF}(undef,0)
   push!(X,1)
   push!(Y,2)
   push!(V,tomyF(lam/(lam+1.0)))
   for i=2:(n-1)
      push!(X,i)
      push!(Y,i+1)
      push!(V,tomyF(((2.0*lam)+i-1)/(2.0*(lam+i))))
      push!(X,i)
      push!(Y,i-1)
      push!(V,tomyF((i-1)/(2.0*(lam+i-2))))
   end
   push!(X,n)
   push!(Y,n-1)
   push!(V,tomyF((n-1)/(2.0*(lam+n-2))))
   return sparse(X,Y,V)
end

"""
   compute_MCp(n::myI, lam::myI, p::myI)

   Computes n×n M_{λ} matrix for polynomial C^{λ}_p.
   For a simple derivation see Sec. 6.3.1 of
      Computing with functions in two dimensions,
      Townsend, Alex
      PhD Thesis, University of Oxford
   This can be accessed at
   https://pi.math.cornell.edu/~ajt/
"""
function compute_MCp(n::myI, lam::myI, p::myI)
   @assert p>=0

   MX0 = sparse(I,n,n)
   MX1 = compute_MX1(n,lam)

   if p==0
      return MX0
   elseif p==1
      return (2.0*lam).*MX1
   else
      MCp1 = deepcopy(sparse(I,n,n))
      MC   = deepcopy((2.0*lam).*MX1)
      MCm1 = deepcopy(MX0)

      for k=1:(p-1)
         MCp1 = (
            (2.0*(k+lam)/(k+1)).*(MX1*MC)
            -
            ((k-1+2*lam)/(k+1)).*MCm1
           )
         MCm1 = deepcopy(MC)
         MC   = deepcopy(MCp1)
      end
      return MCp1
   end
end

"""
   compute_MXp(n::myI, lam::myI, p::myI)

   Computes n×n M_{λ} matrix for polynomial x_p.
   For a simple derivation see Sec. 6.3.1 of
      Computing with functions in two dimensions,
      Townsend, Alex
      PhD Thesis, University of Oxford
   This can be accessed at
   https://pi.math.cornell.edu/~ajt/
"""
function compute_MXp(n::myI, lam::myI, p::myI)
   @assert p>=0

   MCp = compute_MCp(n,lam,p)

   if p==0
      return dropzeros(MCp)
   elseif p==1
      return dropzeros((0.5/lam).*MCp)
   elseif p==2
      MC0 = compute_MCp(n,lam,0)
      return dropzeros(((0.5/lam)/(1.0+lam)).*(MCp + lam.*MC0))
   elseif p==3
      return nothing
   elseif p==4
      return nothing
   end
end

"""
   compute_D(n::myI, lam::myI)

   Computes n×n D_{λ} matrix over the interval [xmin,xmax].
"""
function compute_D(n::myI, lam::myI, xmin::myF, xmax::myF)
   D = compute_D(n,lam)

   a = tomyF(0.5*(xmax-xmin))

   return dropzeros((a^(-lam)).*D)
end

"""
   compute_MXp(n::myI, lam::myI, p::myI, xmin::myF, xmax::myF)

   Computes n×n M_{λ} matrix for polynomial x_p over the interval [xmin,xmax].
"""
function compute_MXp(n::myI, lam::myI, p::myI, xmin::myF, xmax::myF)

   a = tomyF(0.5*(xmax-xmin))
   b = tomyF(0.5*(xmax+xmin))

   if p==0
      return compute_MXp(n,lam,0)
   elseif p==1
      return dropzeros(
               a*compute_MXp(n,lam,1)
               +
               b*compute_MXp(n,lam,0)
              )
   elseif p==2
      return dropzeros(
               (a^2)*compute_MXp(n,lam,2)
               +
               2*a*b*compute_MXp(n,lam,1)
               +
               (b^2)*compute_MXp(n,lam,0)
              )
   end
end

end


text:
module Force

using Reexport

export AbstractForce, potential, total_force, virtual_workrate, erased

include("AbstractForces.jl")

end


text:
#For testing using the National Instruments USB-6002 DAQ.
rig_key = "dummy-6002"
push!(RIGS, rig_key)
#The device doesn't support buffered digital I/O so we use an analog output
#for camera pulse commands and analog input to measure feedback

const d6002_mappings = Dict("AO0"=>"axial piezo",
                      "AO1"=>"camera1",
                      "AI0"=>"axial piezo monitor",
                      "AI1"=>"camera1 frame monitor",
                      "AI2"=>"analogin1",
                      "AI3"=>"analogin2",
                      "AI4"=>"analogin3",
                      "AI5"=>"analogin4",
                      "AI6"=>"analogin5",
                      "AI7"=>"analogin6")

DEFAULT_DAQCHANS_TO_NAMES[rig_key] = d6002_mappings
DEFAULT_NAMES_TO_DAQCHANS[rig_key] = dictmap(reverse, d6002_mappings)

const d6002_aochans = map(x->"AO$(x)", 0:1)
const d6002_aichans = map(x->"AI$(x)", [0:7...])
const d6002_dochans = String[]
const d6002_dichans = String[]
const d6002_pos_ctrl_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["axial piezo"])
const d6002_pos_mon_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["axial piezo monitor"])
const d6002_cam_ctrl_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["camera1"])
const d6002_cam_mon_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["camera1 frame monitor"])
const d6002_laschans = String[]
const d6002_stimchans = String[]
const d6002_fixed_names = ["axial piezo", "axial piezo monitor", "camera1", "camera1 frame monitor"]
const d6002_fixed_daqchans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], d6002_fixed_names)

AO_CHANS[rig_key] = OrderedSet(d6002_aochans)
AI_CHANS[rig_key] = OrderedSet(d6002_aichans)
DO_CHANS[rig_key] = OrderedSet(d6002_dochans)
DI_CHANS[rig_key] = OrderedSet(d6002_dichans)
POS_CONTROL_CHANS[rig_key] = d6002_pos_ctrl_chans
POS_MONITOR_CHANS[rig_key] = d6002_pos_mon_chans
CAM_CONTROL_CHANS[rig_key] = d6002_cam_ctrl_chans
CAM_MONITOR_CHANS[rig_key] = d6002_cam_mon_chans
LAS_CONTROL_CHANS[rig_key] = d6002_laschans
STIM_CHANS[rig_key] =OrderedSet(d6002_stimchans)
FIXED_NAMES[rig_key] = d6002_fixed_names
FIXED_DAQ_CHANS[rig_key] = d6002_fixed_daqchans

RIG_CHIP_SIZES[rig_key] = PCO_EDGE_4_2_CHIP_SIZE
RIG_FRAMERATE_FUNCS[rig_key] = PCO_EDGE_4_2_FRAMERATE_FUNC

LASER_ON_TIME[rig_key] = 100 * Unitful.μs
LASER_OFF_TIME[rig_key] = 100 * Unitful.μs

CAMERA_ON_TIME[rig_key] = 19526.0 * Unitful.ns #This is the worst jitter measured with Edge 5.5 and 4.2 cameras, see ImagineInterface issue #18
CAMERA_OFF_TIME[rig_key] = 19526.0 * Unitful.ns #This is the worst jitter measured with Edge 5.5 and 4.2 cameras, see ImagineInterface issue #18

PIEZO_RANGES[rig_key] = (0.0μm .. 400.0μm, 0.0V .. 10.0V)
PIEZO_MAX_SPEED[rig_key] = 2000*Unitful.μm / Unitful.s
AO_RANGE[rig_key] = -10.0V .. 10.0V
AI_RANGE[rig_key] = AO_RANGE[rig_key] #TODO: make sure this is true.  (true if we are recording -10..10V on analog inputs)

max_forks_repo_name: HolyLab/ImagineInterface


text:
using Plexifications
using Base.Test: @testset

@testset "Cmplex" begin include("cmplex_test.jl") end
@testset "Prplex" begin include("prplex_test.jl") end
@testset "Nlplex" begin include("nlplex_test.jl") end

@testset "BiComplex" begin include("bicomplex_test.jl") end
@testset "BiPerplex" begin include("biperplex_test.jl") end
@testset "BiNilplex" begin include("binilplex_test.jl") end

@testset "DualComplex" begin include("dualcomplex_test.jl") end
@testset "DualPerplex" begin include("dualperplex_test.jl") end

@testset "TriComplex" begin include("tricomplex_test.jl") end
@testset "TriPerplex" begin include("triperplex_test.jl") end
@testset "TriNilplex" begin include("trinilplex_test.jl") end

@testset "HyperComplex" begin include("hypercomplex_test.jl") end
@testset "HyperPerplex" begin include("hyperperplex_test.jl") end


text:
function convert_arguments(T::Type{<:Scatter}, p::Point)
    return convert_arguments(T, coordinates(p))
end

function convert_arguments(T::Type{<:Scatter}, P::Vector{<:Point})
    return convert_arguments(T, [coordinates(p) for p in P])
end

function convert_arguments(T::Type{<:LineSegments}, p::Point)
    return convert_arguments(T, coordinates(p))
end

function convert_arguments(T::Type{<:LineSegments}, P::Vector{<:Point})
    return convert_arguments(T, [coordinates(p) for p in P])
end

# Makie seems to have trouble with SVectors, so just collect them into Vectors
function convert_arguments(T::Type{<:Scatter}, V::SVector{N, <:Point}) where {N}
    return convert_arguments(T, collect(V))
end

function convert_arguments(T::Type{<:LineSegments}, V::SVector{N, <:Point}) where {N}
    return convert_arguments(T, collect(V))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# This script actually solves the WachterDisasterRisk model with a risk-adjusted linearization
# and times the methods, if desired
using BenchmarkTools, RiskAdjustedLinearizations, Test
include("wachter.jl")

# Settings: what do you want to do?
time_methods        = false
numerical_algorithm = :relaxation
autodiff            = false

# Set up
autodiff_method = autodiff ? :forward : :central
m_wachter = WachterDisasterRisk()
m = inplace_wachter_disaster_risk(m_wachter)
z0 = copy(m.z)
y0 = copy(m.y)
Ψ0 = copy(m.Ψ)

# Solve!
solve!(m; algorithm = numerical_algorithm, autodiff = autodiff_method)

if numerical_algorithm == :relaxation
    sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "iterative_sss_output.jld2"), "r")

elseif numerical_algorithm == :homotopy
    sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "homotopy_sss_output.jld2"), "r")
end

@test isapprox(sssout["z"], m.z, atol=1e-4)
@test isapprox(sssout["y"], m.y, atol=1e-4)
@test isapprox(sssout["Psi"], m.Ψ, atol=1e-4)

if time_methods
    println("Deterministic steady state")
    @btime begin
        solve!(m, z0, y0; algorithm = :deterministic, autodiff = autodiff_method, verbose = :none)
    end

    # Use deterministic steady state as guesses
    solve!(m, z0, y0; algorithm = :deterministic, autodiff = autodiff_method, verbose = :none)
    zdet = copy(m.z)
    ydet = copy(m.y)
    Ψdet = copy(m.Ψ)

    println("Relaxation method")
    @btime begin # called the "iterative" method in the original paper
        solve!(m, zdet, ydet, Ψdet; algorithm = :relaxation, autodiff = autodiff_method, verbose = :none)
    end

    println("Relaxation method with Anderson acceleration")
    @btime begin # called the "iterative" method in the original paper
        solve!(m, zdet, ydet, Ψdet; algorithm = :relaxation, use_anderson = true, m = 3, autodiff = autodiff_method, verbose = :none)
    end

    println("Homotopy method")
    @btime begin # called the "continuation" method in the original paper, but is called homotopy in the original code
        solve!(m, zdet, ydet, Ψdet; algorithm = :homotopy, autodiff = autodiff_method, verbose = :none)
    end
end


text:
### Statistics

function countmap(x::RLEVector{T1,T2}) where {T1,T2}
  tally = Dict{T1,T2}()
  for (v,w) in zip(values(x), widths(x))
    tally[v] = get(tally,v,0) + w
  end
  tally
end

function mode(x::RLEVector)
  tally = countmap(x)
  which_max = argmax(collect(values(tally)))
  max_val = collect(keys(tally))[ which_max ]
  max_val
end


text:
@testset "block test" begin
    @testset "flatten test" begin
		circuit = QCircuit(5)
		block2 = QGateBlock()
		block1 = QGateBlock()
		add!(block1, X(1), X(2), X(3), Y(3))
		add!(block2, block1(5,3,2))
		add!(circuit, block2(5,4,3,2,1))
		flatCircuit = QCircuit(5)
		add!(flatCircuit, X(1), X(3), X(4), Y(4))
		@test flatten(circuit) == flatCircuit
	end

	@testset "block extract param test" begin
		circuit = QCircuit(5)
		block2 = QGateBlock()
		block1 = QGateBlock()
		add!(block1, Rx[π](1), Rx[π/2](2), Rx[π/4](3), Ry[π/8](3))
		add!(block2, block1(5,3,2))
		add!(block2, X(1))
		add!(block2, block1(2,3,5))
		params = extractParams(block2)
		@test params ≈ repeat([π, π/2,π/4, π/8],2)
	end

	@testset "block insert param test" begin
	    circuit = QCircuit(5)
		block2 = QGateBlock()
		block1 = QGateBlock()
		add!(block1, Rx(1), X(2), Rx(3), Ry(3))
		addCopy!(block2, block1(5,3,2))
		add!(block2, X(1))
		addCopy!(block2, block1(2,3,5))
		newparams = rand(Float64, 6)
		insertParams!(block2, newparams)
		@test extractParams(block2) ≈ newparams
	end

	@testset "const gate inverse test" begin
	    @test data(inverse(T)) ≈ data(TDAG)
	    @test data(inverse(S)) ≈ data(SDAG)
	end

	@testset "block/circuit inverse test" begin
	    state = MPSState(5)
	    ref = copy(state)
	    circuit = randomQCircuit(5,10)
	    circuit_inv = inverse(circuit)
	    apply!(state, circuit)
	    apply!(state, circuit_inv)
	    @test ref ≈ state
	end
end


text:
module PrintTypesTersely

import Base

# const _terseprint = Ref(true)
const _terseprint = Ref(false)

on() = set(true)
off() = set(false)

function set(a)
    _terseprint[] = a
end

function with_state(f::Function, a)
    orig_val = _terseprint[]
    _terseprint[] = a
    f()
    _terseprint[] = orig_val
end

function base_show_terse(io::IO, @nospecialize(x::Type))
    # print(io, typeof(x))
    if x isa Union
        print(io, "Union")
        Base.show_delim_array(io, Base.uniontypes(x), '{', ',', '}', false)
        return
    end
    x = Base.unwrap_unionall(x)
    print(io, "$(x.name.name){…}")
    return
end

@static if VERSION < v"1.6.0-"
function base_show_full_1_5(io::IO, @nospecialize(x::Type))
    # basically copied from the Julia sourcecode, seems it's one of most robust fixes to Pevňákoviny
    # specifically function show(io::IO, @nospecialize(x::Type))
    if x isa DataType
        Base.show_datatype(io, x)
        return
    elseif x isa Union
        if x.a isa DataType && Core.Compiler.typename(x.a) === Core.Compiler.typename(DenseArray)
            T, N = x.a.parameters
            if x == StridedArray{T,N}
                print(io, "StridedArray")
                Base.show_delim_array(io, (T,N), '{', ',', '}', false)
                return
            elseif x == StridedVecOrMat{T}
                print(io, "StridedVecOrMat")
                Base.show_delim_array(io, (T,), '{', ',', '}', false)
                return
            elseif StridedArray{T,N} <: x
                print(io, "Union")
                Base.show_delim_array(io, vcat(StridedArray{T,N}, Base.uniontypes(Core.Compiler.typesubtract(x, StridedArray{T,N}))), '{', ',', '}', false)
                return
            end
        end
        print(io, "Union")
        Base.show_delim_array(io, Base.uniontypes(x), '{', ',', '}', false)
        return
    end

    x::UnionAll
    if Base.print_without_params(x)
        return show(io, Base.unwrap_unionall(x).name)
    end

    if x.var.name === :_ || Base.io_has_tvar_name(io, x.var.name, x)
        counter = 1
        while true
            newname = Symbol(x.var.name, counter)
            if !Base.io_has_tvar_name(io, newname, x)
                newtv = TypeVar(newname, x.var.lb, x.var.ub)
                x = UnionAll(newtv, x{newtv})
                break
            end
            counter += 1
        end
    end

    show(IOContext(io, :unionall_env => x.var), x.body)
    print(io, " where ")
    show(io, x.var)
end
end

@static if VERSION >= v"1.6.0-"
function base_show_full_1_6(io::IO, @nospecialize(x::Type))
    # basically copied from the Julia sourcecode, seems it's one of most robust fixes to Pevňákoviny
    # specifically function show(io::IO, @nospecialize(x::Type))
    if Base.print_without_params(x)
        Base.show_type_name(io, Base.unwrap_unionall(x).name)
        return
    elseif get(io, :compact, true) && Base.show_typealias(io, x)
        return
    elseif x isa DataType
        Base.show_datatype(io, x)
        return
    elseif x isa Union
        if get(io, :compact, true)
            Base.show_unionaliases(io, x)
        else
            print(io, "Union")
            Base.show_delim_array(io, Base.uniontypes(x), '{', ',', '}', false)
        end
        return
    end

    x = x::UnionAll
    wheres = TypeVar[]
    let io = IOContext(io)
        while x isa UnionAll
            var = x.var
            if var.name === :_ || Base.io_has_tvar_name(io, var.name, x)
                counter = 1
                while true
                    newname = Symbol(var.name, counter)
                    if !Base.io_has_tvar_name(io, newname, x)
                        var = TypeVar(newname, var.lb, var.ub)
                        x = x{var}
                        break
                    end
                    counter += 1
                end
            else
                x = x.body
            end
            push!(wheres, var)
            io = IOContext(io, :unionall_env => var)
        end
        show(io, x)
    end
    Base.show_wheres(io, wheres)
end
end

function Base.show(io::IO, @nospecialize(x::Type))
    # println("_terseprint[] in Base.show: $(_terseprint[])")
    if _terseprint[]
        return base_show_terse(io, x)
    else
        # VERSION < v"1.6.0-" returns true for 1.5.* and older, but for rc and beta of 1.6 it returns true
        @static if VERSION < v"1.6.0-"
            return base_show_full_1_5(io, x)
        else
            return base_show_full_1_6(io, x)
        end
    end
end

end


text:
# Wilcoxon.jl
# Wilcoxon rank sum (Mann-Whitney U) and signed rank tests in Julia
#
# Copyright (C) 2012   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export SignedRankTest, ExactSignedRankTest, ApproximateSignedRankTest

## COMMON SIGNED RANK

# Automatic exact/normal selection
function SignedRankTest{T<:Real}(x::AbstractVector{T})
    (W, ranks, signs, tie_adjustment, n, median) = signedrankstats(x)
    n_nonzero = length(ranks)
    if n_nonzero <= 15 || (n_nonzero <= 50 && tie_adjustment == 0)
        ExactSignedRankTest(x, W, ranks, signs, tie_adjustment, n, median)
    else
        ApproximateSignedRankTest(x, W, ranks, signs, tie_adjustment, n, median)
    end
end
SignedRankTest{T<:Real,S<:Real}(x::AbstractVector{T}, y::AbstractVector{S}) = SignedRankTest(x - y)

# Get W and absolute ranks for signed rank test
function signedrankstats{S<:Real}(x::AbstractVector{S})
   nonzero_x = x[x .!= 0]
   (ranks, tieadj) = tiedrank_adj(abs(nonzero_x))
   W = 0.0
   for i = 1:length(nonzero_x)
       if nonzero_x[i] > 0
           W += ranks[i]
       end
   end
   (W, ranks, nonzero_x .> 0, tieadj, length(x), median(x))
end

## EXACT WILCOXON SIGNED RANK TEST

immutable ExactSignedRankTest{T<:Real} <: HypothesisTest
    vals::Vector{T} # original values
    W::Float64              # test statistic: Wilcoxon rank-sum statistic
    ranks::Vector{Float64}           # ranks without ties (zero values)
    signs::BitArray{1}      # signs of input of ranks
    tie_adjustment::Float64 # adjustment for ties
    n::Int                  # number of observations
    median::Float64         # sample median
end
ExactSignedRankTest{T<:Real}(x::AbstractVector{T}) =
    ExactSignedRankTest(x, signedrankstats(x)...)
ExactSignedRankTest{S<:Real,T<:Real}(x::AbstractVector{S}, y::AbstractVector{T}) =
    ExactSignedRankTest(x - y)

testname(::ExactSignedRankTest) = "Exact Wilcoxon signed rank test"
population_param_of_interest(x::ExactSignedRankTest) = ("Location parameter (pseudomedian)", 0, x.median) # parameter of interest: name, value under h0, point estimate

function show_params(io::IO, x::ExactSignedRankTest, ident)
    println(io, ident, "number of observations:      ", x.n)
    println(io, ident, "Wilcoxon rank-sum statistic: ", x.W)
    println(io, ident, "rank sums:                   ", [sum(x.ranks[x.signs]), sum(x.ranks[!x.signs])])
    println(io, ident, "adjustment for ties:         ", x.tie_adjustment)
end

# Enumerate all possible Wilcoxon rank-sum results for a given vector, determining left-
# and right-tailed p values
function signedrankenumerate(x::ExactSignedRankTest)
    le = 0
    gr = 0
    n = length(x.ranks)
    tot = 2^n
    for i = 0:tot-1
        # Interpret bits of i as signs to generate wp for all possible sign combinations
        Wp = 0
        b = i
        j = 1
        while b != 0
            Wp += (b & 1)*x.ranks[j]
            j += 1
            b >>= 1
        end
        le += Wp <= x.W
        gr += Wp >= x.W
    end
    (le/tot, gr/tot)
end

function pvalue(x::ExactSignedRankTest; tail=:both)
    n = length(x.ranks)
    if n == 0
        1
    elseif x.tie_adjustment == 0
        # Compute exact p-value using method from Rmath, which is fast but cannot account for ties
        if tail == :both
            if x.W <= n * (n + 1)/4
                p = 2 * psignrank(x.W, n, true)
            else
                p = 2 * psignrank(x.W - 1, n, false)
            end
        elseif tail == :left
            psignrank(x.W, n, true)
        elseif tail == :right
            psignrank(x.W - 1, n, false)
        end
    else
        # Compute exact p-value by enumerating all possible ranks in the tied data
        if tail == :both
            min(1, 2 * minimum(signedrankenumerate(x)))
        elseif tail == :left
            singedrankenumerate(x.W, x.ranks)[1]
        elseif tail == :right
            signedrankenumerate(x.W, x.ranks)[2]
        end
    end
end

ci(x::ExactSignedRankTest, alpha::Real=0.05; tail=:both) = calculate_ci(x.vals, alpha; tail=tail)

## APPROXIMATE SIGNED RANK TEST

immutable ApproximateSignedRankTest{T<:Real} <: HypothesisTest
    vals::Vector{T} # original values
    W::Float64              # test statistic: Wilcoxon rank-sum statistic
    ranks::Vector{Float64} # ranks without ties (zero values)
    signs::BitArray{1}      # signs of input of ranks
    tie_adjustment::Float64 # adjustment for ties
    n::Int                  # number of observations
    median::Float64         # sample median
    mu::Float64             # normal approximation: mean
    sigma::Float64          # normal approximation: std
end

function ApproximateSignedRankTest{T<:Real}(x::Vector, W::Float64, ranks::Vector{T}, signs::BitArray{1}, tie_adjustment::Float64, n::Int, median::Float64)
    nz = length(ranks) # num non-zeros
    mu = W - nz * (nz + 1)/4
    std = sqrt(nz * (nz + 1) * (2 * nz + 1) / 24 - tie_adjustment / 48)
    ApproximateSignedRankTest(x, W, ranks, signs, tie_adjustment, n, median, mu, std)
end
ApproximateSignedRankTest{T<:Real}(x::AbstractVector{T}) =
    ApproximateSignedRankTest(x, signedrankstats(x)...)
ApproximateSignedRankTest{S<:Real,T<:Real}(x::AbstractVector{S}, y::AbstractVector{T}) =
    ApproximateSignedRankTest(x - y)

testname(::ApproximateSignedRankTest) = "Approximate Wilcoxon signed rank test"
population_param_of_interest(x::ApproximateSignedRankTest) = ("Location parameter (pseudomedian)", 0, x.median) # parameter of interest: name, value under h0, point estimate

function show_params(io::IO, x::ApproximateSignedRankTest, ident)
    println(io, ident, "number of observations:      ", x.n)
    println(io, ident, "Wilcoxon rank-sum statistic: ", x.W)
    println(io, ident, "rank sums:                   ", [sum(x.ranks[x.signs]), sum(x.ranks[!x.signs])])
    println(io, ident, "adjustment for ties:         ", x.tie_adjustment)
    println(io, ident, "normal approximation (μ, σ): ", (x.mu, x.sigma))
end

function pvalue(x::ApproximateSignedRankTest; tail=:both)
    if x.mu == x.sigma == 0
        1
    else
        if tail == :both
            2 * ccdf(Normal(), abs(x.mu - 0.5 * sign(x.mu))/x.sigma)
        elseif tail == :left
            cdf(Normal(), (x.mu + 0.5)/x.sigma)
        elseif tail == :right
            ccdf(Normal(), (x.mu - 0.5)/x.sigma)
        end
    end
end

ci(x::ApproximateSignedRankTest, alpha::Real=0.05; tail=:both) = calculate_ci(x.vals, alpha; tail=tail)

# implementation method inspired by these notes: http://www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf
function calculate_ci(x::AbstractVector, alpha::Real=0.05; tail=:both)
    check_alpha(alpha)

    if tail == :both
        c = 1 - alpha
    else
        c = 1 - 2 * alpha
    end
    n = length(x)
    m = div(n * (n + 1), 2)
    k_range = 1:div(m, 2)
    l = [1 - 2 * psignrank(i, n, true) for i in k_range]
    k = indmin(abs(l-c))
    vals = Float64[]
    enumerated = enumerate(x)
    for (outer_index, outer_value) in enumerated
        for (inner_index, inner_value) in enumerated
            if outer_index > inner_index
                continue
            end
            push!(vals, (inner_value + outer_value) / 2)
        end
    end
    sort!(vals)
    left = vals[k + 1]
    right = vals[m - k]
    if tail == :both
        return (left, right)
    elseif tail == :left
        return (left, Inf)
    elseif tail == :right
        return (-Inf, right)
    end
end


text:
# We realise early stopping via exceptions

struct StopException{T} <: Exception
  ret::T
end
Stop = StopException(nothing)
stop() = throw(Stop)
stop(ret) = throw(StopException(ret))

"""
  contextmanager handling custom breakpoints with `stop()`

This is usually only used within creating a new continuable from a previous one

# Examples
```julia
@cont stop_at4(continuable) = stoppable(continuable) do x
  x == 4 && stop()
  cont(x)
end
```
"""
function stoppable(func, continuable, default_return = nothing)
  try
    continuable(func)
    default_return  # default returnvalue to be able to handle `stop(returnvalue)` savely
  catch exc
    if !isa(exc, StopException)
      rethrow(exc)
    end
    exc.ret
  end
end

"""
    Continuables.@ifsomething expr
If `expr` evaluates to `nothing`, equivalent to `return nothing`, otherwise the macro
evaluates to the value of `expr`. Not exported, useful for implementing iterators.

# Example

```jldoctest
julia> using Continuables

julia> Continuables.@ifsomething iterate(1:2)
(1, 1)
julia> let elt, state = Continuables.@ifsomething iterate(1:2, 2); println("not reached"); end
```
"""
macro ifsomething(ex)
    quote
        result = $(esc(ex))
        result === nothing && return nothing
        result
    end
end


text:
include("diagnostics.jl")
using DataFrames, Plots

export calc_auc
export roc_curve

"calc the AUC of a function given xs and ys"
function calc_auc(x::Array, y::Array)
    df = rename!(hcat(x, y) |> DataFrame, Symbol.(["x", "y"]))
    df = sort(df, [:x, :y])
    auc = 0
    for i in 2:nrow(df)
        # rectangular component
        x = df[i, :x][1] - df[i-1, :x][1]
        y = df[i-1, :y][1]
        auc += x*y
        # triangular component
        height = df[i, :y][1] - y
        auc += 0.5*x * height
    end
    return auc
end

"generates data necssary for ROC curve, array of true positive and false positive rates + AUC"
function roc_data(actual::Array{Int64}, pred::Array{Float64}, step = 0.01)
    x = 0:step:1
    tps = [true_positives(actual, predict_from_threshold(pred, threshold), rate=true) for threshold in x]
    fps = [false_positives(actual, predict_from_threshold(pred, threshold), rate=true) for threshold in x]
    auc = calc_auc(fps, tps)
    return tps, fps, auc
end

"""
_given array of predictions (0 ≤ x ≤ 1] and array of actuals, return ROC AUC plot_

#### parameters:
    actual : Array{Int}
        actual labels, ∈ [0,1]
    pred : Array{Float64}
        an array of prediction values 0 ≤ x ≤ 1

#### returns: Plots.Plot{Plots.GRBackend}
    plot of ROC AUC
"""
function roc_curve(actual::Array{Int}, pred::Array{Float64})
    # diagonal line
    p1 = plot(
        [0, 1],
        [0, 1],
        label = "",
        line = :dash,
        color = :black,
        ylim = (0, 1),
        xlim = (0, 1),
    )
    y, x, auc = roc_data(actual, pred)
    plot!(x, y, label = "AUC=$(round(auc, digits=4))", color = :blue)

    plot!(legend = :bottomright)
    plot!(
        ylabel = "True Positive Rate",
        xlabel = "False Positive Rate",
        title = "ROC and AUC",
    )
    return p1
end


text:
include("meta.jl")

mutable struct TaskSpec
    input::Dict{String, Any}
    output::Dict{String, Any}
    error::Dict{String, Any}
end

mutable struct TaskStatus
    status::Dict{String, Any}
end

mutable struct Task
    typemeta::TypeMeta
    objectmeta::ObjectMeta
    spec::TaskSpec
    status::TaskStatus
end

function convert(object::Dict{String, Any}, task::Task)

end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
# Julia script to take IR images using the FLIR VUE PRO R and a Raspberry Pi.
# Author: R. Vezy
# Date: 07/03/2021

# First, don't forget to activate and instantiate the project:
# using Pkg; Pkg.activate("."); Pkg.instantiate()

# !NB: please use `sudo pigpiod` in the terminal of your Raspberry Pi before using this
# ! script and after install of the PiGPIO package. See https://github.com/JuliaBerry/PiGPIO.jl

# Import the library:
using PiGPIO

# Select the pin you use for the signal (must be PWM-compatible):
pin = 18

# Instantiate the connexion:
p = Pi()

# Set the mode of the Pin to output
set_mode(p, pin, PiGPIO.OUTPUT)

# Set the PWM frequency at 100 Hz:
PiGPIO.set_PWM_frequency(p, pin, 50)
# NB: the FLIR Vue Pro uses LOW = 1ms/20ms,
# MID = 1.5ms/20ms, HIGH = 2ms/20ms at 3.3-5vdc, 50Hz.
# So we need a 0.002 s pulse. At 50Hz, we need a duty cycle
# of 10%.

# Set the range to 0-100 (easier):
PiGPIO.set_PWM_range(p, pin, 100)

waiting_time = 60 # (seconds) time between two pictures
time_for_transfer = 0.5

# This is the code doing the work:
while true
    try
        PiGPIO.set_PWM_dutycycle(p, pin, 10)
        sleep(time_for_transfer) # wait for the signal to transfer
        PiGPIO.set_PWM_dutycycle(p, pin, 1) # Reset the signal to LOW
        sleep(waiting_time - time_for_transfer) # wait `waiting_time` seconds
    catch e
        println("Error in execution: $e")
    end
end
# The while loop above will execute until the end of time if you let it continue. It also
# catches any issues such as disconnection or anything else, it will then proceed to usual
# whenever the connexion is right again.

# In case you need to stop the script, don't forget to put the pin to INPUT mode again:
set_mode(p, pin, PiGPIO.INPUT) # changed it back to INPUT mode


text:
#=
Author : Ng Yeesian
Modified : Joshua Ong /Guy Farmer / Michael Hilborn
generates the robust deployment model
=#
using AmbulanceDeployment

struct robustGamma
    _single::Vector{Int}
    _local::Vector{Int}
    _regional::Vector{Int}
    _global::Int
end
params = Params(0.01, 0.5, 1e-6, 500, 50)
struct Qrobust
    m::JuMP.Model
    I::UnitRange{Int}
    J::UnitRange{Int}
    d::Vector{JuMP.VariableRef}
    q::Vector{JuMP.VariableRef}
    γ::robustGamma
end

struct RobustDeployment <: DeploymentModel
    m::JuMP.Model
    Q::Qrobust

    I::UnitRange{Int}
    J::UnitRange{Int}

    x::Vector{JuMP.VariableRef}
    y::Vector{Matrix{JuMP.VariableRef}}
    z::Vector{Vector{JuMP.VariableRef}}
    η::JuMP.VariableRef

    scenarios::Vector{Vector{Int}}
    upperbounds::Vector{Float64}
    lowerbounds::Vector{Float64}
    deployment::Vector{Vector{Int}}

    upptiming::Vector{Float64}
    lowtiming::Vector{Float64}
end
deployment(m::RobustDeployment) = m.deployment[end]

function Gamma(p::DeploymentProblem; α=params.α)
    demand = p.demand[p.train,:]
    γ_single = Any[] # vec(maximum(demand,1) + 1 *(maximum(demand,1) .== 0))

    for x in 1:size(demand)[2]
        y = maximum(demand[:,x])
        if y == 0
            y = 1
        end
        push!(γ_single,y)
    end
    total_hours = 300 #size(demand)[1]
    γ_local = [quantile(Poisson(mean(sum(demand[:,vec(p.adjacency[i,:])]) / total_hours )),1-α) for i=1:p.nregions]
    γ_regional = [quantile(Poisson(mean(sum(demand[:,p.coverage[:,i]])/ total_hours)),1-α) for i in 1:p.nlocations]
    γ_global = quantile(Poisson(mean(sum(demand) / total_hours )),1-α) #way to many calls possible so I added to dicide by demand
    print("gamma global " * string(γ_global) )
    robustGamma(γ_single,γ_local,γ_regional,γ_global)
end

function Qrobust(problem::DeploymentProblem; α=params.α, verbose=false, solver=Gurobi.Optimizer(OutputFlag=0)) #, MIPGapAbs=0.9)
    if verbose
        solver=Gurobi.Optimizer(OutputFlag=1) #, MIPGapAbs=0.9)
    end
    γ = Gamma(problem, α=α)
    print("gamma global " * string(γ) )
    upp_bound = maximum(γ._single)
    I = 1:problem.nlocations
    J = 1:problem.nregions

    m = Model(GLPK.Optimizer)
    JuMP.@variable(m, d[1:problem.nregions]>=0, Int)
    JuMP.@variable(m, p[1:problem.nregions], Bin)
    JuMP.@variable(m, q[1:problem.nlocations], Bin)

    for i in I, j in J
        problem.coverage[j, i] && JuMP.@constraint(m, p[j] <= q[i])
    end

    # Uncertainty
    for j in J
        JuMP.@constraint(m, d[j] <= γ._single[j]*p[j])
        adjacent_regions = filter(k->problem.adjacency[k,j],J)
        JuMP.@constraint(m, sum(d[k] for k in adjacent_regions) <= γ._local[j])
    end
    for i in I
        covered_regions = filter(j->problem.coverage[j,i],J)
        JuMP.@constraint(m, sum(d[j] for j in covered_regions) <= γ._regional[i])
    end

    JuMP.@constraint(m, sum(d[j] for j in J) <= γ._global)

    Qrobust(m, I, J, d, q, γ)
end

function evaluate(Q::Qrobust, x::Vector{T}) where {T <: Real}
    JuMP.@objective(Q.m, Max, sum(Q.d[j] for j in Q.J) - sum(x[i]*Q.q[i] for i in Q.I))
    status = JuMP.optimize!(Q.m)
    JuMP.objective_value(Q.m), Int[round(Int,JuMP.value(d)) for d in Q.d]
end

function evaluate_objvalue(Q::Qrobust, x::Vector{T}) where {T <: Real}
    JuMP.@objective(Q.m, Max, sum(Q.d[j] for j in Q.J) - sum(x[i]*Q.q[i] for i in Q.I))
    status = JuMP.solve(Q.m)
    JuMP.objective_value(Q.m)
end

function RobustDeployment(p::DeploymentProblem; α=params.α)
     eps=params.ε
     tol=params.δ
     m = Model(GLPK.Optimizer)
     solver=Gurobi.Optimizer(OutputFlag=0, MIPGapAbs=0.9)
     #JuMP.set_optimizer_attribute(m,MOI.TimeLimitSec(), 30)
     verbose=false
     master_verbose=false
    if master_verbose
        solver=Gurobi.Optimizer(OutputFlag=1, MIPGapAbs=0.9)
        #JuMP.set_optimizer_attribute(m,MOI.TimeLimitSec(), 30)
    end
    I = 1:p.nlocations
    J = 1:p.nregions

    warmstart = naive_solution(p)

    #m = JuMP.Model(solver=solver)
    #m = Model(with_optimizer(GLPK.Optimizer))
    JuMP.@variable(m, x[i=1:p.nlocations] >= 0, Int, start=warmstart[i])
    JuMP.@variable(m, η >= 0)
    y = Vector{Matrix{JuMP.VariableRef}}()
    z = Vector{JuMP.VariableRef}()

    # Initial Restricted Master Problem
    JuMP.@objective(m, Min, η)
    JuMP.@constraint(m, sum(x[i] for i=I) <= p.nambulances)
    for j in J # coverage over all regions
        JuMP.@constraint(m, sum(x[i] for i in filter(i->p.coverage[j,i], I)) >= 1)
    end

    for i in I
        JuMP.@constraint(m, x[i] <= 5)
    end

    RobustDeployment(m, Qrobust(p, α=α, verbose=verbose), I, J, x, y, z, η,
                     Vector{Vector{Int}}(), Vector{Float64}(), Vector{Float64}(),
                     Vector{Int}[warmstart], Vector{Float64}(), Vector{Float64}())
end

function add_scenario(model::RobustDeployment, p::DeploymentProblem, scenario::Vector{T}; tol=params.δ) where {T <: Real}
    # Create variables yˡ
    push!(model.y, Array{VariableRef}(undef,(p.nlocations,p.nregions)))
    l = length(model.y)
    # for i in model.I, j in model.J
    #     model.y[l][i,j] = JuMP.VariableRef(model.m, lower_bound=0, upper_bound=p.nambulances, variable_type=:Int, base_name=String("y[$i,$j,$l]"))
    # end
    for i in model.I, j in model.J
               variable = @variable(model.m)
               set_lower_bound(variable, 0)
               set_upper_bound(variable, p.nambulances)
               set_integer(variable)
               set_name(variable, String("y[$i,$j,$l]"))
               model.y[l][i,j] = variable
           end
    push!(model.z, Array{VariableRef}(undef, p.nregions))
    # for j in model.J
    #     model.z[l][j] = JuMP.VariableRef(model.m, 0, Inf, :Int, String("z[$j,$l]"))
    # end
    for j in model.J
       variable = @variable(model.m)
       set_lower_bound(variable, 0)
       set_upper_bound(variable, Inf)
       set_integer(variable)
       set_name(variable, String("z[$j,$l]"))
       model.z[l][j] = variable
   end

    # (1) η >= 1ᵀ(dˡ + Bᴶyˡ)^+
    JuMP.@constraint(model.m, model.η >= sum(model.z[l][j] for j=model.J) + tol*sum(model.y[l][i,j] for i=model.I, j=model.J))
    for i in model.I
             outflow = JuMP.@expression(model.m, sum(model.y[l][i,j] for j in filter(j->p.coverage[j,i], model.J)))
             JuMP.@constraint(model.m, model.x[i] >= outflow)
           end
    # (2) yˡ ∈ Y(x)
    for j in model.J # shortfall from satisfying demand/calls
        inflow = JuMP.@expression(model.m, sum(model.y[l][i,j] for i in filter(i->p.coverage[j,i], model.I)))
        JuMP.@constraint(model.m, model.z[l][j] >= scenario[j] - inflow)
    end
end

function optimize_robust!(model::RobustDeployment, p::DeploymentProblem; verbose=false, maxiter=params.maxiter, eps=params.ε)
    LB = 0.0
    UB, scenario = evaluate(model.Q, model.deployment[end])
    push!(model.lowerbounds, LB)
    push!(model.upperbounds, UB)
    push!(model.scenarios, scenario)

    for k in 1:maxiter
        verbose && println("iteration $k: LB $LB, UB $UB")
        abs(UB - LB) < eps && break
        verbose && println("  solving Q with $(model.deployment[end])")

        add_scenario(model, p, scenario)
        # tic()
        status = JuMP.optimize!(model.m)
        push!(model.lowtiming)
        #@assert status == :Optimal

        LB = JuMP.objective_value(model.m)

        push!(model.deployment, [round(Int,x) for x in JuMP.value.(model.x)])

        #tic()
        shortfall, scenario = evaluate(model.Q, model.deployment[end])
        push!(model.upptiming)
        UB = min(UB, shortfall)

        # for tracking convergence later
        push!(model.upperbounds, UB)
        push!(model.scenarios, scenario)
        push!(model.lowerbounds, LB)
    end
end

#overriding optimize doesn't seem to work so I commented out this line.
#optimize!(model::RobustDeployment, p::DeploymentProblem) = JuMP.optimize!(model.m)

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
module HALeqO

export haleqo

using SolverCore, NLPModels
using LinearAlgebra, SparseArrays
using HSL, PositiveFactorizations

"""
    normsq(v)
"""
normsq(v) = dot(v, v)

"""
    merit(y, yhat, μ, fx, cx)
"""
merit(y, yhat, μ, fx, cx) =
    fx + 0.25 * μ * normsq(y) + normsq(cx + μ .* (yhat - 0.5 * y)) / μ

"""
    haleqo( nlp )
"""
function haleqo(
    nlp::AbstractNLPModel;
    x::AbstractVector = copy(nlp.meta.x0),
    y::AbstractVector = copy(nlp.meta.y0),
    tol::Real = 1e-8,
    μ::Real = 0.1,
    max_iter::Int = 3000,
    max_time::Real = 300.0,
    max_eval::Int = 100000,
    use_filter::Bool = false,
)

    start_time = time()
    nx, ny = nlp.meta.nvar, nlp.meta.ncon
    T = eltype(x)

    if has_bounds(nlp) || inequality_constrained(nlp)
        error("Problem has inequalities, can't solve it")
    end

    # parameters
    θ = 0.5
    ls_beta = 0.5
    κμminus = 0.1
    ls_eta = 1e-4
    μmin = 1e-16

    # initialization
    iter = 0
    eltime = 0.0
    status = :unknown
    fx = obj(nlp, x)
    cx = cons(nlp, x)
    yhat = copy(y)
    subres = zeros(T, nx + ny)
    subres[1:nx] .= grad(nlp, x) + jtprod(nlp, x, y)
    subres[nx+1:nx+ny] .= cx
    optimality = norm(subres[1:nx], Inf)
    cviolation = norm(cx, Inf)
    cviol__old = cviolation
    residy = cviolation
    is_solved = optimality ≤ tol && cviolation ≤ tol
    is_infeasible = false
    is_tired = iter ≥ max_iter || eltime ≥ max_time || neval_obj(nlp) > max_eval
    dir = zeros(T, nx + ny)
    xold = zeros(T, nx)
    yold = zeros(T, ny)

    rhoBM = (1.0 / μ) * max(1.0, fx) / max(1.0, 0.5*normsq(cx)) # without abs()
    rhoBM = max(1e-4, min(rhoBM, 1e4))
    μ = 1.0 / rhoBM

    if use_filter
        phi_beta = 0.1 # 0 < beta \le 1
        phi_V(cviol, optim) = cviol + phi_beta * optim
        phi_O(cviol, optim) = cviol * phi_beta + optim
        phi_theta = 0.1 # 0 < theta < 1
        phi_V_max = phi_V(cviolation, optimality)
        phi_O_max = phi_O(cviolation, optimality)
        phi_V_max *= phi_theta
        phi_O_max *= phi_theta
    end

    @info log_header(
        [:iter, :fx, :cviol, :optim, :μ, :resy],
        [Int, T, T, T, T, T],
        hdr_override = Dict(:fx => "f(x)", :optim => "‖∇L‖", :cviol => "‖c(x)‖"),
    )
    @info log_row(Any[iter, fx, cviolation, optimality])

    while !(is_solved || is_infeasible || is_tired)

        # search direction
        # compute Newton direction for the regularized sub-problem, possibly with
        # some perturbation of the Hessian matrix to obtain a descent direction
        # for the merit function
        H = hess(nlp, x, y)
        H = Symmetric(SparseMatrixCSC(Matrix(cholesky(Positive, H, Val{false}))))
        J = jac(nlp, x)

        KKT = [H J'; J UniformScaling(-μ)]
        dir .= -subres
        LDLT = Ma57(KKT)
        ma57_factorize(LDLT)
        ma57_solve!(LDLT, dir)

        # gradient of merit function
        # ∇x merit = resx + (2/μ) resy
        # ∇y merit = - resy
        xold .= subres[1:nx] + (2.0 / μ) .* jtprod(nlp, x, subres[nx+1:nx+ny])
        # slope of merit along search direction
        # slope = ∇merit ⋅ dir
        slope = dot(dir[1:nx], xold) - dot(dir[nx+1:nx+ny], subres[nx+1:nx+ny])
        if slope ≥ 0.0
            status = :not_desc
            break
        end

        # line-search
        # backtracking line-search along the primal-dual augmented Lagrangian
        # merit function with Armijo sufficient decrease condition
        xold .= x
        yold .= y
        mold = merit(y, yhat, μ, fx, cx)
        slope *= ls_eta
        τ = one(T)
        x .+= dir[1:nx]
        y .+= dir[nx+1:nx+ny]
        while true
            fx = obj(nlp, x)
            cons!(nlp, x, cx)
            m = merit(y, yhat, μ, fx, cx)
            # check Armijo's condition
            if m ≤ mold + τ * slope
                break
            end
            τ *= ls_beta
            x .= xold + τ .* dir[1:nx]
            y .= yold + τ .* dir[nx+1:nx+ny]
        end

        # evaluate residuals and check termination
        subres[1:nx] .= grad(nlp, x) + jtprod(nlp, x, y)
        optimality = norm(subres[1:nx], Inf)
        cviolation = norm(cx, Inf)
        is_solved = optimality ≤ tol && cviolation ≤ tol

        iter += 1
        if !is_solved
            eltime = time() - start_time
            is_tired = iter ≥ max_iter || eltime ≥ max_time || neval_obj(nlp) > max_eval
            if !is_tired
                is_infeasible =
                    μ < μmin && norm(jtprod(nlp, x, cx), Inf) < cviolation * √tol
                if !is_infeasible
                    subres[nx+1:nx+ny] .= cx + μ .* (yhat - y)
                    residy = norm(subres[nx+1:nx+ny], Inf)
                end
            end
        end

        # sub-problem update
        if optimality ≤ tol && residy ≤ tol
            iter_type = :M
            # check improvement in constraint violation
            if cviolation > θ * cviol__old
                # update dual regularization parameter
                μ *= κμminus
            end
            # update dual estimate
            yhat .= y
            # re-evaluate some quantities
            subres[nx+1:nx+ny] .= cx
            cviol__old = cviolation
            residy = cviolation
            @info log_row(Any[iter, fx, cviolation, optimality, μ, "$(iter_type)"])
        elseif use_filter && phi_V(cviolation, optimality) ≤ phi_V_max
            iter_type = :V
            phi_V_max *= phi_theta
            yhat .= y
            subres[nx+1:nx+ny] .= cx
            residy = cviolation
            @info log_row(Any[iter, fx, cviolation, optimality, μ, "$(iter_type)"])
        elseif use_filter && phi_O(cviolation, optimality) ≤ phi_O_max
            iter_type = :O
            phi_O_max *= phi_theta
            yhat .= y
            subres[nx+1:nx+ny] .= cx
            residy = cviolation
            @info log_row(Any[iter, fx, cviolation, optimality, μ, "$(iter_type)"])
        else
            iter_type = :F
            @info log_row(Any[iter, fx, cviolation, optimality, μ, residy])
        end

    end

    status = if is_solved
        :first_order
    elseif is_infeasible
        :infeasible
    elseif is_tired
        if iter ≥ max_iter
            :max_iter
        elseif eltime ≥ max_time
            :max_time
        elseif neval_obj(nlp) > max_eval
            :max_eval
        end
    end

    # output
    eltime = time() - start_time
    return GenericExecutionStats(
        status,
        nlp,
        solution = x,
        objective = fx,
        multipliers = y,
        dual_feas = optimality,
        primal_feas = cviolation,
        iter = iter,
        elapsed_time = eltime,
    )
end

end # module


text:
__precompile__()

module Secp256k1

using BitConverter
import Base: +, -, *, ^, /, ==, inv, sqrt, show, div
export ∞, Signature, KeyPair, ECDSA

include("lib/errors.jl")
include("lib/FieldElement.jl")
include("lib/Infinity.jl")
include("lib/Point.jl")
include("lib/scheme-types.jl")
include("lib/ECDSA.jl")

end # module


text:
const SparseBlockMap{T,B} = BlockMap{T,2,B,<:SparseMatrixCSC{B},<:NTuple{2,<:AbstractUnitRange{Int}}}

@inline SparseBlockMap(::Type{B}, block_sizes::Vararg{AbstractVector{Int}, 2}) where B =
    SparseBlockMap(B, map(blockedrange, block_sizes))

function isblockzero(A::SparseBlockMap, I::Integer, J::Integer)
    colnz = nzrange(A.blocks, J)
    isempty(colnz) || I ∉ rowvals(A.blocks)[colnz]
end
isblockzero(A::SparseBlockMap, I::BlockIndex, J::BlockIndex) =
    isblockzero(A, Int(block(I)), Int(block(J)))

function Base.getindex(A::SparseBlockMap{T}, i::Integer, j::Integer) where T
    I,J = findblockindex.(axes(A), (i,j))
    isblockzero(A,I,J) && return zero(T)
    bi,bj = Int.(block.((I,J)))
    i,j = blockindex.((I,J))
    block_getindex(A.blocks[bi,bj], i, j)
end

function SparseBlockMap(::Type{B}, block_sizes::NTuple{2,<:AbstractUnitRange{Int}}) where B
    blocks = spzeros(B, length.(block_sizes)...)
    BlockMap(blocks, block_sizes)
end

# * Application

function LinearAlgebra.mul!(w::AbstractVector{T}, A::SparseBlockMap{T}, v::AbstractVector{T},
                            α::Number=true, β::Number=false) where T
    m,n = size(A)
    size(w,1) == m && size(v,1) == n ||
        throw(DimensionMismatch("Sizes do not agree"))

    ax = axes(A)
    M,N = length.(ax)

    if iszero(β)
        w .= false # To clear NaNs
    else
        w .*= β
    end

    rows = rowvals(A.blocks)
    vals = nonzeros(A.blocks)

    for J = 1:N
        x = view(v, getindex.(ax, (1,J))[2])
        Threads.@threads for nz in nzrange(A.blocks, J)
            I = rows[nz]
            b = vals[nz]
            y = view(w, getindex.(ax, (I,J))[1])
            mul!(y, b, x, α, true)
        end
    end

    w
end

export SparseBlockMap