text:
using TextParse
using TextParse: Record, Field, Numeric, tryparsenext

# Helper type to memory map a file as a string
struct MMapString{T} <: AbstractArray{T, 1}
    # using Type parameter, too lazy to write out the result type of mmap,
    # but also don't want to loose performance
    x::T
end
MMapString(path::String) = MMapString(Mmap.mmap(open(path), Vector{UInt8}))
# we only need to overload the iteration protocol for TextParse to work!
Base.getindex(x::MMapString, i) = Char(x.x[i])
Base.next(x::MMapString, i) = Char(x.x[i]), i + 1
Base.done(x::MMapString, i) = i > length(x.x)
Base.start(x::MMapString) = 1
Base.size(x::MMapString) = size(x.x)

# Base.write creates a mutable container (Ref) for its argument to savely pass it to C.
# Even Base.unsafe_write that directly takes a pointer still has a check that allocates 128 bytes.
# Be ready for unsafe_unsafe_write, short uuwrite.
# This saves us 40 seconds && 3.5 gb less memory allocation!
# This is actually a great example of how Julia currently fails, and how one can claw back the performance.
# We just call the C library directly and use the almost deprecated syntax (`&`) to directly take a pointer from the tuple.
# in Julia 0.7 this is likely no issue and we should be able to just use `write`,
# since the compiler got a lot better at stack allocating and eliminating mutable containers.
uuwrite(io, ptr::T) where T = Int(ccall(:ios_write, Csize_t, (Ptr{Cvoid}, Ref{T}, Csize_t), io.ios, ptr, sizeof(T)))

function save2bin(path, n = typemax(Int))
    str = MMapString(path)
    # Parser descriptor of 'num,num\n' which is the format in the csv
    rec = Record((
        Field(Numeric(Float32), delim = ','),
        Field(Numeric(Float32), eoldelim = true)
    ))
    # skip the header...
    pos = findfirst(str, '\n') + 1 # Nice thing is Julia's findfirst works with any iterator
    io = open(homedir()*"/gpspoints.bin", "w")
    i = 0
    while !done(str, pos) && i <= n
        p_or_null, pos = tryparsenext(
            rec, str, pos, length(str)
        )
        isnull(p_or_null) && continue
        p = get(p_or_null)
        uuwrite(io, p)
        i += 1
    end
    close(io)
    i
end

tic()
save2bin(homedir() * "/Downloads/gps-points.csv");
toc()

"""
Transforms from longitude/latitude to pixel on screen, with `dims` refering to
the dimensions of the screen in pixel
"""
@inline function gps2pixel(point, dims)
    lon, lat = point[1], point[2]
    x = ((dims[1] / 180.0) * (90.0 + (lon / 10^7)))
    y = ((dims[2] / 360.0) * (180.0 - (lat / 10^7)))
    (x, y)
end

function to_image_inner!(img, points, start, stop)
    dims = size(img)
    for i = start:stop
        @inbounds begin
            p0 = gps2pixel(points[i], dims)
            idx = Int.(round.(p0))
            xidx, yidx = dims[1] - idx[1], dims[2] - idx[2]
            if checkbounds(Bool, img, xidx, yidx)
                # we should give each point a radius and then add the coverage to each pixel
                # for a smoother image
                # this does well enough for this short example:
                img[xidx, yidx] += 0.001f0
            end
        end
    end
end
function to_image!(img, points, range)
    N = length(range)
    NT = Threads.nthreads()
    slices = floor(Int, N / NT)
    offset = minimum(range)
    Threads.@threads for i = 1:NT
        # @threads creates a closure, which sometimes introduces type stabilities.
        # this is why it's a good practise to move the loop body behind a function barrier
        # (https://docs.julialang.org/en/latest/manual/performance-tips/#kernal-functions-1)
        to_image_inner!(img, points, offset + ((i-1) * slices + 1), offset + (i * slices))
    end
    return img
end

# Simply saving the image
using FileIO, Colors
img = zeros(Float32, 600, 960)
io = open(homedir() * "/gpspoints.bin")
points = Mmap.mmap(io, Vector{NTuple{2, Float32}})
tic()
to_image!(img, points, 1:length(points))
toc()
FileIO.save("gps.png", Gray.(clamp.(1f0 .- img, 0, 1)))
close(io)

# Or an interactive version
using Makie, Images
io = open(homedir() * "/gpspoints.bin")
# Now that we have the data as a binary blob, we can just memory map
# it as a Vector of points (NTuple{2, Float32})
points = Mmap.mmap(io, Vector{NTuple{2, Float32}})
resolution = (600, 960)
scene = Scene(resolution = reverse(resolution))
img = zeros(Float32, resolution)
imviz = heatmap(img, colornorm = (0, 1))
center!(scene)
fill!(img, 0f0)
slice = 10^7
range = slice:slice:length(points)
stop = 1
vio = VideoStream(scene, homedir()*"/Desktop/", "gps")
while true
    start = stop
    stop = min(stop + slice, length(points))
    to_image!(img, points, start:stop)
    imviz[:heatmap] = img # update image in place
    recordframe!(vio)
    stop == length(points) && break
end
finish(vio, "gif") # finish streaming and export as gif!

using Matcha

function test(str, pattern)
    res = matchat(str, 1, pattern)
    for i = 2:length(str)
        match, res = matchat(str, i, pattern)
        match && return res
    end
    res
end

mf(x) = x in ('h', 'y', 'e')

str = MMapString(homedir() * "/gpspoints.bin")

pattern = (Greed(x-> x in ('h', 'y', 'e'), 3:3),)
@time matchone(str, pattern)


text:
"""
    StaticBool(x::Bool) -> True/False

A statically typed `Bool`.
"""
abstract type StaticBool{bool} <: Integer end

struct True <: StaticBool{true} end

struct False <: StaticBool{false} end

StaticBool{true}() = True()
StaticBool{false}() = False()
StaticBool(x::StaticBool) = x
function StaticBool(x::Bool)
    if x
        return True()
    else
        return False()
    end
end

StaticInt(x::False) = Zero()
StaticInt(x::True) = One()
Base.Bool(::True) = true
Base.Bool(::False) = false

Base.:(~)(::True) = False()
Base.:(~)(::False) = True()
Base.:(!)(::True) = False()
Base.:(!)(::False) = True()

Base.:(==)(::True, ::True) = True()
Base.:(==)(::True, ::False) = False()
Base.:(==)(::False, ::True) = False()
Base.:(==)(::False, ::False) = True()

Base.:(|)(x::StaticBool, y::StaticBool) = _or(x, y)
_or(::True, ::False) = True()
_or(::False, ::True) = True()
_or(::True, ::True) = True()
_or(::False, ::False) = False()
Base.:(|)(x::Bool, y::True) = y
Base.:(|)(x::Bool, y::False) = x
Base.:(|)(x::True, y::Bool) = x
Base.:(|)(x::False, y::Bool) = y

Base.:(&)(x::StaticBool, y::StaticBool) = _and(x, y)
_and(::True, ::False) = False()
_and(::False, ::True) = False()
_and(::True, ::True) = True()
_and(::False, ::False) = False()
Base.:(&)(x::Bool, y::True) = x
Base.:(&)(x::Bool, y::False) = y
Base.:(&)(x::True, y::Bool) = y
Base.:(&)(x::False, y::Bool) = x

Base.xor(y::StaticBool, x::StaticBool) = _xor(x, y)
_xor(::True, ::True) = False()
_xor(::True, ::False) = True()
_xor(::False, ::True) = True()
_xor(::False, ::False) = False()
Base.xor(x::Bool, y::StaticBool) = xor(x, Bool(y))
Base.xor(x::StaticBool, y::Bool) = xor(Bool(x), y)

Base.sign(x::StaticBool) = x
Base.abs(x::StaticBool) = x
Base.abs2(x::StaticBool) = x
Base.iszero(::True) = False()
Base.iszero(::False) = True()
Base.isone(::True) = True()
Base.isone(::False) = False()

Base.:(<)(x::StaticBool, y::StaticBool) = _lt(x, y)
_lt(::False, ::True) = True()
_lt(::True, ::True) = False()
_lt(::False, ::False) = False()
_lt(::True, ::False) = False()

Base.:(<=)(x::StaticBool, y::StaticBool) = _lteq(x, y)
_lteq(::False, ::True) = True()
_lteq(::True, ::True) = True()
_lteq(::False, ::False) = True()
_lteq(::True, ::False) = False()

Base.:(+)(x::True) = One()
Base.:(+)(x::False) = Zero()
Base.:(-)(x::True) = -One()
Base.:(-)(x::False) = Zero()

Base.:(+)(x::StaticBool, y::StaticBool) = StaticInt(x) + StaticInt(y)
Base.:(-)(x::StaticBool, y::StaticBool) = StaticInt(x) - StaticInt(y)
Base.:(*)(x::StaticBool, y::StaticBool) = x & y

# from `^(x::Bool, y::Bool) = x | !y`
Base.:(^)(x::StaticBool, y::False) = True()
Base.:(^)(x::StaticBool, y::True) = x
Base.:(^)(x::Integer, y::False) = one(x)
Base.:(^)(x::Integer, y::True) = x
Base.:(^)(x::BigInt, y::False) = one(x)
Base.:(^)(x::BigInt, y::True) = x

Base.div(x::StaticBool, y::False) = throw(DivideError())
Base.div(x::StaticBool, y::True) = x

Base.rem(x::StaticBool, y::False) = throw(DivideError())
Base.rem(x::StaticBool, y::True) = False()
Base.mod(x::StaticBool, y::StaticBool) = rem(x, y)

Base.promote_rule(::Type{<:StaticBool}, ::Type{<:StaticBool}) = StaticBool
Base.promote_rule(::Type{<:StaticBool}, ::Type{Bool}) = Bool
Base.promote_rule(::Type{Bool}, ::Type{<:StaticBool}) = Bool

@generated _get_tuple(::Type{T}, ::StaticInt{i}) where {T<:Tuple, i} = T.parameters[i]

Base.all(::Tuple{Vararg{True}}) = true
Base.all(::Tuple{Vararg{Union{True,False}}}) = false
Base.all(::Tuple{Vararg{False}}) = false

Base.any(::Tuple{Vararg{True}}) = true
Base.any(::Tuple{Vararg{Union{True,False}}}) = true
Base.any(::Tuple{Vararg{False}}) = false

ifelse(::True, x, y) = x

ifelse(::False, x, y) = y

Base.show(io::IO, ::StaticBool{bool}) where {bool} = print(io, "static($bool)")

ext: jl


text:
using MultiStochGrad, Test

using Logging
using Base.CoreLogging

logger = ConsoleLogger(stdout, CoreLogging.Info)
global_logger(logger)

include("TestLinearRegression.jl")

@testset "linreg" begin
@test test_linear_regression_scsg()
@test test_linear_regression_svrg()
end

# possibly  after download of mnist files
include("TestLogisticRegression.jl")

@testset "logistic_reg" begin
@test mnist_logistic_regression_svrg()
@test mnist_logistic_regression_scsg()
end


text:
using DeepForest, Flux
using StatsBase: mode, mean
using Test

@testset "Forest init" begin
    f = Forest(10, 2, 10, 10, 0.25)
    @test length(f.trees) == 10
    @test length(f.trees[1].subset) == 3
end

@testset "Forest precompute" begin
    x = rand(Float32, 3, 4)
    y = [1, 2, 2, 2]
    f = Forest(10, 2, 10, 3, 2/3)
    DeepForest.precompute!(f, x, y)
    @test f.trees[1].best.first != 1 || f.trees[1].best.second != 1
end

@testset "Forest predict" begin
    x = rand(Float32, 3, 4)
    y = [2, 2, 2, 2]
    f = Forest(10, 2, 10, 3, 2/3)
    DeepForest.precompute!(f, x, y)
    @test predict(f, x) == [2, 2, 2, 2]
end

@testset "Forest loss" begin
    x = rand(Float32, 3, 4)
    y = [1, 2, 2, 2]
    f = Forest(10, 2, 10, 3, 2/3)
    DeepForest.precompute!(f, x, y)

    @test loss!(f, x, y) != 0

    gs = gradient(() -> loss!(f, x, y), params(f))
    @test gs.grads != IdDict()

    FS = DeepForest.ForestSlice(f, Array(1:5))
    gs = gradient(() -> loss!(FS, x, y), params(f))
    @test gs.grads != IdDict()
end

@testset "Forest train" begin
    x = rand(Float32, 2, 1000);
    x[1, :] .*= 2*π;
    x[1, :] .-= π;
    x[2, :] .*= 3;
    x[2, :] .-= 1.5;
    y = Int64.(x[2, :] .< sin.(x[1, :])) .+ 1;

    model = Forest(25, 2, 10, 2, 1.0)

    forest_train!(1500, model, x, y)
    @test mean(predict(model, x) .== y) > .9
end

@testset "Forest importance" begin
    x = rand(Float32, 2, 1000);
    x[1, :] .*= 2*π;
    x[1, :] .-= π;
    x[2, :] .*= 3;
    x[2, :] .-= 1.5;
    y = Int64.(x[2, :] .< sin.(x[1, :])) .+ 1;

    model = Forest(25, 2, 10, 2, 1.0)

    forest_train!(1500, model, x, y)
    imp = importance(model, x)
    @test imp[2] > 0
    @test imp[1] > 0
end


text:
export Origami
module Origami

export evaluate, Point, Fold

using Base.Iterators

struct Point{T}
    x::T
    y::T

    function Point(x::S, y::T) where {S,T}
        R = promote_type(S, T)
        R = typeof(R(1) / R(1))
        new{R}(convert(R, x), convert(R, y))
    end
    Point{S}(x::T, y::T) where {S,T} = Point(convert(S, x), convert(S, y))
end
import Base: -
-(a::Point, b::Point) = Point(a.x - b.x, a.y - b.y)

const ArrayOrTuple{X} = Union{AbstractArray{X},NTuple{N,X} where N}
const Shape{T} = ArrayOrTuple{Point{T}}
const Fold{T} = Tuple{Point{T},Point{T}}
const Folds{T} = ArrayOrTuple{Fold{T}}

function evaluate(shape::Shape{T}, folds, probe::Point{T}) where T
    if isempty(folds)
        Int(inshape(probe, shape))
    else
        fold = first(folds)
        shape1, shape2 = split_and_fold(shape, fold)
        folds2 = drop(folds, 1)
        evaluate(shape1, folds2, probe) + evaluate(shape2, folds2, probe)
    end
end

function inshape(point::Point, shape::Shape)
    length(shape) <= 3 && return false
    y = point.y
    p0 = shape[1]
    r = 0
    for k = 2:length(shape)
        p1 = shape[k]
        if min(p0.y, p1.y) <= y <= max(p1.y, p0.y)
            x = atsideof(point, (p0, p1))
            (x == 0 || 0 != r != x) && return false
            r = x
        end
        p0 = p1
    end
    r != 0
end

split_and_fold(shape::Shape, fold::Fold) = foldifright.(split(shape, fold), Ref(fold))
foldifright(shape::Shape, fold::Fold) = isrightof(shape, fold) ? folded(shape, fold) : shape

function isrightof(shape::Shape, fold)
    all(atsideof.(shape, Ref(fold)) .>= 0)
end

# returns -1: left, 1: right, 0: on line
function atsideof(point::Point, fold::Fold)
    p1, p2 = fold
    p= p2 - p1
    q = point - p1
    a = p.y * q.x - p.x * q.y
    a * a * 10^12 > (p.x^2 + p.y^2) ? (a < 0 ? -1 : 1) : 0
end

function folded(shape::Shape, fold::Fold)
    isempty(shape) ? shape : folded.(shape, Ref(fold))
end

function folded(point::Point, fold::Fold)
    A, p = matrix_vector(fold)
    v = A * ([point.x; point.y] - p) + p
    Point(v...)
end

function matrix_vector(fold::Fold)
    p1, p2 = fold
    dp = p2 - p1
    x, y = dp.x, dp.y
    c = x^2 + y^2
    a = (x^2 - y^2) / c
    b = (x * y * 2) / c
    [a b; b -a], [p1.x; p1.y]
end

function split(shape::Shape, fold::Fold)
    isempty(shape) && return shape, shape
    v = shape[1]
    side0 = atsideof(v, fold)
    shape1 = empty(shape)
    shape2 = empty(shape)
    switch(sh) = sh === shape1 ? shape2 : shape1
    sh = shape1
    push!(sh, v)
    for k = 2:length(shape)
        w = shape[k]
        side1 = atsideof(w, fold)
        if side1 == side0
            push!(sh, w)
        elseif side0 == 0
            side0 = side1
            push!(sh, w)
        elseif side1 == 0
            push!(sh, w)
            sh = switch(sh)
            push!(sh, w)
            side0 = -side0
        else
            u = intersect(v, w, fold)
            push!(sh, u)
            sh = switch(sh)
            push!(sh, u)
            push!(sh, w)
            side0 = side1
        end
        v = w
    end
    normal!(shape1), normal!(shape2)
end

function normal!(shape::AbstractArray)
    if length(shape) > 1 && shape[1] !== shape[end]
        push!(shape, shape[1])
    end
    if length(shape) <= 3
        empty!(shape)
    end
    shape
end

function intersect(x0::Point, x1::Point, fold::Fold)
    p0, p1 = fold
    dp = p1 - p0
    dx = x1 - x0
    b = p0 - x0
    a, _ = [dx.x dp.x; dx.y dp.y] \ [b.x; b.y]
    Point(a * dx.x + x0.x, a * dx.y + x0.y)
end

end # module


text:
b₀(ξ) = abs(ξ) < 0.5 ? 1 : 0
b₁(ξ) = (ξ >= -1 && ξ <= 1) ? (1 - abs(ξ)) : 0
b₂(ξ) = (ξ >= -1.5 && ξ <= 1.5) ?
    ((ξ >= -0.5 && ξ <= 0.5) ? (0.75 - ξ^2) : (0.5 * (1.5 - abs(ξ))^2)) : 0

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
HTTP.@register(ROUTER, "GET", "/fileio/qe/parse/output/pw", req -> FileIO.qe_read_pw_output(JSON3.read(req.body, String)))
HTTP.@register(ROUTER, "GET", "/fileio/qe/parse/output/projwfc", req -> FileIO.qe_read_projwfc_output(JSON3.read(req.body, String)))
HTTP.@register(ROUTER, "GET", "/fileio/qe/parse/output/hp", req -> FileIO.qe_read_hp_output(JSON3.read(req.body, String)))
HTTP.@register(ROUTER, "GET", "/fileio/qe/parse/input", req -> FileIO.qe_read_calculation(JSON3.read(req.body, String)))

HTTP.@register(ROUTER, "GET", "/fileio/wannier90/parse/input", req -> FileIO.wan_read_calculation(JSON3.read(req.body, String)))
HTTP.@register(ROUTER, "GET", "/fileio/wannier90/parse/output", req -> FileIO.wan_read_output(JSON3.read(req.body, String)))


text:
"""
         simulate(p    :: Drip;
                  T    :: Integer = 500,                    # Optional: length of simulation
                  burn :: Integer = 100,                    # Optional: length of initial burn (in addition to T)
                  N    :: Union{Int4,Nothing} = nothing,    # Optional: number of simulated agents (returns the average beliefs of a large set of agents by default)
                  seed :: Union{Integer,Nothing} = nothing  # Optional: seed number for fundamental shocks
                  ) -> Path
Returns a `Path` structure with a simulated path of the fundamental (`x`), beliefs (`x_hat`)
    and actions (`a`)
    **under the steady state information structure**. In particular, if `n` is the
    dimension of `x` and `m` is the dimension of `a`, then

* `x` has dimension `n*T` where `x(:,t)` is the simulated value of `x` at time `t`.
* `x_hat` has dimension `n*T*N` where `x_hat(:,i,t)` is the simulated value of `x_hat` of agent `i` at time `t`.
* `a` has dimension `m*T*N` where `a(:,i,t)` is the the simulated value of `x_hat` of agent `i` at time `t`.
"""
function simulate(p    :: Drip;
                  T    :: Integer = 500,                    # Optional: length of simulation
                  burn :: Integer = 100,                    # Optional: length of initial burn (in addition to T)
                  N    :: Union{Integer,Nothing} = nothing, # Optional: number of simulated agents (returns the average beliefs of a large set of agents by default)
                  seed :: Union{Integer,Nothing} = nothing  # Optional: seed number for fundamental shocks
                  )
    # create flag for if N is nothing and change N to 1
    average_belief = (N == nothing)
    if average_belief; N = 1; end
    # get dimensions
    (n,m)    = size(p.H)
    (_,k)    = size(p.Q)
    # set seed for fundamental shocks
    seed != nothing ? Random.seed!(seed) : Random.seed!()
    x_Shocks = randn(k,burn+T);

    Random.seed!()
    a_Shocks = average_belief ? zeros(m,burn+T,N) : randn(m,burn+T,N);
    x        = zeros(n,burn+T);
    x_hat    = zeros(n,burn+T,N);
    a        = zeros(m,burn+T,N);

    SqΣz = average_belief ? 0 : sqrt(p.ss.Σ_z);

    for t in 2:T+burn
        x[:,t] = p.A*x[:,t-1]+p.Q*x_Shocks[:,t];
        for i in 1:N
            x_hat[:,t,i] = p.A*x_hat[:,t-1,i]+(p.ss.K*p.ss.Y')*(x[:,t]-p.A*x_hat[:,t-1,i])+p.ss.K*SqΣz*a_Shocks[:,t,i];
            a[:,t,i]    .= p.H'*x_hat[:,t,i];
        end
    end
    x     = x[:,burn+1:T+burn];
    x_hat = average_belief ? x_hat[:,burn+1:T+burn] : x_hat[:,burn+1:T+burn,:];
    a = average_belief ? a[:,burn+1:T+burn] : a[:,burn+1:T+burn,:];
    return(Path(T,x,x_hat,a))
end


text:
#InspectDR: Top/high-level functionnality of Gtk layer
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
#Hack to keep text height in status label from changing as exponents start
#being included in coordinate values (Noticed on Windows platforms):
#(Causes display glitches/slowdowns)
const TEXTH_HACK_STR = "⁻" #Superscript `-` has bigger height than other chars.
#TODO: find a better way to keep height constant in status bar.

const COORDSTAT_NONE = "$TEXTH_HACK_STR(x, y) = ( , )"

#==Callback wrapper functions (PlotWidget-level)
===============================================================================#
#=COMMENTS
 -Define callback wrapper functions with concrete types to assist with precompile.
 -Remove unecessary arguments/convert to proper types
=#
@guarded function cb_keypress(w::Ptr{Gtk.GObject}, event::Gtk.GdkEventKey, pwidget::PlotWidget)
	handleevent_keypress(pwidget.state, pwidget, event)
	nothing #Known value
end
@guarded function cb_scalechanged(w::Ptr{Gtk.GObject}, pwidget::PlotWidget)
	handleevent_scalechanged(pwidget.state, pwidget)
	nothing #Known value
end
@guarded function cb_mousepress(w::Ptr{Gtk.GObject}, event::Gtk.GdkEventButton, pwidget::PlotWidget)
	handleevent_mousepress(pwidget.state, pwidget, event)
	nothing #Known value
end
@guarded function cb_mouserelease(w::Ptr{Gtk.GObject}, event::Gtk.GdkEventButton, pwidget::PlotWidget)
	handleevent_mouserelease(pwidget.state, pwidget, event)
	nothing #Known value
end
@guarded function cb_mousemove(w::Ptr{Gtk.GObject}, event::Gtk.GdkEventMotion, pwidget::PlotWidget)
	handleevent_mousemove(pwidget.state, pwidget, event)
	nothing #Known value
end
@guarded function cb_mousescroll(w::Ptr{Gtk.GObject}, event::Gtk.GdkEventScroll, pwidget::PlotWidget)
	handleevent_mousescroll(pwidget.state, pwidget, event)
	nothing #Known value
end

#==Callback wrapper functions (GtkPlot-level)
===============================================================================#
@guarded function cb_wnddestroyed(w::Ptr{Gtk.GObject}, gplot::GtkPlot)
	gplot.destroyed = true
	nothing #Known value
end
@guarded function cb_mnufileexport(w::Ptr{Gtk.GObject}, gplot::GtkPlot)
	filepath = Gtk.save_dialog("Export plot...", _Gtk.Null(),
		(_Gtk.FileFilter("*.png,*.svg,*.eps", name="All supported formats"), "*.png", "*.svg", "*.eps")
	)
	if isempty(filepath); return nothing; end

	ext = splitext(filepath)[end]
	mime = get(MAPEXT2MIME, ext, nothing)
	if nothing == mime
		Gtk.warn_dialog("Unrecognized file type: '$ext'")
		return nothing
	end

	try
		_write(filepath, mime, gplot)
	catch
		Gtk.warn_dialog("Write failed: '$filepath'")
	end
	nothing #Known value
end
@guarded function cb_mnufileclose(w::Ptr{Gtk.GObject}, gplot::GtkPlot)
	window_close(gplot.wnd)
	nothing #Known value
end

@guarded function cb_mnudatatraces(w::Ptr{Gtk.GObject}, gplot::GtkPlot)
	tracedialog_show(gplot)
	nothing #Known value
end

#==Higher-level event handlers
===============================================================================#
#plothover event: show plot coordinates under mouse.
#-------------------------------------------------------------------------------
function handleevent_plothover(gplot::GtkPlot, pwidget::PlotWidget, x::Float64, y::Float64)
	plotinfo = pwidget.plotinfo #WANTCONST

	istrip = pwidget.mouseover.istrip
	pos = pwidget.mouseover.pos
	if istrip > 0 && pos != nothing
		xfmt = hoverfmt(plotinfo.xfmt)
		yfmt = hoverfmt(plotinfo.strips[istrip].yfmt)
		xstr = formatted(pos.x, xfmt)
		ystr = formatted(pos.y, yfmt)
		statstr = "$TEXTH_HACK_STR(x, y) = ($xstr, $ystr)"
	else
		statstr = COORDSTAT_NONE
	end

	set_gtk_property!(gplot.status, :label, statstr)
	nothing
end

#==Menu builders:
===============================================================================#
function Gtk_addmenu(parent::Union{_Gtk.Menu, _Gtk.MenuBar}, name::String)
	item = _Gtk.MenuItem(name)
	mnu = _Gtk.Menu(item)
	push!(parent, item)
	return mnu
end
function Gtk_addmenuitem(mnu::_Gtk.Menu, name::String)
	item = _Gtk.MenuItem(name)
	push!(mnu, item)
	return item
end

#=="Constructors"
===============================================================================#

#-------------------------------------------------------------------------------
function PlotWidget(plot::Plot)
	vbox = _Gtk.Box(true, 0)
		can_focus(vbox, true)
#		set_gtk_property!(vbox, "focus-on-click", true)
#		set_gtk_property!(vbox, :focus_on_click, true)
	canvas = Gtk.Canvas()
		set_gtk_property!(canvas, :vexpand, true)
	w_xscale = _Gtk.Scale(false, 1:XAXIS_SCALEMAX)
		xscale = _Gtk.Adjustment(w_xscale)
		set_gtk_property!(xscale, :value, 1)
#		draw_value(w_xscale, false)
		value_pos(w_xscale, Int(GtkPositionType.GTK_POS_RIGHT))
	w_xpos = _Gtk.Scale(false, -.5:XAXIS_POS_STEPRES:.5)
		xpos = _Gtk.Adjustment(w_xpos)
		set_gtk_property!(xpos, :value, 0)
#		draw_value(w_xpos, false)
		value_pos(w_xpos, Int(GtkPositionType.GTK_POS_RIGHT))

	push!(vbox, canvas)
	#push!(vbox, w_xpos)
	#push!(vbox, w_xscale)

	w = width(canvas); h = height(canvas)
	plotbuf = CairoBufferedPlot(
		Cairo.CairoARGBSurface(w, h), Cairo.CairoARGBSurface(w, h)
	)
	#TODO: how do we get a maximum surface for all monitors?
	#TODO: or can we resize in some intelligent way??
	#plotbuf = Cairo.CairoRGBSurface(1920,1200) #Appears slow for average monitor size???
#	plotbuf = Gtk.cairo_surface_for(canvas) #create similar - does not work here
	curstrip = 1 #TODO: Is this what is desired?

	plotinfo = Plot2DInfo(plot)
	pwidget = PlotWidget(vbox, canvas, plot, plotinfo, ISNormal(),
		w_xscale, xscale, w_xpos, xpos,
		plotbuf, curstrip, GtkMouseOver(),
		CtrlMarkerGroup(plot.layout[:font_annotation]), nothing, true, true,
		#Event handlers:
		nothing
	)

	#Take control of parentannot:
	plot.parentannot = [pwidget.markers]

	#Register callback functions:
	signal_connect(cb_scalechanged, xscale, "value-changed", Nothing, (), false, pwidget)
	signal_connect(cb_scalechanged, xpos, "value-changed", Nothing, (), false, pwidget)
	signal_connect(cb_keypress, vbox, "key-press-event", Nothing, (Ref{Gtk.GdkEventKey},), false, pwidget)
	signal_connect(cb_mousepress, vbox, "button-press-event", Nothing, (Ref{Gtk.GdkEventButton},), false, pwidget)
	signal_connect(cb_mouserelease, vbox, "button-release-event", Nothing, (Ref{Gtk.GdkEventButton},), false, pwidget)
	signal_connect(cb_mousemove, vbox, "motion-notify-event", Nothing, (Ref{Gtk.GdkEventMotion},), false, pwidget)
	signal_connect(cb_mousescroll, vbox, "scroll-event", Nothing, (Ref{Gtk.GdkEventScroll},), false, pwidget)

	#Register event: draw function
	Gtk.@guarded Gtk.draw(pwidget.canvas) do canvas
		#=NOTE:
		PlotWidget should probably be subclassed from GtkCanvas (the draw event
		would then have a reference to the PlotWidget...), but Julia does not
		make this easy.  Instead, this function generates an annonymous function
		that implicitly has a reference to the appropriate PlotWidget instance.
		=#
		if invalidbuffersize(pwidget)
			render(pwidget)
		end
		ctx = getgc(pwidget.canvas)
			w = width(pwidget.canvas); h = height(pwidget.canvas)
			bb = BoundingBox(0, w, 0, h)
		wipe(ctx, bb, COLOR_WHITE)
		Cairo.set_source_surface(ctx, pwidget.plotbuf.surf, 0, 0)
		Cairo.paint(ctx) #Applies contents of plotbuf.surf
		selectionbox_draw(ctx, pwidget) #TODO: use more generic method name??
		#TODO: Can/should we explicitly Cairo.destroy(ctx)???
	end

	return pwidget
end
#Build a PlotWidget & register event handlers for GtkPlot object
function PlotWidget(gplot::GtkPlot, plot::Plot)
	pwidget = PlotWidget(plot)
	pwidget.eh_plothover = HandlerInfo(gplot, handleevent_plothover)
	return pwidget
end

#Synchronize with gplot.src.subplots
function sync_subplots(gplot::GtkPlot)
	wlist = gplot.subplots #widget list
	plist = gplot.src.subplots #(sub)plot list
	resize!(wlist, length(plist))

	for (i, s) in enumerate(plist)
		if !isassigned(wlist, i)
			wlist[i] = PlotWidget(gplot, s)
		else
			if wlist[i].src != s
				Gtk.destroy(wlist[i].widget)
				wlist[i] = PlotWidget(gplot, s)
			end
		end
	end

	#Blindly re-construct grid:
	for i in length(gplot.grd):-1:1
		Gtk.delete!(gplot.grd, gplot.grd[i]) #Does not destroy existing child widgets
	end
	ncols = gplot.src.layout.values.ncolumns #WANTCONST
	for (i, w) in enumerate(wlist)
		row = div(i-1, ncols)+1
		col = i - ((row-1)*ncols)
		gplot.grd[col,row] = w.widget

		#FIXME/HACK: rebuilding grd appears to inhibit the redraw mechanism.
		#Toggling w.canvas -> visible unclogs refresh algorithm somehow.
		set_gtk_property!(w.canvas, :visible, false)
		set_gtk_property!(w.canvas, :visible, true)
	end
	return
end

#-------------------------------------------------------------------------------
function GtkPlot(mp::Multiplot)
	#Generate graphical elements:
	mb = _Gtk.MenuBar()
	mnufile = Gtk_addmenu(mb, "_File")
		mnuexport = Gtk_addmenuitem(mnufile, "_Export")
		push!(mnufile, _Gtk.SeparatorMenuItem())
		mnuquit = Gtk_addmenuitem(mnufile, "_Quit")
	mnudata = Gtk_addmenu(mb, "_Data")
		mnutraces = Gtk_addmenuitem(mnudata, "_Traces")
	grd = Gtk.Grid() #Main grid with different subplots.
		set_gtk_property!(grd, :column_homogeneous, true)
		#set_gtk_property!(grd, :column_spacing, 15) #Gap between
	status = _Gtk.Label("")#"(x,y) =")
		set_gtk_property!(status, :hexpand, true)
		set_gtk_property!(status, :ellipsize, PANGO_ELLIPSIZE_END)
		set_gtk_property!(status, :xalign, 0.0)
		set_gtk_property!(status, :label, COORDSTAT_NONE)
		sbar_frame = _Gtk.Frame(status)
			set_gtk_property!(sbar_frame, "shadow-type", GtkShadowType.GTK_SHADOW_ETCHED_IN)

	vbox = _Gtk.Box(true, 0)
		push!(vbox, mb) #Menu bar
		push!(vbox, grd) #Subplots
		push!(vbox, sbar_frame) #status bar
	wnd = Gtk.Window(vbox, "", 640, 480, true)

	gplot = GtkPlot(false, wnd, grd, [], mp, status)
	refresh_title(gplot)
	sync_subplots(gplot)

	if length(gplot.subplots) > 0
		set_focus(wnd, gplot.subplots[end].widget)
	end

	Gtk.showall(wnd)
	signal_connect(cb_wnddestroyed, wnd, "destroy", Nothing, (), false, gplot)
	signal_connect(cb_mnufileexport, mnuexport, "activate", Nothing, (), false, gplot)
	signal_connect(cb_mnufileclose, mnuquit, "activate", Nothing, (), false, gplot)
	signal_connect(cb_mnudatatraces, mnutraces, "activate", Nothing, (), false, gplot)

	return gplot
end

#-------------------------------------------------------------------------------
function GtkPlot(plot::Plot, args...; kwargs...)
	mp = Multiplot(args...; kwargs...)
	_add(mp, plot)
	return GtkPlot(mp)
end
GtkPlot(args...; kwargs...) = GtkPlot(Plot2D(), args...; kwargs...)

#==High-level interface
===============================================================================#
function clearsubplots(gplot::GtkPlot)
	for s in gplot.subplots
		Gtk.destroy(s.widget)
	end
	gplot.subplots = []
	gplot.src.subplots = []
	return gplot
end

function refresh(gplot::GtkPlot)
	if !gplot.destroyed
		refresh_title(gplot)
		active = active_subplot(gplot)
		set_gtk_property!(gplot.grd, :visible, false) #Suppress gliching
			sync_subplots(gplot)
			map(refresh, gplot.subplots) #Is this necessary?
		set_gtk_property!(gplot.grd, :visible, true)
		set_focus(gplot.subplots[active].widget) #Restore focus after modifying gplot.grd
		Gtk.showall(gplot.grd)
		#TODO: find a way to force GUI to updates here... Animations don't refresh...
		sleep(eps(0.0)) #Ugly Hack: No guarantee this works... There must be a better way.
	end
	return
end

function Base.display(d::GtkDisplay, mp::Multiplot)
	return GtkPlot(mp)
end
function Base.display(d::GtkDisplay, p::Plot)
	return GtkPlot(p)
end

#Last line


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

using .GeoStatsBase
import .GeoStatsBase: preprocess, solve_single

export TuringPat, PARAMS1

const PARAMS1 = [
  Params(2,  4,   0.01),
  Params(5,  10,  0.02),
  Params(10, 20,  0.03),
  Params(20, 40,  0.04),
  Params(50, 100, 0.05)
]

"""
    TuringPat(var₁=>param₁, var₂=>param₂, ...)

Turing pattern simulation.

## Parameters

* `params` - basic parameters (default to `PARAMS1`)
* `blur`   - blur algorithm (default to `BoxBlur`)
* `edge`   - edge condition (default to `Clamp`)
* `iter`   - number of iterations (default to `100`)

### References

Turing 1952. *The chemical basis of morphogenesis.*
"""
@simsolver TuringPat begin
  @param params = PARAMS1
  @param blur = BoxBlur
  @param edge = Clamp
  @param iter = 100
end

function preprocess(problem::SimulationProblem, solver::TuringPat)
  # retrieve domain size
  sz = size(domain(problem))

  # result of preprocessing
  preproc = Dict{Symbol,NamedTuple}()

  for (var, V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = TuringPatParam()
    end

    # determine simulation parameters
    params = varparams.params
    blur = varparams.blur(sz)
    edge = varparams.edge()
    iter = varparams.iter

    # construct patterns from parameters
    patterns = [SimplePattern(param, sz) for param in params]

    preproc[var] = (patterns=patterns,blur=blur,edge=edge,iter=iter)
  end

  preproc
end

function solve_single(problem::SimulationProblem, var::Symbol,
                      solver::TuringPat, preproc)
  # retrieve domain size
  sz = size(domain(problem))

  # unpack preprocessed parameters
  patterns, blur, edge, iter = preproc[var]

  # determine result type
  V = variables(problem)[var]

  # perform simulation
  sim = Sim(rand(V, sz), patterns, edge, blur)
  for i in 1:iter
    step!(sim)
  end
  real = scale01(sim.fluid)

  # flatten result
  vec(real)
end


text:
@testset "Transformed cubes" begin
    @testset "Simple map" begin
        data = collect(reshape(1:20,4,5))
        axlist = [RangeAxis("XVals",1.0:4.0), CategoricalAxis("YVals",[1,2,3,4,5])]
        props = Dict("att1"=>5, "att2"=>"Hallo")
        a = YAXArray(axlist, data, props)
        a2 = map(i->i+1, a)
        @test all(a2.data .== a.data .+1)
        @test a2.axes == a.axes
        @test a2.properties==a.properties
        #TODO: Test that cleaners get promoted properly
        a3 = map(+,a,a2)
        @test all(a3.data .== a.data .+ a2.data)
        @test a3.axes == a.axes
        @test a3.properties==a.properties
    end
    @testset "Cube concatenation" begin
        data = [rand(4,5) for i=1:3]
        axlist = [RangeAxis("XVals",1.0:4.0), CategoricalAxis("YVals",[1,2,3,4,5])]
        props = [Dict("att$i"=>i) for i=1:3]
        singlecubes = [YAXArray(axlist, data[i], props[i]) for i=1:3]
        newcube = concatenatecubes(singlecubes, CategoricalAxis("ZVals",["A","B","C"]))
        @test caxes(newcube) == [axlist;CategoricalAxis("ZVals",["A","B","C"])]
        @test ndims(newcube) == 3
        @test size(newcube) == (4,5,3)
        @test newcube[:,:,:] == cat(data...,dims=3)
        @test getattributes(newcube) == reduce(merge,props)
    end
end


text:
module AnytimeWeightedAStar

export awastar_search, rwastar_search, awastar_search_with_scheduled_weights, AWAStar, simulated_time, wall_time, expansion_rate, quality

include("search_problem/SearchProblem.jl")
include("graph_search/GraphSearch.jl")
include("example_problems/ExampleProblems.jl")

include("utils.jl")
include("priority_dict.jl")
include("multiple_open_lists.jl")
include("awastar.jl")
include("rwastar.jl")
include("awastar_ws.jl")

end # module


text:
@testset "LSA" begin
    # Documents
    doc1 = StringDocument("This is a text about an apple. There are many texts about apples.")
    doc2 = StringDocument("Pears and apples are good but not exotic. An apple a day keeps the doctor away.")
    doc3 = StringDocument("Fruits are good for you.")
    doc4 = StringDocument("This phrase has nothing to do with the others...")
    doc5 = StringDocument("Simple text, little info inside")
    # Corpus
    crps = Corpus(AbstractDocument[doc1, doc2, doc3, doc4, doc5])
    prepare!(crps, strip_punctuation)
    ngram_complexity = 2
    lex = create_lexicon(crps, ngram_complexity)
    lex_keys = collect(keys(lex))
    n = length(crps)
    # Retrieval
    query = StringDocument("Apples and an exotic fruit.")
    for k in [1, 6]
        for stats in [:count, :tf, :tfidf, :bm25]
            for T in [Float32, Float64]
                dtm = DocumentTermMatrix{T}(crps, lex_keys, ngram_complexity=ngram_complexity)
                model = lsa(dtm, k=k, stats=stats, ngram_complexity=ngram_complexity)
                @test model isa LSAModel{String, T, Matrix{T}, Int}
                @test length(model.Σinv) == min(k, n)
                idxs, corrs = cosine(model, dtm, query)
                @test length(idxs) == length(corrs) == length(crps)
                sim = similarity(model, crps[rand(1:n)], query)
                @test -1.0 <= sim <= 1.0
            end
        end
    end

    # Tests for the rest of the functions
    K = 2
    T = Float32
    # Vocabulary
    model = lsa(crps, T, k=K, ngram_complexity=ngram_complexity)
    @test model isa LSAModel{String, T, Matrix{T}, Int}
    @test all(in_vocabulary(model, word) for word in lex_keys)
    @test vocabulary(model) == lex_keys
    @test size(model) == (length(lex_keys), K)
    # Document, corpus embedding
    dtm = DocumentTermMatrix{T}(crps, lex_keys, ngram_complexity=ngram_complexity)
    V = Matrix(embed_document(model, crps))
    @test eltype(V) == T
    @test eltype(embed_document(model, crps[1])) == T
    for i in 1:n
        @test all(V[:,i] .≈ embed_document(model, crps[i]))
    end
    # Index, Similarity
    idx = 2
    word = model.vocab[idx]
    @test index(model, word) == model.vocab_hash[word]
    @test get_vector(model, word) == model.Uᵀ[:, idx]
    @test similarity(model, crps[1], crps[2]) isa T
    @test similarity(model, crps[1], crps[2]) == similarity(model, crps[2], crps[1])
    @test_throws ErrorException LSAModel(DocumentTermMatrix{Int}(crps), k=K)
    @test_throws ErrorException StringAnalysis.embed_word(model, "word")
    # Test saving and loading an LSA model
    T = Float32
    vocab = split("a random string")
    vocab_hash = OrderedDict("a"=>1, "random"=>2, "string"=>3)
    k=2
    stats = :tf
    κ = 2
    β = 0.71
    tol = 1e-5
    model = lsa(crps, T, k=k, stats=stats, tol=T(tol), κ=κ, β=β)
    file = "./_lsa_model.txt"
    save_lsa_model(model, file)
    # Model 1
    loaded_model_1 = load_lsa_model(file, Float64, sparse=true)
    @test loaded_model_1 isa LSAModel{String, Float64, SparseMatrixCSC{Float64, Int}, Int}
    @test all(loaded_model_1.Uᵀ .≈ sparse(model.Uᵀ))
    @test all(loaded_model_1.Σinv .≈ sparse(model.Σinv))
    @test loaded_model_1.vocab == model.vocab
    @test loaded_model_1.vocab_hash == model.vocab_hash
    @test loaded_model_1.stats == model.stats
    @test all(loaded_model_1.idf .≈ model.idf)
    @test loaded_model_1.nwords ≈ model.nwords
    @test loaded_model_1.κ == model.κ
    @test loaded_model_1.β == model.β
    @test loaded_model_1.tol ≈ model.tol
    # Model 2
    loaded_model_2 = load_lsa_model(file, Float32, sparse=false)
    @test loaded_model_2 isa LSAModel{String, Float32, Matrix{Float32}, Int}
    @test all(loaded_model_2.Uᵀ .≈ model.Uᵀ)
    @test all(loaded_model_2.Σinv .≈ model.Σinv)
    @test loaded_model_2.vocab == model.vocab
    @test loaded_model_2.vocab_hash == model.vocab_hash
    @test loaded_model_2.stats == model.stats
    @test all(loaded_model_2.idf .≈ model.idf)
    @test loaded_model_2.nwords ≈ model.nwords
    @test loaded_model_2.κ == model.κ
    @test loaded_model_2.β == model.β
    @test loaded_model_2.tol ≈ model.tol
    rm(file)
end


text:
"""
    AdditionOperation <: AbstractGroupOperation

Group operation that consists of simple addition.
"""
struct AdditionOperation <: AbstractGroupOperation end

Base.:+(e::Identity{AdditionOperation}) = e
Base.:+(e::Identity{AdditionOperation}, ::Identity{AdditionOperation}) = e
Base.:+(::Identity{AdditionOperation}, p) = p
Base.:+(p, ::Identity{AdditionOperation}) = p

Base.:-(e::Identity{AdditionOperation}) = e
Base.:-(e::Identity{AdditionOperation}, ::Identity{AdditionOperation}) = e
Base.:-(::Identity{AdditionOperation}, p) = -p
Base.:-(p, ::Identity{AdditionOperation}) = p

Base.:*(e::Identity{AdditionOperation}, p) = e
Base.:*(p, e::Identity{AdditionOperation}) = e
Base.:*(e::Identity{AdditionOperation}, ::Identity{AdditionOperation}) = e

const AdditionGroupTrait = TraitList{<:IsGroupManifold{AdditionOperation}}

adjoint_action(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, X) = X

function adjoint_action!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Y, p, X)
    return copyto!(G, Y, p, X)
end

identity_element(::AdditionGroupTrait, G::AbstractDecoratorManifold, p::Number) = zero(p)

function identity_element!(::AdditionGroupTrait, G::AbstractDecoratorManifold, p) where {𝔽}
    return fill!(p, zero(eltype(p)))
end

Base.inv(::AdditionGroupTrait, G::AbstractDecoratorManifold, p) = -p
Base.inv(::AdditionGroupTrait, G::AbstractDecoratorManifold, e::Identity) = e

inv!(::AdditionGroupTrait, G::AbstractDecoratorManifold, q, p) = copyto!(G, q, -p)
function inv!(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    q,
    ::Identity{AdditionOperation},
)
    return identity_element!(G, q)
end
function inv!(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    q::Identity{AdditionOperation},
    e::Identity{AdditionOperation},
)
    return q
end

function is_identity(::AdditionGroupTrait, G::AbstractDecoratorManifold, q; kwargs...)
    return isapprox(G, q, zero(q); kwargs...)
end
# resolve ambiguities
function is_identity(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    ::Identity{AdditionOperation};
    kwargs...,
)
    return true
end
function is_identity(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    ::Identity;
    kwargs...,
)
    return false
end

compose(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, q) = p + q

function compose!(::AdditionGroupTrait, G::AbstractDecoratorManifold, x, p, q)
    x .= p .+ q
    return x
end

function translate_diff(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    p,
    q,
    X,
    ::ActionDirection,
)
    return X
end

function translate_diff!(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    Y,
    p,
    q,
    X,
    ::ActionDirection,
)
    return copyto!(G, Y, p, X)
end

function inverse_translate_diff(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    p,
    q,
    X,
    ::ActionDirection,
)
    return X
end

function inverse_translate_diff!(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    Y,
    p,
    q,
    X,
    ::ActionDirection,
)
    return copyto!(G, Y, p, X)
end

exp_lie(::AdditionGroupTrait, G::AbstractDecoratorManifold, X) = X

exp_lie!(::AdditionGroupTrait, G::AbstractDecoratorManifold, q, X) = copyto!(G, q, X)

log_lie(::AdditionGroupTrait, G::AbstractDecoratorManifold, q) = q
function log_lie(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    ::Identity{AdditionOperation},
)
    return zero_vector(G, identity_element(G))
end
log_lie!(::AdditionGroupTrait, G::AbstractDecoratorManifold, X, q) = copyto!(G, X, q)
function log_lie!(
    ::AdditionGroupTrait,
    G::AbstractDecoratorManifold,
    X,
    ::Identity{AdditionOperation},
)
    return zero_vector!(G, X, identity_element(G))
end

lie_bracket(::AdditionGroupTrait, G::AbstractDecoratorManifold, X, Y) = zero(X)

lie_bracket!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Z, X, Y) = fill!(Z, 0)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "gmsh"
version = v"4.8.4"

# Collection of sources required to build Gmsh
sources = [
    ArchiveSource("https://gmsh.info/src/gmsh-$version-source.tgz", "760dbdc072eaa3c82d066c5ba3b06eacdd3304eb2a97373fe4ada9509f0b6ace")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gmsh-4.8.4-source
install_license LICENSE.txt
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      -DENABLE_BUILD_DYNAMIC=1 \
      ..
make -j${nproc}
make install
mv ${prefix}/lib/gmsh.jl ${prefix}/lib/gmsh.jl.bak
sed ${prefix}/lib/gmsh.jl.bak \
  -e 's/^\(import Libdl\)/#\1/g' \
  -e 's/^\(const lib.*\)/#\1/g' \
  -e 's/^\(module gmsh\)$/\1\nusing gmsh_jll: libgmsh\nconst lib = libgmsh/g' \
  > ${prefix}/lib/gmsh.jl
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgmsh", "gmsh"], :libgmsh),
    FileProduct("lib/gmsh.jl",:gmsh_api)
]

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

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,  preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Geometry which displays points at given (x, y) positions.

struct PointGeometry <: Gadfly.GeometryElement
    tag::Symbol
end
PointGeometry(; tag=empty_tag) = PointGeometry(tag)

"""
    Geom.point

Draw scatter plots of the `x` and `y` aesthetics.

# Optional Aesthetics
- `color`: Categorical data will choose maximally distinguishable colors from
  the LCHab color space.  Continuous data will map onto LCHab as well.  Colors
  can also be specified explicitly for each data point with a vector of colors of
  length(x).  A vector of length one specifies the color to use for all points.
  Default is `Theme.default_color`.
- `shape`: Categorical data will cycle through `Theme.point_shapes`.  Shapes
  can also be specified explicitly for each data point with a vector of shapes of
  length(x).  A vector of length one specifies the shape to use for all points.
  Default is `Theme.point_shapes[1]`.
- `size`: Categorical data and vectors of `Ints` will interpolate between
  `Theme.point_size_{min,max}`.  A continuous vector of `AbstractFloats` or
  `Measures` of length(x) specifies the size of each data point explicitly.  A
  vector of length one specifies the size to use for all points.  Default is
  `Theme.point_size`.
- `alpha`: Categorical data will use the alpha palette in `Theme.alphas`.
  Continuous data will remap from 0-1. Default is `Theme.alphas[1]`.
"""
const point = PointGeometry

element_aesthetics(::PointGeometry) = [:x, :y, :size, :color, :shape, :alpha]

# Generate a form for a point geometry.
#
# Args:
#   geom: point geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::PointGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.point", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.point", aes, :x, :y)

    default_aes = Gadfly.Aesthetics()
    default_aes.shape = Function[Shape.circle]
    default_aes.color = discretize_make_ia(RGBA{Float32}[theme.default_color])
    default_aes.size = Measure[theme.point_size]
    default_aes.alpha = [theme.alphas[1]]
    aes = inherit(aes, default_aes)

    if eltype(aes.size) <: Int
      size_min, size_max = extrema(aes.size)
      size_range = size_max - size_min
      point_size_range = theme.point_size_max - theme.point_size_min
      interpolate_size(x) = theme.point_size_min + (x-size_min) / size_range * point_size_range
    end

    aes_alpha = eltype(aes.alpha) <: Int ? theme.alphas[aes.alpha] : aes.alpha

    ctx = context()

    for (x, y, color, size, shape, alpha) in Compose.cyclezip(aes.x, aes.y, aes.color, aes.size, aes.shape, aes_alpha)
        shapefun = typeof(shape) <: Function ? shape : theme.point_shapes[shape]
        sizeval = typeof(size) <: Int ? interpolate_size(size) : size
        strokecolor = aes.color_key_continuous != nothing && aes.color_key_continuous ?
                    theme.continuous_highlight_color(color) :
                    theme.discrete_highlight_color(color)
        class = svg_color_class_from_label(aes.color_label([color])[1])
        compose!(ctx, (context(),
              (context(), shapefun([x], [y], [sizeval]), svgclass("marker")),
              fill(color), stroke(strokecolor), fillopacity(alpha),
              svgclass(class)))
    end

    compose!(ctx, linewidth(theme.highlight_width))

    return compose!(context(order=4), svgclass("geometry"), ctx)
end


text:
#### Layer dimensions bug

cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)

## Bug as noticed in v0.6.0
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.6.0")
using SimpleSDMLayers

t1 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...) # 331 x 571, as in earlier versions
t2 = SimpleSDMPredictor(WorldClim, BioClim, 1)[coords] # 332 x 571, one extra line!!
isequal(t1, t2) # not equal

## Verification in v0.4.8
# Version used in betadiversity-hotspots before jump to v0.6.0
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.4.8")
using SimpleSDMLayers

t1 = worldclim(1; coords...).grid # 331 x 571
t2 = worldclim(1)[coords].grid # 332 x 571!!
isequal(t1, t2) # Nope

## Verification in v0.3.6
# Version used in betadiversity-hotspots before jump to v0.4.8
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.3.6")
using SimpleSDMLayers

t1 = worldclim(1; coords...) # 331 x 571
t2 = worldclim(1)[coords] # 331 x 571
isequal(size(t1), size(t2))

## Verification in v0.4.0
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.4.0")
using SimpleSDMLayers

t1 = worldclim(1; coords...)
t2 = worldclim(1)[coords]
size(t1) # 331 x 571
size(t2) # 331 x 571
isequal(t1, t2)
# Still works...

## Verification in v0.4.4
# Trying a middle version to find out
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.4.4")
using SimpleSDMLayers

t1 = worldclim(1; coords...)
t2 = worldclim(1)[coords]
isequal(size(t1), size(t2))
# Doesn't work

## Verification in v0.4.2
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.4.2")
using SimpleSDMLayers

t1 = worldclim(1; coords...)
t2 = worldclim(1)[coords]
isequal(size(t1), size(t2)) # Equal!

isequal(t1, t2) # Nope
# Doesn't work

isequal(t1.grid, t2.grid) # grids are equal
lims1 = t1.left, t1.right, t1.bottom, t1.top
lims2 = t2.left, t2.right, t2.bottom, t2.top
isequal(lims1, lims2) # equal too
# Grids are equal, limits are equal too...
# What's going on??

## Verification in v0.4.3
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.4.3")
using SimpleSDMLayers

t1 = worldclim(1; coords...)
t2 = worldclim(1)[coords]
isequal(size(t1), size(t2)) # nope
# So dimensions problem occurred between v0.4.2 and v0.4.3
# But what about equality problem?

## Re-Verification in v0.3.6
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.3.6")
using SimpleSDMLayers

t1 = worldclim(1; coords...)
t2 = worldclim(1)[coords]
isequal(size(t1), size(t2))
isequal(t1, t2)
# Ok some equality problem was still around then
# Let's forget about it for now and focus on the dimensions problem
# So back to 0.4.2 and 0.4.3

# Back on v0.4.3
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.add(name = "SimpleSDMLayers", version = "v0.4.3")
using SimpleSDMLayers

t1 = worldclim(1; coords...)
t2 = worldclim(1)[coords]
isequal(size(t1), size(t2))

t1.grid
t2.grid

t1.grid[end, end]
t2.grid[end, end]

t1.grid[331, 571]
t1[coords.right, coords.top]
t2[coords.right, coords.top]

layer = worldclim(1)
lat_stride = stride(layer, 2)

# _match_longitude
last(findmin(abs.(lon .- longitudes(layer))))

i = SimpleSDMLayers._match_latitude(layer, coords.top)
j = SimpleSDMLayers._match_longitude(layer, coords.right)
layer.grid[i, j]
layer.grid[i-1, j]

last(findmin(abs.(coords.top .- latitudes(layer))))
coords.top .- latitudes(layer)
collect(coords.top .- latitudes(layer))[i]
collect(coords.top .- latitudes(layer))[i-1]
coords.top .- latitudes(layer) |>
    x -> filter(y -> -0.5 < y < 0.5, x)
coords.top .- 0.5*lat_stride .- latitudes(layer) |>
    x -> filter(y -> -0.5 < y < 0.5, x)

i2 = last(findmin(abs.(coords.top .-0.5*lat_stride .- latitudes(layer))))
layer.grid[i2, j]
t1.grid[end, end]

last(findmin(abs.(coords.bottom .- 0.5*lat_stride .- latitudes(layer))))
last(findmin(abs.(coords.bottom .- latitudes(layer))))
coords.bottom .- latitudes(layer) |>
    x -> filter(y -> -0.5 < y < 0.5, x)
coords.bottom .- 0.5*lat_stride .- latitudes(layer) |>
    x -> filter(y -> -0.5 < y < 0.5, x)

last(findmin(abs.(coords.right .- 0.5*lat_stride .- longitudes(layer))))
last(findmin(abs.(coords.right .- longitudes(layer))))
last(findmin(abs.(coords.right .- longitudes(layer))))

## Dev some changes
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers

t1 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...) # 331 x 571, as in earlier versions
t2 = SimpleSDMPredictor(WorldClim, BioClim, 1)[coords] # 332 x 574, one extra line!!
isequal(t1, t2) # not equal

using Test
S = SimpleSDMResponse(rand(Bool, 4, 4), 0.3, 0.7, 0.5, 0.9)
for (i,l) in enumerate(latitudes(S))
    @test SimpleSDMLayers._match_latitude(S, l) == i
end
[SimpleSDMLayers._match_latitude(S, l) for l in latitudes(S)]

for (i,l) in enumerate(longitudes(S))
    @test SimpleSDMLayers._match_longitude(S, l) == i
end
[SimpleSDMLayers._match_latitude(S, l) for l in latitudes(S)]


text:
module Monotonicers

using Dates
using DataFrames
using Statistics

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils
import AutoMLPipeline.AbsTypes: fit!, transform!

export fit!,transform!
export Monotonicer,ismonotonic,dailyflips

# Transforms instances with nominal features into one-hot form
# and coerces the instance matrix to be of element type Float64.

"""
    Monotonicer()

Monotonic filter to detect and normalize two types of dataset:
- daily monotonic
- entirely non-decreasing/non-increasing data

Example:

    fname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
    csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
    valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
    valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
    stfier = Statifier(Dict(:processmissing=>true))
    mono = Monotonicer(Dict())

    mypipeline = @pipeline csvfilter |> valgator |> mono |> stfier
    result = fit_transform!(mypipeline)

Implements: `fit!`, `transform!`

"""
mutable struct Monotonicer <: Transformer
  model
  args

  function Monotonicer(args=Dict())
    default_args = Dict(
    )
    new(nothing, mergedict(default_args, args))
  end
end

"""
    fit!(st::Monotonicer,features::T, labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}

A function that checks if `features` are two-column data of  Dates and Values
"""
function fit!(st::Monotonicer, features::DataFrame, labels::Vector=[])
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  st.model = st.args
end

"""
    transform!(st::Monotonicer, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Normalize monotonic or daily monotonic data by taking the diffs and counting the flips.
"""
function transform!(st::Monotonicer, features::DataFrame)
  features != DataFrame() || return DataFrame()
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  sum(names(features) .== (:Date,:Value))  == 2 || error("wrong column names")
  mfeatures=features
  # double check monotonic
  # based on flips and daily flips
  if ismonotonic(features.Value)
    mfeatures = antimonotonize(features)
  end
  # daily mono reset everynight at most 2x
  ndailyflips = dailyflips(mfeatures)
  if ndailyflips < 0.5
    antimono = antimonotonize(mfeatures)
  elseif 0.5 <= ndailyflips <= 2.0
    antimono = antimonotonizedaily(mfeatures)
  else
    antimono = mfeatures
  end
  return antimono
end

function antimonotonizedaily(dat::DataFrame)
  # get diffs, negative values represent reset
  # replace negative values by succeeding neighbor
  df = deepcopy(dat)
  df.Value .= [ -1.0; diff(dat.Value) ]
  valuelength = length(df.Value)
  # find locations of negatives
  negatives = findall(x->x<0,df.Value)
  # if last element is negative, zero it
  if negatives[end] == valuelength
    df.Value[negative[end]] = 0.0
    pop!(negatives)
  end
  # replace negatives with immediate neighbor
  for val in reverse(negatives)
    df.Value[val] = df.Value[val+1]
  end
  return df
end

function antimonotonize(data::DataFrame)
  dat = deepcopy(data)
  typeof(dat) <: DataFrame || error("input must be a dataframe")
  ncol(dat) == 2 || error("input must have two columns")
  sum(names(dat) .== (:Date,:Value))  == 2 || error("wrong column names")
  vals = dat.Value
  fvals = diff(vals)
  dvals = [fvals[1];fvals]
  dat.Value = dvals
  return dat
end

function dailyflips(dat::DataFrame)
  dates = dat.Date
  values = dat.Value
  datemin = minimum(dates)
  datemax = maximum(dates)
  ndays = round(datemax-datemin,Dates.Day(1))
  nflips::Float64 = sum((skipmissing(values) |> collect |> diff) .< 0 )
  avg::Float64 = nflips/ndays.value
end

function ismonotonic(dat::Vector)
  #check if monotonic increasing
  cdata=dat |> skipmissing |> collect
  maxflips = minimum([10;div(length(cdata),2)]) # flip cutoff to identify monotonic
  nflips = length(cdata) - sum(cdata .== accumulate(max,cdata))
  nflips < maxflips ? true : false
end

end


text:
# This code is based on a routine originally copyright Chris Sims.
# See http://sims.princeton.edu/yftp/gensys/

# eye(n::Integer) deprecated in Julia v0.7.0 onwards
@inline eye(n::Integer) = Matrix{Float64}(I, n, n)

"""
```
gensys(Γ0, Γ1, c, Ψ, Π)
gensys(Γ0, Γ1, c, Ψ, Π, div)
gensys(F::LinearAlgebra.GeneralizedSchur, c, Ψ, Π)
gensys(F::LinearAlgebra.GeneralizedSchur, c, Ψ, Π, div)
```

Generate state-space solution to canonical-form DSGE model.

System given as
```
Γ0*y(t) = Γ1*y(t-1) + c + Ψ*z(t) + Π*η(t),
```
with z an exogenous variable process and η being endogenously
determined one-step-ahead expectational errors.

Returned system is
```
y(t) = G1*y(t-1) + C + impact*z(t) + ywt*inv(I-fmat*inv(L))*fwt*z(t+1)
```
Returned values are
```
G1, C, impact, fmat, fwt, ywt, gev, eu, loose
```

If `z(t)` is i.i.d., the last term drops out.

If `div` is omitted from argument list, a `div`>1 is calculated.

### Return codes

* `eu[1] = 1` for existence
* `eu[2] = 1` for uniqueness
* `eu[1] = -1` for existence only with not-s.c. z
* `eu = [-2, -2]` for coincident zeros
* `eu = [-3, -3]` if a LAPACKException is thrown while computing the Schur
  decomposition

### Notes

We constrain Julia to use the complex version of the `schurfact` routine regardless of the
types of `Γ0` and `Γ1`, to match the behavior of Matlab.  Matlab always uses the complex version
of the Schur decomposition, even if the inputs are real numbers.
"""
function gensys(Γ0::Array{Float64, 2}, Γ1::Array{Float64, 2}, c::Array{Float64, 1}, Ψ::Array{Float64, 2}, Π::Array{Float64, 2}, div::Float64 = 0.0; verbose::Symbol = :high)
    F = try
        schur!(complex(Γ0), complex(Γ1))
    catch ex
        if isa(ex, LinearAlgebra.LAPACKException)
            Base.@info "LAPACK exception thrown while computing Schur decomposition of Γ0 and Γ1."
            eu = [-3, -3]

            G1 = Array{Float64, 2}(undef,0,0)
            C = Array{Float64, 1}(undef,0)
            impact = Array{Float64, 2}(undef,0,0)
            #fmat = Array{Complex{Float64}, 2}(undef,0,0)
            #fwt = Array{Complex{Float64}, 2}(undef,0,0)
            #ywt = Vector{Complex{Float64}}(undef,0)
            #gev = Vector{Complex{Float64}}(undef,0)
            #loose = Array{Float64, 2}(undef,0,0)

            return G1, C, impact, eu #G1, C, impact, fmat, fwt, ywt, gev, eu, loose
        else
            rethrow(ex)
        end
    end
    if div == 0.0
        gensys(F, c, Ψ, Π; verbose = verbose)
    else
        gensys(F, c, Ψ, Π, div; verbose = verbose)
    end
end

function gensys(F::LinearAlgebra.GeneralizedSchur, c::Array{Float64, 1}, Ψ::Array{Float64, 2}, Π::Array{Float64, 2}; verbose::Symbol = :low)
    gensys(F, c, Ψ, Π, new_div(F), verbose = verbose)
end

# Method that does the real work. Work directly on the decomposition F
function gensys(F::LinearAlgebra.GeneralizedSchur, c::Array{Float64, 1}, Ψ::Array{Float64, 2}, Π::Array{Float64, 2}, div::Float64; verbose::Symbol = :low)
    eu      = [0, 0]
    ϵ       = 1e-6  # small number to check convergence
    nunstab = 0
    zxz     = 0
    a, b,   = F.S, F.T
    n       = size(a, 1)

    select = BitArray(undef, n)
    for i in 1:n
        # nunstab is the variable name used by Chris Sims, but it seems
        # that nunstab should actually correspond to the number of stable λs
        # i.e. nunstab += 1/div > abs(a[i,i])/abs(b[i,i]), which is basically
        # 1 - a small number > abs(a[i,i])/abs(b[i,i])
        select[i] = !(abs(b[i, i]) > div * abs(a[i, i]))
        if (abs(a[i, i]) < ϵ) && (abs(b[i, i]) < ϵ)
            zxz = 1
        end
    end
    nunstab = n - sum(select)

    if zxz == 1
        @warn "Coincident zeros. Indeterminacy and/or nonexistence."
        eu=[-2, -2]

        G1     = Array{Float64, 2}(undef,0, 0)
        C      = Array{Float64, 1}(undef,0)
        impact = Array{Float64, 2}(undef,0)
        #fmat   = Array{Complex{Float64}, 2}(undef,0,0)
        #fwt    = Array{Complex{Float64}, 2}(undef,0,0)
        #ywt    = Vector{Complex{Float64}}(undef,0)
        #gev    = Vector{Complex{Float64}}(undef,0)
        #loose  = Array{Float64, 2}(undef,0,0)

        return G1, C, impact, eu #G1, C, impact, fmat, fwt, ywt, gev, eu, loose
    end

    FS = ordschur!(F, select)
    a, b, qt, z = FS.S, FS.T, FS.Q, FS.Z
    gev   = hcat(diag(a), diag(b))
    qt1   = qt[:, 1:(n - nunstab)]
    qt2   = qt[:, (n - nunstab + 1):n]
    etawt = adjoint(qt2) * Π
    neta  = size(Π, 2)

    # branch below is to handle case of no stable roots, rather than quitting with an error
    # in that case.
    if nunstab == 0
        etawt = zeros(0, neta)
        ueta  = zeros(0, 0)
        deta  = zeros(0, 0)
        veta  = zeros(neta, 0)
        bigev = 0
    else
        etawtsvd = svd!(etawt)
        bigev = (LinearIndices(etawtsvd.S))[findall(etawtsvd.S .> ϵ)]
        ueta  = etawtsvd.U[:, bigev]
        veta  = etawtsvd.V[:, bigev]
        deta  = Matrix(Diagonal(etawtsvd.S[bigev]))
    end

    existence = length(bigev) >= nunstab
    if existence
        eu[1] = 1
    else
        if VERBOSITY[verbose] >= VERBOSITY[:high]
            @warn "Nonexistence: number of unstable roots exceeds number of jump variables"
        end
    end

    # Note that existence and uniqueness are not just matters of comparing
    # numbers of roots and numbers of endogenous errors.  These counts are
    # reported below because usually they point to the source of the problem.

    # branch below to handle case of no stable roots
    if nunstab == n
        etawt1 = zeros(0, neta)
        bigev  = 0
        ueta1  = zeros(0, 0)
        veta1  = zeros(neta, 0)
        deta1  = zeros(0, 0)
    else
        etawt1 = adjoint(qt1) * Π
        ndeta1 = min(n - nunstab, neta)
        etawt1svd = svd!(etawt1)
        bigev = (LinearIndices(etawt1svd.S))[findall(etawt1svd.S .> ϵ)]
        ueta1 = etawt1svd.U[:, bigev]
        veta1 = etawt1svd.V[:, bigev]
        deta1 = Matrix(Diagonal(etawt1svd.S[bigev]))
    end
    if isempty(veta1)
        unique = true
    else
        loose = veta1 - (veta * adjoint(veta)) * veta1
        loosesvd = svd!(loose)
        nloose = sum(abs.(loosesvd.S) .> ϵ * n)
        unique = (nloose == 0)
    end

    if unique
        eu[2] = 1
    else
        if VERBOSITY[verbose] >= VERBOSITY[:high]
            @warn "Indeterminacy: $(nloose) loose endogeneous error(s)"
        end
    end

    tmat = hcat(eye(n - nunstab), -(ueta * (deta \ veta') * veta1 * (deta1 * adjoint(ueta1)))')
    G0   = vcat(tmat * a, hcat(zeros(nunstab, n - nunstab), eye(nunstab)))
    G1   = vcat(tmat * b, zeros(nunstab, n))

    # G0 is always non-singular because by construction there are no zeros on
    # the diagonal of a(1:n-nunstab,1:n-nunstab), which forms G0's ul corner.
    G0I    = inv(G0)
    G1     = G0I * G1
    usix   = (n - nunstab + 1):n
    Busix  = b[usix,usix]
    Ausix  = a[usix,usix]
    C      = G0I * vcat(tmat * (adjoint(qt) * c), (Ausix - Busix) \ (adjoint(qt2) * c))
    impact = G0I * vcat(tmat * (adjoint(qt) * Ψ), zeros(nunstab, size(Ψ, 2)))
    #fmat   = Busix \ Ausix
    #fwt    = -Busix \ (adjoint(qt2) * Ψ)
    #ywt    = G0I[:, usix]

    #loose  = G0I * vcat(etawt1 * (eye(neta) - (veta * adjoint(veta))), zeros(nunstab, neta))

    G1     = real(z * (G1 * adjoint(z)))
    C      = real(z * C)
    impact = real(z * impact)
    #loose  = real(z * loose)

    #ywt = z * ywt

    return G1, C, impact, eu #G1, C, impact, fmat, fwt, ywt, gev, eu, loose
end

function new_div(F::LinearAlgebra.GeneralizedSchur)
    ϵ    = 1e-6  # small number to check convergence
    n    = size(F.T, 1)
    a, b = F.S, F.T
    div  = 1.01
    for i in 1:n
        if abs(a[i, i]) > 0
            divhat = abs(b[i, i]) / abs(a[i, i])
            if 1 + ϵ < divhat && divhat <= div
                div = .5 * (1.0 + divhat)
            end
        end
    end
    return div
end


text:
import REPL
using REPL.TerminalMenus

println("No juliab.so found, follow instructions to create")

categories_to_packages = Dict(
    "Plotting" => ["Plots"],
    "Data" => ["CSV", "DataFrames", "DataFramesMeta","PrettyTables"],
    "Utils" => ["Revise", "OhMyREPL"],
    "StatsKit" => ["StatsKit"],
    "Arrays" => ["OffsetArrays", "StaticArrays"],
    "Maths" => ["Tau", "LinearAlgebra"],
    "Pluto" => ["Pluto"],
)
all_packages = unique(mapreduce(last, vcat, categories_to_packages))

# Show menu with the categories
all_category_names = collect(keys(categories_to_packages))
category_menu = MultiSelectMenu(all_category_names)
chosen_category_indices = request("Select the categories you care about:", category_menu)
chosen_category_names = all_category_names[collect(chosen_category_indices)]

# Collect all the packages from selected categories
chosen_category_packages = unique(mapreduce(c->categories_to_packages[c], vcat, chosen_category_names))

# Preselect packages from the categories
category_menu = MultiSelectMenu(all_packages; pagesize=30)
for (i, package) in enumerate(category_menu.options)
    if package in chosen_category_packages
        push!(category_menu.selected, i)
    end
end

# Show menu for packages
chosen_package_indices = request("For finer grain control, select individual packages you care about, or just continue:", category_menu)
chosen_package_names = all_packages[collect(chosen_package_indices)]

package_precompile_execution_files=["package_precompile_execution_files/$(p).jl" for p in chosen_package_names]
package_precompile_statements_files=["package_precompile_statements_files/$(p).jl" for p in chosen_package_names]
category_precompile_execution_files=["category_precompile_execution_files/$(p).jl" for p in chosen_category_names]

precompile_execution_files = vcat(package_precompile_execution_files, category_precompile_execution_files)
filter!(isfile, precompile_execution_files)

filter!(isfile, package_precompile_statements_files)

if isempty(chosen_package_names)
    error("No packages specified, no point building .so file. Just use julia instead")
end

import Pkg
println("Adding packages: ", join(chosen_package_names, ", "))
Pkg.add(chosen_package_names)
Pkg.update()

Pkg.add("PackageCompiler")
using PackageCompiler
println("Compiling")
chosen_packages = Symbol.(chosen_package_names)
create_sysimage(chosen_packages, sysimage_path=joinpath(@__DIR__(), "juliab.so"),
    precompile_statements_file=package_precompile_statements_files,
    precompile_execution_file=precompile_execution_files)

max_stars_repo_name: jondea/juliab

max_issues_repo_name: jondea/juliab

max_forks_repo_name: jondea/juliab


text:
need_variablify(p::Sets.Projection) = need_variablify(p.set)
Sets.perspective_variable(p::Sets.Projection) = Sets.perspective_variable(p.set)
clear_spaces(p::Sets.Projection) = clear_spaces(p.set)

need_variablify(lm::Sets.LinearImage) = need_variablify(lm.set)
function variablify(lm::Sets.LinearImage)
    return Sets.LinearImage(variablify(lm.set), lm.A)
end

"""
    apply_map(li::Sets.LinearImage{<:Sets.Polar}, new_vars)

The set ``(AS)^\\circ``, the polar of the set ``AS``, is ``A^{-\\top}S^\\circ``.
"""
function apply_map(li::Sets.LinearImage{<:Sets.Polar}, new_vars)
    return Sets.polar(apply_map(Sets.LinearPreImage(Sets.polar(li.set), li.A'), new_vars))
end

function apply_map(li::Sets.LinearPreImage{<:Sets.PolarPoint}, new_vars)
    return Sets.PolarPoint(li.A' * li.set.a)
end

"""
    apply_map(li::Sets.LinearPreImage{<:Sets.Ellipsoid}, new_vars)

Given ``S = \\{\\, x \\mid x^\\top Q x \\le 1\\,\\}``, we have
``A^{-1}S = \\{\\, x \\mid x^\\top A^\\top Q A x \\le 1\\,\\}``.
"""
function apply_map(li::Sets.LinearPreImage{<:Sets.Ellipsoid}, new_vars)
    return Sets.Ellipsoid(Symmetric(li.A' * li.set.Q * li.A))
end

function apply_map(li::Sets.LinearPreImage{<:Sets.Piecewise}, new_vars)
    return Sets.Piecewise(
        [apply_map(Sets.LinearPreImage(set, li.A), new_vars) for set in li.set.sets],
        li.A \ li.set.polytope,
        [li.A \ piece for piece in li.set.pieces],
        li.set.graph # /!\ FIXME The nij are now incorrect
    )
    return Sets.Ellipsoid(Symmetric(li.A' * li.set.Q * li.A))
end

"""
    apply_map(li::Sets.LinearPreImage{<:Sets.PolySet}, new_vars)

Given ``S = \\{\\, x \\mid p(x) \\le 1\\,\\}``, we have
``A^{-1}S = \\{\\, x \\mid p(Ax) \\le 1\\,\\}``.
"""
function apply_map(li::Sets.LinearPreImage{<:Sets.PolySet}, new_vars)
    deg = li.set.degree
    @assert iseven(deg)
    q = apply_matrix(li.set.p, li.A, new_vars, div(deg, 2))
    return Sets.PolySet(deg, q)
end
function apply_map(li::Sets.LinearPreImage{<:Sets.ConvexPolySet}, new_vars)
    deg = li.set.degree
    @assert iseven(deg)
    q = apply_matrix(li.set.p, li.A, new_vars, div(deg, 2))
    return Sets.ConvexPolySet(deg, q, nothing)
end

"""
    apply_map(li::Sets.LinearImage{<:Sets.PerspectiveDualOf{<:Union{Sets.PerspectiveEllipsoid,
                                                                    Sets.PerspectiveConvexPolynomialSet}}}, new_vars)

The set ``(AS)^\\circ``, the polar of the set ``AS``, is ``A^{-\\top}S^\\circ``
and given ..., we have
...
"""
function apply_map(li::Sets.LinearImage{<:Sets.PerspectiveDualOf{<:Sets.Householder{T}}}, new_vars) where T
    old_vars = Sets.space_variables(li.set)
    p = MP.subs(Sets.perspective_gauge0(li.set.set), old_vars => li.A' * new_vars)
    dual = Sets.Householder(Sets.UnknownSet{T}(), p, li.A * li.set.set.h,
                            Sets.perspective_variable(li.set), new_vars)
    return Sets.perspective_dual(dual)
end

abstract type SymbolicVariable end

# FIXME, for Sets.AbstractSet, we should apply it directly
function Base.:(*)(A::AbstractMatrix, set::Union{SetVariableRef, Sets.AbstractSet, SymbolicVariable})
    return Sets.LinearImage(set, A)
end


text:
@testset "releg" begin

    @testset "fuse/split" begin
        space = (
            U1Space(0=>2, 1=>3),
            U1Space(0=>1, 1=>1),
            dual(U1Space(0=>3, 1=>4)),
        )

        # Note that all the possible sectors are given!
        sects = [Sector{U1}(0, 0, 0),
                 Sector{U1}(0, 1, 1),
                 Sector{U1}(1, 0, 1)]
        nzblks = [reshape(collect(1:6), 2,1,3),
                  reshape(collect(7:14), 2,1,4),
                  reshape(collect(15:26), 3,1,4)]

        A0 = SymTensor(zero(U1), space, SortedDict(zip(sects,nzblks)))
        #A0 = SymTensor(zero(U1), space,
        #SortedDict{Sector{U1,3}, Array{Int,3}}(zip(sects,nzblks)))

        A1 = fuselegs(A0, 1, 2)
        A2 = splitleg(A1, 1, space[1:2])
        @test A0 == A2

        A1 = fuselegs(A0, 2, 2, true)
        A2 = splitleg(A1, 2, space[2:3])
        @test A0 == A2
    end

    @testset "dual" begin
        A = eye(U1Space(0=>1, 1=>2, 2=>1))
        A[Sector{U1}(2, 2)] *= -1
        A[Sector{U1}(1, 1)] *= 1/sqrt(2) .* [-1. 1; 1 1]

        B = splitleg(A, 1, U1Space(0=>1, 1=>1), U1Space(0=>1, 1=>1))
        B = splitleg(A, 1, U1Space(0=>1, 1=>1), U1Space(0=>1, 1=>1))
        fuselegs(dual(B), 1, 2, false) == dual(A)
    end

    @testset "U1 unitary" begin
        A = eye(Float64, U1Space(0=>1, 1=>2, 2=>1))
        A[Sector{U1}(1, 1)] =  [2. 3;4 5]
        d = U1Space(0=>1, 1=>1)
        cod = (d, d)
        dom = (dual(d), dual(d))
        B = splitleg(splitleg(A, 2, dom), 1, cod)
        @test sectors(B) == sort([Sector{U1}(1,1,1,1), Sector{U1}(1,0,0,1),
                                  Sector{U1}(0,1,0,1), Sector{U1}(1,0,1,0),
                                  Sector{U1}(0,1,1,0), Sector{U1}(0,0,0,0)])
        i = ones(1,1,1,1)
        @test collect(values(B.blocks)) == [i,2*i, 4*i, 3*i,5*i, 1*i]
    end

    @testset "Associativity" begin
        d = U1Space(-1=>1, 0=>2, 1=>1)
        space = (d, dual(d), d, dual(d))

        A = fill_linearindex(zero(U1),  space[1:4])
        A1 = fuselegs(A, 2, 3)
        A2 = fuselegs(A, 3, 2)
        A3 = fuselegs(A, 1, 3)
        A4 = fuselegs(A, 1, 2)
        A5 = fuselegs(A, 2 ,2)

        #@test A2 ≈ splitleg(A1, 2, space[2], fuse(space[3], space[4]))
        #@test A4 ≈ splitleg(A3, 1, fuse(space[1], space[2]), space[3])
        #@test A1 ≈ fuselegs(A2, 2, 2)
        #@test A1 ≈ fuselegs(A5, 2, 2)
        #@test A3 ≈ fuselegs(A4, 1, 2)
    end
end

@testset "contract" begin

    @testset "converts" begin
        vlist = (
            U1Space(0=>2, 1=>3, 2=>4),
            U1Space(-2=>4, -1=>3, 0=>2, 1=>1),
            U1Space(-1=>1, 0=>2, 1=>3, 2=>4),
            U1Space(-2=>3, -1=>2, 0=>1),
        )
        A = fill_linearindex(vlist)
        @test SymMatrix(A, [1,2], [3,4]) == fuselegs(fuselegs(A, 3, 2, true), 1, 2)
        @test SymMatrix(A, [1,2], [3,4]) == fuselegs(fuselegs(A, 1, 2), 2, 2, true)
        #@test SymVector(A, [1,2,3,4]) == fuselegs(A,1,4)

        perm = randperm(4)
        @test SymMatrix(A, perm[1:2], perm[3:4]) ==
            fuselegs(fuselegs(permutelegs(A, perm), 3, 2, true), 1, 2)
        #@test SymVector(A, perm) == fuselegs(permutelegs(A, perm), 1,  4)
    end

    @testset "matrices" begin
        vlist = (
            U1Space(0=>2, 1=>3, 2=>4),
            dual(U1Space(0=>4, 1=>2, 2=>2, 3=>3)),
            U1Space(0=>4, 1=>2, 2=>2, 3=>3),
            dual(U1Space(0=>3, 1=>5, 2=>2)),
        )
        A = fill_linearindex(vlist[1:2])
        B = fill(1.0, zero(U1), vlist[3:4])

        A_ = array(A)
        B_ = array(B)

        C = contract(A, (1, -1), B, (-1, 2))

        @test A_*B_ == array(C)
    end

    @testset "multi leg tensors" begin
        vlist = (
            U1Space(0=>2, 1=>2, 2=>2),
            U1Space(0=>2, -1=>2, -2=>2, -3=>2),
            U1Space(0=>2, 1=>2, 2=>2),
            dual(U1Space(0=>2, 1=>2)),
            U1Space(0=>2, 1=>2),
            U1Space(0=>2, 1=>2, 2=>2, 3=>2, 4=>2),
            U1Space(0=>2, 1=>2),
            U1Space(0=>2, -1=>2, -2=>2)
        )
        A = rand(Float64, zero(U1), vlist[1:4])
        B = rand(Float64, zero(U1), vlist[5:8])

        @test splitleg(fuselegs(A, 1, 3), 1, vlist[1:3]) ≈ A
        @test splitleg(fuselegs(A, 2, 2), 2, vlist[2:3]) ≈ A
        @test splitleg(fuselegs(A, 2, 3), 2, vlist[2:4]) ≈ A

        A_ = array(A)
        B_ = array(B)

        C = contract(A, (1,2,3,-1), B, (-1,4,5,6))
        @tensor C_[a,b,c,e,f,g] := A_[a,b,c,d] * B_[d,e,f,g]

        @test C_ ≈ array(C);
    end
end

@testset "LA fns" begin
    vlist = (
        U1Space(0=>2, 1=>3, 2=>4),
        U1Space(0=>4, -1=>2, -2=>2, -3=>3),
        U1Space(0=>3, 1=>5, 2=>2),
        U1Space(0=>6, -1=>7),
        U1Space(0=>6, 1=>7),
        U1Space(0=>2, 1=>3, 2=>4, 3=>5, 4=>6),
        U1Space(0=>2, 1=>3),
        U1Space(0=>3, -1=>4, -2=>5)
    )
    @testset "dot" begin
        A = rand(ComplexF64, zero(U1), vlist[1:4])
        A_ = array(A)
        @test dot(A_,A_) ≈ dot(A, A)
    end
end

@testset "factorization" begin

    V1 = U1Space(0=>2, 1=>3, 2=>4)
    V2 = U1Space(0=>4, -1=>2, -2=>2, -3=>3)
    V3 = U1Space(0=>4, 1=>2, 2=>2, 3=>3)
    V4 = U1Space(0=>3, -1=>5, -2=>2)

    @testset "svd" begin
        A1 = rand(Float64, zero(U1), (V1,V2))
        A2 = rand(Float64, U1(1), (V1,V2))
        A3 = rand(Float64, U1(-1), (V1,V2))
        u,s,v = _svd_(A1)
        @test u*s*v ≈ A1
        u,s,v = _svd_(A2)
        @test u*s*v ≈ A2
        u,s,v = _svd_(A3)
        @test u*s*v ≈ A3

        B1 = rand(ComplexF64, zero(U1), (V1,dual(V3)))
        B2 = rand(ComplexF64, U1(1), (V1,dual(V3)))
        B3 = rand(ComplexF64, U1(-1), (V1,dual(V3)))
        u,s,v = _svd_(B1)
        @test u*s*v ≈ B1
        u,s,v = _svd_(B2)
        @test u*s*v ≈ B2
        u,s,v = _svd_(B3)
        @test u*s*v ≈ B3
    end

    @testset "svdtrunc" begin
        A1 = SymTensor(rand, ComplexF64, zero(U1), (V1,V2))
        @test _svd_(A1) ≈ svdtrunc(A1)
    end
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
# Use the access token to access the Spotify Web API
# https://developer.spotify.com/documentation/general/guides/authorization-guide/#client-credentials-flow

"""
    spotify_request(url_ext::String, method::String= "GET"; scope = "client-credentials")
     -> (r::JSON3.Object, retry_in_seconds::Int)

Access the Spotify Web API.
Error results return an empty Object.
Errors are written to 'stderr', expect for 'API rate limit exceeded', as
the output would typically occur in the middle of recursive calls.
"""
function spotify_request(url_ext::String, method::String= "GET"; scope = "client-credentials")
    # Not so cool, but this is how we get rid of spaces in vectors of strings:
    url_ext = replace(url_ext, ' ' => "")
    url = "https://api.spotify.com/v1/$url_ext"
    headers = header_maker(;scope)
    isnothing(headers) && return()
    resp = HTTP.Messages.Response()
    try
        resp = HTTP.request(method, url, headers)
    catch e
        if  e isa HTTP.ExceptionRequest.StatusError && e.status ∈ [400, 401, 403, 404, 429]
            response_body = e.response.body |> String
            if e.status == 429 # API rate limit temporarily exceeded.
                retry_in_seconds =  HTTP.header(e.response, "retry-after")
                return JSON3.Object(), parse(Int, retry_in_seconds)

            elseif e.status == 404 # Not found, e.g. when a track/playlist/album ID is incorrect
                @info "404 Not Found - Check if input arguments are okay"
                return JSON3.Object(), 0

            elseif e.status == 403
                @info "403 - Bad OAuth request (wrong consumer key, bad nonce, expired timestamp...)"
                return JSON3.Object(), 0

            elseif e.status == 400 # e.g. when a search query is empty
                @info "Check if the search query is present"
                return JSON3.Object(), 0

            else
                @info "You should try to re-authenticate the user"
                @error response_body
                return JSON3.Object(), 0
            end
        else
            @warn "HTTP.request call: method = $method\n headers = $headers \n $url_ext"
            @error e
            @error response_body
            return JSON3.Object(), 0
        end
    end
    response_body = resp.body |> String
    response_body |> JSON3.read , 0
end

#=
"""
Access the Spotify web api through the 'Implicit grant flow'.
Implicit grant flow is for clients that are implemented entirely
using JavaScript and running in the resource owner’s browser.
Julia: No prob.!
..but we do launch a browser for login..
"""
function other_request()
    error("ouch - not used")
    url = OAUTH_AUTHORIZE_URL
    headers = header_maker("Bearer")
    resp = HTTP.Messages.Response()
    method = "GET"
    try
        resp = HTTP.request(method, url, headers)
    catch e
        request = "HTTP.request call: method = $method\n  url = $url\n  headers = $headers \n"
        @error request
        @error e
        return nothing
    end
    response_body = resp.body |> String
    response_body |> JSON3.read
end
=#


text:
export
    LaneChangeChoice,
    LaneChangeModel,
    get_lane_offset,

    DIR_RIGHT,
    DIR_MIDDLE,
    DIR_LEFT,

    MOBIL,
    TimLaneChanger

include("lane_change_models.jl")
include("MOBIL.jl")
include("tim_lane_changers.jl")


text:
"""
$(TYPEDEF)

Adjacency struct. Essentially, this is the sparsity pattern of a matrix whose
nonzero elements all have the same value in the CSC format.
"""
struct VariableTargetAdjacency{T}
    """
    Column entries: row indices of this column
    """
    colentries::Vector{T}

    """
    Start indices of respective column data in `colentries`
    """
    colstart::Vector{T}
end

"""
    $(SIGNATURES)

    Comparison of two adjacencies
"""
function Base.:(==)(a::VariableTargetAdjacency{Ta}, b::VariableTargetAdjacency{Tb}) where {Ta,Tb}
    Ta==Tb && a.colentries==b.colentries &&  a.colstart==b.colstart
end

"""
$(TYPEDSIGNATURES)

Create an empty VariableTargetAdjacency
"""
VariableTargetAdjacency(t::Type{T}) where T=VariableTargetAdjacency{T}(Vector{T}(undef,0),[one(T)])

"""
$(TYPEDSIGNATURES)

Create an empty VariableTargetAdjacency with default type
"""
VariableTargetAdjacency()=VariableTargetAdjacency(Int64)

"""
$(TYPEDSIGNATURES)

Create a VariableTargetAdjacency from Matrix
"""
VariableTargetAdjacency(m::Matrix{T}) where T=VariableTargetAdjacency{T}(vec(m),collect(1:size(m,1):size(m,1)*size(m,2)+1))

"""
$(TYPEDSIGNATURES)

Show adjacency (in trasposed form; preliminary)
"""
function Base.show(io::IO,adj::VariableTargetAdjacency)
    for isource=1:num_sources(adj)
        for itarget=1:num_targets(adj,isource)
            print(io,adj[itarget,isource], " ")
        end
        println(io, " ");
    end
end

"""
$(TYPEDSIGNATURES)

Access adjacency as if it is a 2D Array
"""
Base.getindex(adj::VariableTargetAdjacency,i,isource)=adj.colentries[adj.colstart[isource]+i-1]
Base.getindex(adj::VariableTargetAdjacency,::Colon,isource)=adj.colentries[adj.colstart[isource]:adj.colstart[isource+1]-1]

"""
$(TYPEDSIGNATURES)

Number of targets for given source
"""
num_targets(adj::VariableTargetAdjacency,isource)=adj.colstart[isource+1]-adj.colstart[isource]

"""
$(TYPEDSIGNATURES)

Number of sources in adjacency
"""
num_sources(adj::VariableTargetAdjacency)=length(adj.colstart)-1

"""
$(TYPEDSIGNATURES)

Number of targeta
"""
num_targets(adj::VariableTargetAdjacency)=maximum(adj.colentries)

"""
$(TYPEDSIGNATURES)

Number of links
"""
num_links(adj::VariableTargetAdjacency)=length(adj.colentries)

"""
$(TYPEDSIGNATURES)

Maximum number of targets per source
"""
max_num_targets_per_source(adj::VariableTargetAdjacency)=maximum(adj.colstart[2:end].-adj.colstart[1:end-1])

Base.size(adj::VariableTargetAdjacency)=(max_num_targets_per_source(adj),num_sources(adj))

function Matrix(adj::VariableTargetAdjacency{T}) where T
    m=zeros(T,size(adj)...)
    for isrc=1:num_sources(adj)
        for itgt=1:num_targets(adj,isrc)
            m[itgt,isrc]=adj[itgt,isrc]
        end
    end
    m
end
"""
$(TYPEDSIGNATURES)

Append a column to adjacency.
"""
function Base.append!(adj::VariableTargetAdjacency,column)
    for i=1:length(column)
        push!(adj.colentries,column[i])
    end
    push!(adj.colstart,length(adj.colentries)+1)
end

"""
$(TYPEDEF)

Use Matrix to store fixed target adjacency
"""
const FixedTargetAdjacency{T} = Matrix{T}

"""
$(TYPEDSIGNATURES)

Number of targets per source if adjacency is a matrix
"""
num_targets(adj::FixedTargetAdjacency,isource)=size(adj)[1]

"""
$(TYPEDSIGNATURES)

Number of sources in adjacency
"""
num_sources(adj::FixedTargetAdjacency)=size(adj)[2]

"""
$(TYPEDSIGNATURES)

Overall number of targets
"""
num_targets(adj::FixedTargetAdjacency)=maximum(vec(adj))

"""
$(TYPEDSIGNATURES)

Number of entries
"""
num_links(adj::FixedTargetAdjacency)=length(adj)

"""
$(TYPEDSIGNATURES)

Maximum number of targets per source
"""
max_num_targets_per_source(adj::FixedTargetAdjacency)=size(adj,1)

"""
$(TYPEDEF)

Adjacency type as union of FixedTargetAdjacency and VariableTargetAdjacency
"""
const Adjacency{T}=Union{FixedTargetAdjacency{T},VariableTargetAdjacency{T}}

"""
$(TYPEDSIGNATURES)

Constructors for Adjacency
"""
Adjacency{T}(a::FixedTargetAdjacency{T}) where T =a
Adjacency{T}(a::VariableTargetAdjacency{T}) where T =a

"""

Transpose adjacency
"""
function atranspose(adj::Adjacency{T}) where T
    # 0th pass: calculate number of rows !!! todo: how to call ?
    t_adj=VariableTargetAdjacency(zeros(T,num_links(adj)),zeros(T,num_targets(adj)+1))

    # 1st pass: calculate new column sizes, store them in t_adj.colstart
    for isource=1:num_sources(adj)
        for itarget=1:num_targets(adj,isource)
            t_adj.colstart[adj[itarget,isource]]+=1
        end
    end

    # 2nd pass: calculate initial assembly
    # addresses in adj.ja by summing up deltas
    # store them a the ends of  the columns:
    # in t_adj.ja[t_adj.colstart[isource+1]-1].
    # it will be overwritten when the last element of
    # the column is assembled.
    # We get automatically increasing  column indices.

    delta=t_adj.colstart[1];
    t_adj.colstart[1]=1;
    for isource=2:num_sources(t_adj)
	save=t_adj.colstart[isource];
	t_adj.colstart[isource]=t_adj.colstart[isource-1]+delta;
	if delta>0
            t_adj.colentries[t_adj.colstart[isource]-1]=t_adj.colstart[isource-1];
        end
	delta=save;
    end
    isource=num_sources(t_adj)+1
    t_adj.colstart[isource]=t_adj.colstart[isource-1]+delta;
    if delta>0
        t_adj.colentries[t_adj.colstart[isource]-1]=t_adj.colstart[isource-1];
    end

    # 3rd pass: assemble new columns
    for isource=1:num_sources(adj)
        for itarget=1:num_targets(adj,isource)
	    asm_idx=t_adj.colstart[adj[itarget,isource]+1]-1;
	    asm_loc=t_adj.colentries[asm_idx];
	    t_adj.colentries[asm_idx]+=1;
	    t_adj.colentries[asm_loc]=isource;
        end
    end
    return t_adj
end

"""
$(TYPEDSIGNATURES)

Try to turn variable target adjacency into fixed target adjacency
"""
function tryfix(a::Adjacency{T}) where T
    ntargets=num_targets(a,1)
    for i=1:num_sources(a)
        if num_targets(a,i)!=ntargets
            return a
        end
    end
    FixedTargetAdjacency(reshape(a.colentries,(ntargets,num_sources(a))))
end

"""
$(TYPEDSIGNATURES)

Turn fixed target adjacency into variable target adjacency
"""
function makevar(a::FixedTargetAdjacency{T}) where T
    ntargets=num_targets(a,1)
    nsources=num_sources(a)
    colstart=T[i*ntargets+1 for i=0:nsources]
    VariableTargetAdjacency(vec(a),colstart)
end

"""
$(TYPEDSIGNATURES)

Create sparse incidence matrix from adjacency
"""
function asparse(a::VariableTargetAdjacency)
    n=num_sources(a)
    m=num_targets(a)
    nzval=ones(Int,length(a.colentries))
    SparseMatrixCSC(m,n,a.colstart,a.colentries,nzval)
end

"""
$(TYPEDSIGNATURES)

Create sparse incidence matrix from adjacency
"""
asparse(a::FixedTargetAdjacency)=asparse(makevar(a))

"""
$(TYPEDSIGNATURES)

Create variable target adjacency from adjacency matrix
"""
VariableTargetAdjacency(m::SparseMatrixCSC{Tv,Ti}) where {Tv<:Integer, Ti<:Integer} = VariableTargetAdjacency{Ti}(m.rowval,m.colptr)


text:
using DeferredAcceptance
using Plots

m = 6

qualities = round.(rand(m), digits=3)
sort!(qualities, rev=true)
capacities = 0.8 * ones(m) / m
cutoffs_container = zeros(m, 1)

"""
Example of a custom demand function. This is actually
equivalent to demands_MNL_onetest(), provided in the DeferredAcceptance module.
"""
function MNL_single_score_demand(cutoffs)
    global cutoffs_container
    cutoffs_container = hcat(cutoffs_container, cutoffs)

    sort_order = sortperm(cutoffs)
    cutoffs[sort_order]

    γ = exp.(qualities)
    demands = zeros(m)

    prob_of_th = diff([cutoffs[sort_order]; 1])

    for c in 1:m, d in c:m     # For each score threshold
        demands[sort_order[c]] += prob_of_th[d] *
                                  γ[sort_order[c]] / sum(γ[sort_order[1:d]])
    end

    return demands
end

nonatomictatonnement(MNL_single_score_demand, capacities, maxit=15)

p = plot(cutoffs_container'[2:end, :],
         label = reshape(["p$i, δ = $(qualities[i])" for i in 1:m], 1, :),
         ls = [:dashdot :dot :solid],
         c = [:crimson :dodgerblue :rebeccapurple :olivedrab :goldenrod :teal],
         title = "Tâtonnement process with $m schools, \none score, equal capacities, MNL choice",
         titlefontsize = 12,
         xlabel = "iteration",
         ylabel = "score cutoff",
         ylim = (0, 1))

# savefig(p, string("examples/MNLuniform/plot.pdf"))
# savefig(p, string("examples/MNLuniform/plot.png"))

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
# -*- coding: utf-8 -*-
# ## Agent Based Modeling
# This script is a self contained agent based model, we define a modeling microframework and then use it to implement an agent based model. This agent based model represents the SIR model of infectious disease with discrete time. This model is not intended as a realistic model of disease, but as an example of how to reverse engineer a model structure from code that implements it. See the agentgraft.jl script for an example of model augmentation on this class of models.
#
# Agent based models are useful for representing computations that cannot be done in closed form and need to be simulated. These simulations arise when analyzing complex systems, such as road networks, social behavior, or very small scale physical systems. The basic principles of ABM are to have agents that can be in states, and transitions between those states, the model then advances the simulation forward in time to get from an initial condition to a final condition. This framework is based on synchronous ABM where there is a global clock and all agents update their state once per clock tick.

module AgentModels

import Base: count

# """    AgentModel

# the root type for agent based models.

# See also: StateModel
# """

abstract type AgentModel end

# """     StateModel

# holds the components of an agent based simulation using finite state machines.

# - states: a collection of distinct states an agent can occupy
# - agents: a collection of states `aᵢ = agents[i] ∈ states` indicating that agent `i` is in state `aᵢ`
# - transitions: the functions `f: states -> states`
# """

mutable struct StateModel <: AgentModel
    states
    agents
    transitions
    loads
end

# +
# Determine number of agents in a given state
function count(sm::StateModel, state)
    return length(filter(a->(a==state), sm.agents))
end

# counts the number of agents
function count(sm::StateModel)
    return length(sm.agents)
end
# -

# stateload computes the fraction of agents in each state
# it is used by tick! to update the statemodel for computing the probability of infection.
function stateload(sm::StateModel, state::Symbol)
    return (count(sm, state)+1)/(count(sm)+1)
end

#     tick!(sm::StateModel)
#
# performs the operations that need to happen once per time step. Things like caching values that need to be computed at the beginning of the timestep like the distribution of the agents' states.

function tick!(sm::StateModel)
    sm.loads = map(s->stateload(sm, s), sm.states)
    return sm.loads
end

#     step!(sm::StateModel, n)
#
# advance the simulation by `n` ticks of time.
# This is an in-place operation that modifies the current state of the simulation.
#

function step!(sm::StateModel, n)
    for s in 1:n
      tick!(sm)
      for (i, a) in enumerate(sm.agents)
          sm.agents[i] = sm.transitions[a](sm, i,a)
      end
        describe(sm)
    end
    return sm
end

# """    describe(sm::StateModel)

# summarize the state of the simulation for presentation or analysis.
# """

function describe(sm::StateModel)
    counts = zeros(Int, size(sm.states))
    d = Dict{eltype(sm.states), Int}()
    for (i,s) in enumerate(sm.states)
        d[s] = i
    end
    for a in sm.agents
        i = d[a]
        counts[i] += 1
    end
    return collect(map(x->Pair(x...), zip(sm.states, counts)))
end

# ## Run the model
#
# This script defines a basic agent based model of disease spread called SIRS. Each agent is in one of 3 states
#
# 1. $S$ Susceptible
# 2. $I$ Infected
# 3. $R$ Recovered
#
# <img src="https://docs.google.com/drawings/d/e/2PACX-1vSeA7mAQ-795lLVxCWXzbkFQaFOHMpwtB121psFV_2cSUyXPyKMtvDjssia82JvQRXS08p6FAMr1hj1/pub?w=1031&h=309">
#
# The agents go from `S->I`, `I-R`, and `R->S` based on random numbers. The probability of S-> is dependent on the fraction of agents in state :I. The probability of recovering is a constant ρ, and the disease confers some temporary immunity with probability μ.
#
# This script has an entrypoint to call it so that you can include this file and run as many simulations as you want. The intended use case is to repeatedly call `main` and accumulate the return values into an array for later analysis.

function main(nsteps)
    n = 20
    a = fill(:S, n)
    ρ = 0.5 + randn(Float64)/4 # chance of recovery
    μ = 0.5 # chance of immunity
    T = Dict(
        :S=>(x...)->rand(Float64) < stateload(x[1], :I) ? :I : :S,
        :I=>(x...)->rand(Float64) < ρ ? :I : :R,
        :R=>(x...)->rand(Float64) < μ ? :R : :S,
    )

    sam = StateModel([:S, :I, :R], a, T, zeros(Float64,3))
    newsam = step!(deepcopy(sam), nsteps)
    @show newsam.agents
    counts = describe(newsam)
    return newsam, counts
end

# +
# # An example of how to run this thing.
# main(10)
# -

end


text:
using ConstrainedDynamics
using Rotations

box = Box(rand(4)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
Body(box)
@test true
cylinder = Cylinder(rand(3)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
Body(cylinder)
@test true
sphere = Sphere(rand(2)...; color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
Body(sphere)
@test true
mesh = Mesh("test.obj", rand(), rand(3,3); color = RGBA(rand(3)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
Body(mesh)
@test true