text:
import Pkg # stdlib

function version_codename_string(x::TomlFile)::String
    toml_file_parsed::Dict{String, Any} = parse_toml_file(x)
    version_codename_string::String = toml_file_parsed["version_codename"]
    return version_codename_string
end

function version_codename_string()::String
    predictmd_toml_file::TomlFile = TomlFile(
        package_directory("Project.toml")
        )
    resultversion_codename_string::String = version_codename_string(
        predictmd_toml_file
        )
    return resultversion_codename_string
end

# function version_codename_string(m::Method)::String
#     m_package_directory::String = package_directory(m)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_codename_string::String = version_codename_string(
#         m_toml_file
#         )
#     return resultversion_codename_string
# end

# function version_codename_string(f::Function)::String
#     m_package_directory::String = package_directory(f)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_codename_string::String = version_codename_string(
#         m_toml_file
#         )
#     return resultversion_codename_string
# end

# function version_codename_string(f::Function, types::Tuple)::String
#     m_package_directory::String = package_directory(f, types)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_codename_string::String = version_codename_string(
#         m_toml_file
#         )
#     return resultversion_codename_string
# end

# function version_codename_string(m::Module)::String
#     m_package_directory::String = package_directory(m)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_codename_string::String = version_codename_string(
#         m_toml_file
#         )
#     return resultversion_codename_string
# end

"""
    version_codename()::String

Return the version code name of PredictMD.
"""
function version_codename()::String
    resultversion_codename_string::String = version_codename_string()
    return resultversion_codename_string
end

"""
    version_codename(m::Method)::String

If method `m`
is part of a Julia package, returns the version code name of that package.

If method `m`
is not part of a Julia package, throws an error.
"""
# function version_codename(m::Method)::String
#     resultversion_codename_string::String = version_codename_string(m)
#     return resultversion_codename_string
# end

"""
    version_codename(f::Function)::String

If function `f`
is part of a Julia package, returns the version code name of
that package.

If function `f`
is not part of a Julia package, throws an error.
"""
# function version_codename(f::Function)::String
#     resultversion_codename_string::String = version_codename_string(f)
#     return resultversion_codename_string
# end

"""
    version_codename(f::Function, types::Tuple)::String

If function `f` with type signature `types`
is part of a Julia package, returns the version code name of
that package.

If function `f` with type signature `types`
is not part of a Julia package, throws an error.
"""
# function version_codename(f::Function, types::Tuple)::String
#     resultversion_codename_string::String = version_codename_string(
#         f,
#         types,
#         )
#     return resultversion_codename_string
# end

"""
    version_codename(m::Module)::String

If module `m` is part of a Julia package, returns the version code name of
that package.

If module `m` is not part of a Julia package, throws an error.
"""
# function version_codename(m::Module)::String
#     resultversion_codename_string::String = version_codename_string(m)
#     return resultversion_codename_string
# end


text:
#for every different type of linkfunction, internally julia will find the right apply_inverse_link function

# This super type of all response distribution types
abstract type InverseLinkFunction end

##LINK FUNCTIONS##

"""inverse cauchit link."""

struct CauchitLink <: InverseLinkFunction
end

function cauchit_inverse_link(x)
  return atan(x) / pi + one(x) / 2
end

"""inverse cloglog link."""

struct CloglogLink <: InverseLinkFunction
end

function cloglog_inverse_link(x)
  return one(x) - exp(-exp(x))
end

"""inverse identity link."""

struct IdentityLink <: InverseLinkFunction
end

function identity_inverse_link(x)
  return x
end

"""inverse inverse link."""

struct InverseLink <: InverseLinkFunction
end

function inverse_inverse_link(x)
  return one(x) / x
end

"""inverse logit link."""

struct LogitLink <: InverseLinkFunction
end

function logit_inverse_link(x)
  return one(x) / (one(x) + exp(-x))
end

"""inverse log link."""

struct LogLink <: InverseLinkFunction
end

function log_inverse_link(x)
  return exp(x)
end

"""inverse probit link."""

struct ProbitLink <: InverseLinkFunction
end

function probit_inverse_link(x)
  return (one(x) + erf(x / sqrt(2 * one(x)))) / 2
end

"""inverse sqrt link."""

struct SqrtLink <: InverseLinkFunction
end

function sqrt_inverse_link(x)
  return x * x
end

##APPLY INVERSE LINK FUNCTIONS

apply_inverse_link(μ, link::LogLink) = log_inverse_link.(μ)

apply_inverse_link(μ, link::IdentityLink) = identity_inverse_link.(μ)

apply_inverse_link(μ, link::SqrtLink) = sqrt_inverse_link.(μ)

apply_inverse_link(μ, link::ProbitLink) = probit_inverse_link.(μ)

apply_inverse_link(μ, link::LogitLink) = logit_inverse_link.(μ)

apply_inverse_link(μ, link::InverseLink) = inverse_inverse_link.(μ)

apply_inverse_link(μ, link::CauchitLink) = cauchit_inverse_link.(μ)

apply_inverse_link(μ, link::CloglogLink) = cloglog_inverse_link.(μ)

ext: jl

max_stars_repo_name: sarah-ji/Tutorials

max_issues_repo_name: sarah-ji/Tutorials

max_forks_repo_name: sarah-ji/Tutorials


text:
using DrWatson
@quickactivate "GCN_HM_GRN-Integration"
using Turing, PyCall

pushfirst!(PyVector(pyimport("sys")."path"), "/gpfs_home/spate116/singhlab/GCN_Integration/notebooks/VI/");
train_model = pyimport("gcn").train_model
run_model = pyimport("gcn").run_model
extract_subgraph = pyimport("gcn").extract_subgraph
pyimport("sys")."stdout" = PyTextIO(stdout)
pyimport("sys")."stderr" = PyTextIO(stderr)

using Zygote
Zygote.@adjoint function pycall(f, x...; kw...)
    x = map(py, x)
    y = pycall(f, x...; kw...)
    y.detach().numpy(), function (ȳ)
        y.backward(gradient = py(ȳ))
        (nothing, map(x->x.grad.numpy(), x)...)
    end
end

node = 10
model, x, y, edge_index = train_model();
edge_index, node_idx = extract_subgraph(node, 2, edge_index)
edge_mask = repeat([1], get(edge_index."shape", 1));

using Flux
Turing.setadbackend(:zygote)
@model graph_mask(y) = begin
    N = get(edge_index."shape", 1)
    edge_mask = Vector{Int8}(undef, N)
    edge_mask ~ filldist(Bernoulli(0.1), N)
    pred = run_model(edge_mask, edge_index, model, node_idx)
    y = y/sum(y)
    y ~ MvNormal(pred/sum(pred), repeat([0.0001], length(y)))
end

@time s = sample(graph_mask(run_model(repeat([1], get(edge_index."shape", 1)), edge_index, model, node_idx)),
        PG(5), 1000)

using DataFrames
@show DataFrame(describe(s)[1])

ei = edge_index."detach"()."cpu"()."numpy"()

base = run_model(repeat([1], get(edge_index."shape", 1)), edge_index, model, node_idx)
@show base/sum(base)

ŷ = run_model([x > 0.1 ? 1 : 0 for x ∈ DataFrame(describe(s)[1])[!, 2]], edge_index, model, node_idx)
@show ŷ/sum(ŷ)

using GraphRecipes
using Plots
using LightGraphs
theme(:default)
gr(fmt=:svg)

uni = unique(ei);
conv = Dict(uni[i] => i for i ∈ 1:length(uni));
conv_rev = Dict(i => uni[i] for i ∈ 1:length(uni));
vals = DataFrame(describe(s)[1])[!, 2];
weigh = Dict((ei[1, i], ei[2, i]) => vals[i] for i ∈ 1:length(vals))

ei_conv = map(x -> conv[x], ei)
inp = ei_conv
classes = y."detach"()."cpu"()."numpy"()[unique(map(x -> conv_rev[x], inp))]

g = DiGraph(Edge.(zip(inp[1, :], inp[2, :])))

plt = plot(graphplot(g, names=uni, arrow=true, nodecolor=classes, edgewidth=(s, d, w)->
    2*weigh[(conv_rev[s],conv_rev[d])], fontsize=10), size=(650, 650))
display(plt)
savefig("/gpfs_home/spate116/singhlab/GCN_Integration/notebooks/VI/out.png")


text:
module TestRecode
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType

if VERSION < v"0.7.0-"
    using CategoricalArrays: replace!
end

const ≅ = isequal

## Test recode!, used by recode

# Test both recoding into x itself and into an uninitialized vector
# (since for CategoricalVectors possible bugs can happen when working in-place)

@testset "Recoding from $(typeof(x)) to $(typeof(y))" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y === z
    @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [6, 7, 8, 100, 0, -1]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(y)) with duplicate recoded values" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, 1=>100, 2:4=>100, [5; 9:10]=>-1)
    @test y === z
    @test y == [100, 100, 100, 100, -1, 6, 7, 8, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [6, 7, 8, 100, -1]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(y)) with unused level" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1, 100=>1)
    @test y === z
    @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [6, 7, 8, 100, 0, -1, 1]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(y)) with duplicate default" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, 100, 1=>100, 2:4=>100, [5; 9:10]=>-1)
    @test y === z
    @test y == [100, 100, 100, 100, -1, 100, 100, 100, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [100, -1]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(y)) with default" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, -10, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y === z
    @test y == [100, 0, 0, 0, -1, -10, -10, -10, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [100, 0, -1, -10]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(y)) with first value being Float64" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, 1.0=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [6, 7, 8, 100, 0, -1]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(y)) with overlapping pairs" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1, 1:10=>0)
    @test y === z
    @test y == [100, 0, 0, 0, -1, 0, 0, 0, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [100, 0, -1]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(y)) with changes to levels order" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Int}(undef, size(x)),
          CategoricalArray{Int}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y === z
    @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(y, CategoricalArray)
        @test levels(y) == [6, 7, 8, 100, 0, -1]
        @test !isordered(y)
    end
end

@testset "Recoding from $(typeof(x)) to categorical array with missing values" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"]))

    # check that error is thrown
    y = Vector{String}(undef, 4)
    @test_throws MissingException recode!(y, x, "a", "c"=>"b")

    y = CategoricalVector{String}(undef, 4)
    @test_throws MissingException recode!(y, x, "a", "c"=>"b")
end

@testset "Recoding array with missings and default from $(typeof(x)) to $(typeof(y))" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"])),
    y in (similar(x), Array{Union{String, Missing}}(undef, size(x)),
          CategoricalArray{Union{String, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, "a", "c"=>"b")
    @test y === z
    @test y ≅ ["a", missing, "b", "a"]
    if isa(y, CategoricalArray)
        @test levels(y) == ["b", "a"]
        @test !isordered(y)
    end
end

@testset "Recoding array with missings and no default from $(typeof(x)) to $(typeof(y))" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"])),
    y in (similar(x), Array{Union{String, Missing}}(undef, size(x)),
          CategoricalArray{Union{String, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, "c"=>"b")
    @test y === z
    @test y ≅ ["a", missing, "b", "d"]
    if isa(y, CategoricalArray)
        @test levels(y) == ["a", "d", "b"]
        @test !isordered(y)
    end
end

@testset "Collection in LHS recoding array with missings and no default from $(typeof(x)) to $(typeof(y))" for
    x in (["1", missing, "3", "4", "5"], CategoricalArray(["1", missing, "3", "4", "5"])),
    y in (similar(x), Array{Union{String, Missing}}(undef, size(x)),
          CategoricalArray{Union{String, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, ["3","4"]=>"2")
    @test y === z
    @test y ≅ ["1", missing, "2", "2", "5"]
    if isa(y, CategoricalArray)
        @test levels(y) == ["1", "5", "2"]
        @test !isordered(y)
    end
end

@testset "Recoding array with missings, default and with missing as a key pair from $(typeof(x)) to $(typeof(y))" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"])),
    y in (similar(x), Array{Union{String, Missing}}(undef, size(x)),
          CategoricalArray{Union{String, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, "a", "c"=>"b", missing=>"d")
    @test y === z
    @test y == ["a", "d", "b", "a"]
    if isa(y, CategoricalArray)
        @test levels(y) == ["b", "d", "a"]
        @test !isordered(y)
    end
end

@testset "Collection with missing in LHS recoding array with missings, default from $(typeof(x)) to $(typeof(y))" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"])),
    y in (similar(x), Array{Union{String, Missing}}(undef, size(x)),
          CategoricalArray{Union{String, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, "a", [missing, "c"]=>"b")
    @test y === z
    @test y == ["a", "b", "b", "a"]
    if isa(y, CategoricalArray)
        @test levels(y) == ["b", "a"]
        @test !isordered(y)
    end
end

@testset "Recoding array with missings, no default and with missing as a key pair from $(typeof(x)) to $(typeof(y))" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"])),
    y in (similar(x), Array{Union{String, Missing}}(undef, size(x)),
          CategoricalArray{Union{String, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, "c"=>"b", missing=>"d")
    @test y === z
    @test y == ["a", "d", "b", "d"]
    if isa(y, CategoricalArray)
        @test levels(y) == ["a", "b", "d"]
        @test !isordered(y)
    end
end

@testset "Collection with missing in LHS recoding array with missings, no default from $(typeof(x)) to $(typeof(y))" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"])),
    y in (similar(x), Array{Union{String, Missing}}(undef, size(x)),
          CategoricalArray{Union{String, Missing}}(undef, size(x)), x)

    z = @inferred recode!(y, x, ["c", missing]=>"b")
    @test y === z
    @test y == ["a", "b", "b", "d"]
    if isa(y, CategoricalArray)
        @test levels(y) == ["a", "d", "b"]
        @test !isordered(y)
    end
end

@testset "Recoding into an array of incompatible size from $(typeof(x)) to $(typeof(y))" for
    x in (["a", missing, "c", "d"], CategoricalArray(["a", missing, "c", "d"])),
    y in (similar(x, 0), Array{Union{String, Missing}}(undef, 0),
          CategoricalArray{Union{String, Missing}}(undef, 0))

    @test_throws DimensionMismatch recode!(y, x, "c"=>"b", missing=>"d")
end

@testset "Recoding into an array with incompatible eltype from $(typeof(x)) to $(typeof(y))" for
    x in ([1:10;], CategoricalArray(1:10)),
    y in (similar(x, String), Array{String}(undef, size(x)),
          CategoricalArray{String}(undef, size(x)))

    @test_throws ArgumentError recode!(y, x, 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
end

@testset "Recoding into an array with incompatible eltype from $(typeof(x)) to $(typeof(y))" for
    x in ((Union{Int, Missing})[1:10;], CategoricalArray{Union{Int, Missing}}(1:10)),
    y in (similar(x), Array{Union{Int, Missing}}(undef, size(x)),
          CategoricalArray{Union{Int, Missing}}(undef, size(x)))

    @test_throws MethodError recode!(y, x, 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
end

## Test in-place recode!()

@testset "Recoding from $(typeof(x)) to $(typeof(x)) without default" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10))

    z = @inferred recode!(x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test x === z
    @test x == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(x, CategoricalArray)
        @test levels(x) == [6, 7, 8, 100, 0, -1]
        @test !isordered(x)
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(x)) without default" for
    x in ([1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10))

    z = @inferred recode!(x, 1, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test x === z
    @test x == [100, 0, 0, 0, -1, 1, 1, 1, -1, -1]
    if isa(x, CategoricalArray)
        @test levels(x) == [100, 0, -1, 1]
        @test !isordered(x)
    end
end

@testset "recode() promotion for $(typeof(x))" for
    x in (1:10, [1:10;], CategoricalArray(1:10), CategoricalArray{Union{Int, Missing}}(1:10))
    T = eltype(x) >: Missing ? Missing : Union{}

    # Recoding from Int to Float64 due to a second value being Float64
    y = @inferred recode(x, 1=>100.0, 2:4=>0, [5; 9:10]=>-1)
    @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Union{Float64, T}, DefaultRefType})
        @test levels(y) == [6, 7, 8, 100, 0, -1]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Union{Float64, T}}
    end

    # Recoding from Int to Float64, with Float64 default and all other values Int
    y = @inferred recode(x, -10.0, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y == [100, 0, 0, 0, -1, -10, -10, -10, -1, -1]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Union{Float64, T}, DefaultRefType})
        @test levels(y) == [100, 0, -1, -10]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Union{Float64, T}}
    end

    # Recoding from Int to Any
    y = @inferred recode(x, 1=>"a", 2:4=>0, [5; 9:10]=>-1)
    @test y == ["a", 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Any, DefaultRefType})
        @test levels(y) == [6, 7, 8, "a", 0, -1]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Any}
    end

    # Recoding from Int to String, with String default
    y = @inferred recode(x, "d", 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
    @test y == ["a", "b", "b", "b", "c", "d", "d", "d", "c", "c"]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Union{String, T}, DefaultRefType})
        @test levels(y) == ["a", "b", "c", "d"]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Union{String, T}}
    end

    # Recoding from Int to String, with all original levels recoded
    y = @inferred recode(x, 1:4=>"a", [5; 9:10]=>"b", 6:8=>"c")
    @test y == ["a", "a", "a", "a", "b", "c", "c", "c", "b", "b"]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Union{String, T}, DefaultRefType})
        @test levels(y) == ["a", "b", "c"]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Union{String, T}}
    end

    # Recoding from Int to Int/String (i.e. Any), with default String and other values Int
    y = @inferred recode(x, "x", 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y == [100, 0, 0, 0, -1, "x", "x", "x", -1, -1]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Any, DefaultRefType})
        @test levels(y) == [100, 0, -1, "x"]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Any}
    end

    # Recoding from Int to Int/String, without any Int value in pairs
    # and keeping some original Int levels
    # This must fail since we cannot take into account original eltype (whether original
    # levels are kept is only known at run time)
    res = @test_throws ArgumentError recode(x, 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
    @test sprint(showerror, res.value) ==
        "ArgumentError: cannot `convert` value 6 (of type $Int) to type of recoded levels ($(Union{String, T})). " *
        "This will happen with recode() when not all original levels are recoded " *
        "(i.e. some are preserved) and their type is incompatible with that of recoded levels."

    # Recoding from Int to Union{Int, Missing} with missing default
    y = @inferred recode(x, missing, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y ≅ [100, 0, 0, 0, -1, missing, missing, missing, -1, -1]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Union{Int, Missing}, DefaultRefType})
        @test levels(y) == [100, 0, -1]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Union{Int, Missing}}
    end

    # Recoding from Int to Union{Int, Missing} with missing RHS
    y = @inferred recode(x, 1=>missing, 2:4=>0, [5; 9:10]=>-1)
    @test y ≅ [missing, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Union{Int, Missing}, DefaultRefType})
        @test levels(y) == [6, 7, 8, 0, -1]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Union{Int, Missing}}
    end

    # Recoding from Int to Union{Int, Missing} with single missing RHS
    y = @inferred recode(x, 1=>missing)
    @test y ≅ [missing, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Union{Int, Missing}, DefaultRefType})
        @test levels(y) == collect(2:10)
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Union{Int, Missing}}
    end
end

@testset "Recoding from $(typeof(x)) to $(typeof(x))" for
    x in (["a", "c", "b", "a"], CategoricalArray(["a", "c", "b", "a"]))

    y = @inferred recode(x, "c"=>"x", "b"=>"y", "a"=>"z")
    @test y == ["z", "x", "y", "z"]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{String, DefaultRefType})
        @test levels(y) == ["x", "y", "z"]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{String}
    end
end

@testset "Recoding a matrix $(typeof(x))" for
    x in (['a' 'c'; 'b' 'a'], CategoricalArray(['a' 'c'; 'b' 'a']))

    y = @inferred recode(x, 'c'=>'x', 'b'=>'y', 'a'=>'z')
    @test y == ['z' 'x'; 'y' 'z']
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalMatrix{Char, DefaultRefType})
        @test levels(y) == ['x', 'y', 'z']
        @test !isordered(y)
    else
        @test typeof(y) === Matrix{Char}
    end
end

@testset "Recoding from $(typeof(x)) to Int/String (i.e. Any), with index and levels in different orders" for
    x in (10:-1:1, CategoricalArray(10:-1:1))

    y = @inferred recode(x, 0, 1=>"a", 2:4=>"c", [5; 9:10]=>"b")
    @test y == ["b", "b", 0, 0, 0, "b", "c", "c", "c", "a"]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{Any, DefaultRefType})
        @test levels(y) == ["a", "c", "b", 0]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{Any}
    end

    # Recoding from Int to String via default, with index and levels in different orders
    y = @inferred recode(x, "x", 1=>"a", 2:4=>"c", [5; 9:10]=>"b")
    @test y == ["b", "b", "x", "x", "x", "b", "c", "c", "c", "a"]
    if isa(x, CategoricalArray)
        @test isa(y, CategoricalVector{String, DefaultRefType})
        @test levels(y) == ["a", "c", "b", "x"]
        @test !isordered(y)
    else
        @test typeof(y) === Vector{String}
    end
end

@testset "Recoding CategoricalArray with custom reftype" begin
    x = CategoricalVector{Int, UInt8}(1:10)
    y = @inferred recode(x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    @test isa(y, CategoricalVector{Int, UInt8})
    @test levels(y) == [6, 7, 8, 100, 0, -1]
    @test !isordered(y)
end

@testset "Recoding ordered CategoricalArray and merging two categories" begin
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "c"=>"a")
    @test y == ["a", "a", "b", "a"]
    @test isa(y, CategoricalVector{String, DefaultRefType})
    @test levels(y) == ["a", "b"]
    @test isordered(y)
end

@testset "Recoding ordered CategoricalArray and merging one category into another (contd.)" begin
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "b"=>"c")
    @test y == ["a", "c", "c", "a"]
    @test isa(y, CategoricalVector{String, DefaultRefType})
    @test levels(y) == ["a", "c"]
    @test isordered(y)
end

@testset "Recoding ordered CategoricalArray with new level which cannot be ordered" begin
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "b"=>"d")
    @test y == ["a", "c", "d", "a"]
    @test isa(y, CategoricalVector{String, DefaultRefType})
    @test levels(y) == ["a", "c", "d"]
    @test !isordered(y)
end

@testset "Recoding ordered CategoricalArray with conflicting orders" begin
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "b"=>"b", "a"=>"a")
    @test y == ["a", "c", "b", "a"]
    @test isa(y, CategoricalVector{String, DefaultRefType})
    @test levels(y) == ["c", "b", "a"]
    @test !isordered(y)
end

@testset "Recoding ordered CategoricalArray with default already in levels" begin
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "a", "c"=>"b")
    @test y == ["a", "b", "a", "a"]
    @test isa(y, CategoricalVector{String, DefaultRefType})
    @test levels(y) == ["b", "a"]
    @test !isordered(y)
end

@testset "Recoding ordered CategoricalArray with default not in levels" begin
    x = CategoricalArray(["d", "c", "b", "d"])
    ordered!(x, true)
    y = @inferred recode(x, "a", "c"=>"b")
    @test y == ["a", "b", "a", "a"]
    @test isa(y, CategoricalVector{String, DefaultRefType})
    @test levels(y) == ["b", "a"]
    @test !isordered(y)
end

@testset "Recoding CategoricalArray with missings and no default" begin
    x = CategoricalArray{Union{String, Missing}}(["a", "b", "c", "d"])
    x[2] = missing
    y = @inferred recode(x, "c"=>"b")
    @test y ≅ ["a", missing, "b", "d"]
    @test isa(y, CategoricalVector{Union{String, Missing}, DefaultRefType})
    @test levels(y) == ["a", "b", "d"]
    @test !isordered(y)
end

@testset "Recoding CategoricalArray with missings and non-missing default" begin
    x = CategoricalArray{Union{String, Missing}}(["a", "b", "c", "d"])
    x[2] = missing
    y = @inferred recode(x, "a", "c"=>"b")
    @test y ≅ ["a", missing, "b", "a"]
    @test isa(y, CategoricalVector{Union{String, Missing}, DefaultRefType})
    @test levels(y) == ["b", "a"]
    @test !isordered(y)
end

@testset "replace with CategoricalArray" begin
    function testf(replacef, T, x::CategoricalArray{S,R}, pairs::Pair...) where {S, R}
        ca = replacef(x, pairs...)
        @test ca isa CategoricalArray{T, R}

        if VERSION >= v"0.7.0-"
            a = replacef(Array(x), pairs...)
            @test ca ≅ a
        end

        ca
    end

    @testset "strings with missings" begin
        x = categorical(["a", "b", missing, "a"])

        testf(replace, String, x, missing => "")
        testf(replace, Union{String, Missing}, x, "b" => "c")
        testf(replace, Any, x, "a" => 1, "b" => 2)
        testf(replace, Any, x, "a" => 1, "b" => 2, missing => 3)
        y = testf(replace!, Union{String, Missing}, x, "b" => "c")
        @test y === x

        y = testf(replace!, Union{String, Missing}, x, missing => "")
        @test y === x
    end

    @testset "strings without missings" begin
        x = categorical(["a", "b", "", "a"])
        testf(replace, Union{String, Missing}, x, "" => missing)
        @test_throws MethodError replace!(x, "" => missing)

        x = categorical(Union{String, Missing}["a", "b", "", "a"])
        testf(replace!, Union{String, Missing}, x, "" => missing)
    end

    @testset "int to float" begin
        x = categorical([0, 2, 0])

        @testset "replace" begin
            testf(replace, Float64, x, 2 => 1.5, 0 => 0.5)
            testf(replace, Float64, x, 2 => 1.0)
        end

        @testset "replace!" begin
            @test_throws InexactError replace!(x, 2 => 1.5)

            y = testf(replace!, Int, x, 2 => 1.0)
            @test y === x
        end
    end

    @testset "float to int" begin
        x = categorical([0.5, 2.0, 0.5])

        @testset "replace" begin
            testf(replace, Float64, x, 2 => 1)
            testf(replace, Float64, x, 0.5 => 3)
        end

        @testset "replace!" begin
            y = testf(replace!, Float64, x, 2 => 1)
            @test x === y

            y = testf(replace!, Float64, x, 0.5 => 3)
            @test x === x
        end
    end
end

end


text:
@testset "Grid" begin

    nx, ny = 8, 8
    grid = Grid( 1e6, nx, 1e6, ny )

    tx = SurfaceQuasiGeostrophic.fct_unity_approx5(nx)
    ty = SurfaceQuasiGeostrophic.fct_unity_approx5(ny)

    @show tx
    @show ty

    @test true

end


text:
# Dependencies
push!(LOAD_PATH,"./")
using CpelNano

# IO
dataDir = dirname(dirname(pathof(CpelNano)))
nano = "$(dataDir)/examples/full_example/nanopolish/full_example_noise2.0_methylation_calls.sorted.tsv"
fasta = "$(dataDir)/examples/full_example/reference/hg38_chr17_43023997_43145780.fa"
outdir = "$(dataDir)/examples/full_example/cpelnano/"

# Output files
mml_path = "$(outdir)/full_example_mml.cpelnano"
nme_path = "$(outdir)/full_example_nme.cpelnano"
theta_path = "$(outdir)/full_example_theta.cpelnano"

# Remove output in outdir
isfile(mml_path) && rm(mml_path)
isfile(nme_path) && rm(nme_path)
isfile(theta_path) && rm(theta_path)

# Configuration
min_cov = 2.5; max_size_subreg = 500; size_est_reg = 3000; max_em_init = 1; max_em_iters = 20;
config = CpelNano.CpelNanoConfig(min_cov, max_size_subreg, size_est_reg, max_em_init, max_em_iters);
config.out_dir = outdir; config.out_prefix = "full_example"; config.verbose = false;

# Analyze each region in chr
CpelNano.analyze_nano(nano,fasta,config)

# Print files
println("*******************************MML*******************************")
println(String(read(mml_path)))
println("*****************************************************************")
readline()
println("*******************************NME*******************************")
println(String(read(nme_path)))
println("*****************************************************************")
readline()
println("*******************************THETA******************************")
println(String(read(theta_path)))
println("*****************************************************************")
readline()


text:
using JSON, DelimitedFiles
import XLSX

# Read a file with fixed width columns and a header row.
function readfixedwidth(x::String; opts...)
	data, headers = readdlm(x; header=true, opts...)
	return Dict{Symbol,Vector{Float64}}(Symbol(headers[i]) => data[:,i] for i=1:length(headers))
end

function readssp(unit)
	path = joinpath(@__DIR__, "..")
	file = "$path/SSP/SSP database.xlsx"
	models = ["AIM/CGE", "GCAM4", "IMAGE", "MESSAGE-GLOBIOM", "REMIND-MAGPIE", "WITCH-GLOBIOM"]
	gases = Dict(
		"TotalCO2" => "Emissions - CO2 total",
		"OtherCO2" => "Emissions - CO2 land use",
		"CH4" => "Emissions - CH4",
		"N2O" => "Emissions - N2O",
		"Population" => "Population - total",
		"GDP" => "GDP (PPP)"
	)
	regions = Dict(
		"R5.2ASIA" => "Asia",
		"R5.2LAM" => "LAM",
		"R5.2MAF" => "MAF",
		"R5.2OECD" => "OECD",
		"R5.2REF" => "REF",
		"World" => "Global"
	)
	scenarios = ["SSP1-26", "SSP1-34", "SSP1-45", "SSP1-Baseline",
				"SSP2-26", "SSP2-34", "SSP2-45", "SSP2-60", "SSP2-Baseline",
				"SSP3-34", "SSP3-45", "SSP3-60", "SSP3-Baseline",
				"SSP4-26", "SSP4-34", "SSP4-45", "SSP4-60", "SSP4-Baseline",
				"SSP5-26", "SSP5-34", "SSP5-45", "SSP5-60", "SSP5-Baseline"
	]

	dict(n) = n == 0 ? Vector{Float64} : Dict{String, dict(n-1)}
	ssp = dict(4)()
	for (gas, sheet) in gases
		nrows = if gas == "Population"
					781
				elseif gas == "GDP"
					721
				else
					631
				end
		arr = XLSX.readdata(file, sheet, "A1:P$nrows")
		unitfactor = if gas == "TotalCO2" || gas == "OtherCO2"
						(unit == :GtC ? 12/44 : 1.0) / 1000
					elseif gas == "N2O"
						28/44/1000
					elseif gas == "Population"
						0.001
					else
						1.0
					end
		for r=2:nrows
			model, scen, reg = arr[r,1], arr[r,2], regions[arr[r,3]]
			#ssp[gas][model][scen][reg] = float.(arr[r, 6:16]) * unitfactor
			if model in models
				get!(dict(1), get!(dict(2), get!(dict(3), ssp, gas), model), scen)[reg] = signif.(float.(arr[r, 6:16]) * unitfactor, 6)
			end
		end
		for model in models, scen in scenarios
			!haskey(ssp[gas][model], scen) && continue
			sspdict = ssp[gas][model][scen]
			sspdict["Global"] = signif.(sspdict["OECD"] + sspdict["Asia"] + sspdict["LAM"] + sspdict["MAF"] + sspdict["REF"], 6)
		end
	end
	open("SSPscenarios.js", "w") do f
		write(f, JSON.json(ssp, 4))
	end
	return ssp		# the return vector (in the dict) is for [2005, 2010, 2020, ..., 2100]
end

#=
	then find/replace in sublime:
		\s+(\-?\d+\.\d+,*)\n	=> 	$1
		\s+]					=>	]
	finally add var assignment to beginning
=#

# Usage: ClimateCalculator_new.readrcp("RCP3PD", ["FossilCO2","OtherCO2"])
function readrcp(rcpscen,gases,unit)
	path = joinpath(dirname(@__FILE__), "..")
	file = "$path/RCP/$rcpscen.SCEN"
	lines = parse(Int, readline(file))

	# note: readdlm() skips blank lines
	data = readdlm(file, skipstart=6)
	regions = data[1:lines+3:end-1,1]
	headers = data[2,:]
	firstyear = data[4,1]
	endyear = data[4+lines-1,1]
	#rcp = Dict{String, Dict{String, Vector{Float64}}}()
	rcp = Dict(gas => Dict{String, Vector{Float64}}() for gas in [gases; "TotalCO2"])
	for gas in gases, (i,reg) in enumerate(regions)
		vect = Array{Float64}(endyear-firstyear+1)
		colindex = findfirst(headers, gas)
		colindex == 0 && error("Column $gas not found in $headers.")
		lastrow = 4 + (i-1)*(lines+3)
		for row = (5:4+lines-1) + (i-1)*(lines+3)
			years = data[lastrow,1]:data[row,1]
			#println(i," ",reg," ",lastrow," ",row," ",years)
			vect[years-firstyear+1] = linspace(data[lastrow,colindex],data[row,colindex],length(years))
			lastrow = row
		end
		rcp[gas][reg] = unit == :GtCO2 ? vect*44/12 : vect
	end
	for reg in regions
		rcp["TotalCO2"][reg] = rcp["FossilCO2"][reg] + rcp["OtherCO2"][reg]
	end
	rcp
end

function readallrcp()
	rcplist = ["RCP3PD", "RCP45", "RCP6", "RCP85"]
	gases = ["FossilCO2", "OtherCO2", "CH4", "N2O"]
	allrcp = Dict(r => readrcp(r,gases,:GtCO2) for r in rcplist)
	open("RCPscenarios.js", "w") do f
		write(f, JSON.json(allrcp, 4))
	end
	#=open("testtest.txt", "w") do f
		show(IOContext(f, :limit => false), MIME("text/plain"), allrcp)
	end=#
	nothing
end

function historicdata()
	rcplist = ["RCP3PD", "RCP45", "RCP6", "RCP85"]
	path = joinpath(dirname(@__FILE__), "..")
	conc_RCP = Dict(r => readfixedwidth("$path/RCP/$(r)_MIDYEAR_CONCENTRATIONS.DAT", skipstart=38) for r in rcplist)
	forcing_RCP = Dict(r => readfixedwidth("$path/RCP/$(r)_MIDYEAR_RADFORCING.DAT", skipstart=59) for r in rcplist)
	temperatureAnomaly = readfixedwidth("$path/AnnualTemperatures.dat")

	forcing_conc_RCP = Dict(r => Dict(g => Vector{Float64}(undef, length(YEARS)) for g in GAS) for r in rcplist)
	RadiativeForcing = Dict{Symbol,Float64}(g => 0.0 for g in GAS)

	for r in rcplist
		for i=1:length(YEARS)
			conc = Dict{Symbol,Float64}(g => conc_RCP[r][g][i] for g in GAS3)
			radiativeforcing!(conc, RadiativeForcing)
			for g in GAS
				forcing_conc_RCP[r][g][i] = RadiativeForcing[g]
			end
		end
	end

	aerosolforcing = Dict(r => forcing_RCP[r][:TOTAER_DIR_RF] + forcing_RCP[r][:CLOUD_TOT_RF] +
								forcing_RCP[r][:BCSNOW_RF] for r in rcplist)

	# solar RF:  average over one cycle is 0.103756
	solarRF_1::Dict{String,Vector{Float64}} = Dict(r => [forcing_RCP[r][:SOLAR_RF][YEARS .< 2010]; fill(0.103756, YEARS[end]-2010+1)] for r in rcplist)

	# about 50% of RF from tropospheric O3 is due to CH4, the rest is here
	otherforcing = Dict(r => solarRF_1[r] + forcing_RCP[r][:VOLCANIC_ANNUAL_RF] + forcing_RCP[r][:LANDUSE_RF] +
								forcing_RCP[r][:FGASSUM_RF] + forcing_RCP[r][:MHALOSUM_RF] + forcing_RCP[r][:STRATOZ_RF] +
								0.5*forcing_RCP[r][:TROPOZ_RF] for r in rcplist)

	#histNonCO2forcing_1 = forcing_RCP[:CH4_RF] + forcing_RCP[:N2O_RF] + forcing_RCP[:CH4OXSTRATH2O_RF] + 0.5*forcing_RCP[:TROPOZ_RF]

	histTempGISS = zeros(length(YEARS))
	histTempGISS[iyear(1880):iyear(2010)] = temperatureAnomaly[:GISS] .+ 0.25

	return aerosolforcing, otherforcing, conc_RCP, forcing_conc_RCP, histTempGISS
end

const aerosolforcing, otherforcing, conc_RCP, forcing_conc_RCP, histTempGISS = historicdata()

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
import Phylo.API: _src, _dst, _getlength, _getbranchname

"""
    Branch

    A branch connecting two AbstractNodes of a phylogenetic tree
"""
mutable struct Branch{RT, NL} <: AbstractBranch{RT, NL}
    name::Int
    source::NL
    destination::NL
    length::Float64

    function Branch{RT}(name::Int, source::NL, destination::NL,
                        length::Float64 = NaN) where {RT, NL}
        length >= 0.0 || isnan(length) ||
            error("Branch length must be positive or NaN (no recorded length), not $length")
        new{RT, NL}(name, source, destination, length)
    end
end

_src(::AbstractTree, branch::Branch) = branch.source
_dst(::AbstractTree, branch::Branch) = branch.destination
_getlength(::AbstractTree, branch::Branch) = branch.length
_getbranchname(::AbstractTree, branch::Branch) = branch.name

function checkbranch(id::Int, branch::Branch{<: Rooted}, tree::AbstractTree)
    return id > 0 &&
        hasbranch(tree, id) &&
        src(tree, branch) != dst(tree, branch) &&
        hasnode(tree, src(tree, branch)) &&
        hasnode(tree, dst(tree, branch)) &&
        indegree(tree, getnode(tree, dst(tree, branch))) > 0 &&
        outdegree(tree, src(tree, branch)) > 0 &&
        getbranch(tree, getinbound(tree, dst(tree, branch))) === branch
end


text:
function dict(data::Vector{DataAttrs})

    data_len=length(data[1].data)
    @assert unique(map(x->length(x.data),data))==[data_len] "All Axis should have the same length!"

    ret_data=[]
    axis_len=length(data)

    for i in 1:data_len
        dot=[]
        for j in 1:axis_len
            push!(dot,data[j].data[i])
        end
        push!(ret_data,dot)
    end

    xaxis=Dict()
    yaxis=Dict()
    for j in 1:axis_len

        if j==1
           el_type=eltype(data[j].data)

           if el_type<:Number
                xaxis["type"]="value"
           elseif el_type<:AbstractString
                xaxis["type"]="category"
           elseif el_type<:Date
                xaxis["type"]="time"
           end
        end

        if j==2
           el_type=eltype(data[j].data)
           if el_type<:Number
                yaxis["type"]="value"
           elseif el_type<:AbstractString
                yaxis["type"]="category"
           elseif el_type<:Date
                yaxis["type"]="time"
           end
        end

    end

    return (data=ret_data,xaxis=xaxis,yaxis=yaxis)
end

function dict(series::Series)

    nt=dict(series.data)
    xaxis=nt.xaxis
    yaxis=nt.yaxis

    series_options=Dict("type"=>series.plot_type,"name"=>series.name,"data"=>nt.data)
    merge!(series_options,series.options)
   return (xaxis=xaxis,yaxis=yaxis,series_options=series_options)

end

function dict_any(d::Dict)
    d=convert(Dict{String,Any},d)
    for (k,v) in d
        v isa Dict ? d[k]=dict_any(v) : nothing
    end
    return d
end

function dict(p::Plot)

    options=p.options
    series_d=Dict("series"=>[])
    for s in p.series
        nt=dict(s)
        push!(series_d["series"],nt.series_options)

        if length(nt.xaxis)!=0
            options["xAxis"]=[nt.xaxis]
        end
        if length(nt.yaxis)!=0
            options["yAxis"]=[nt.yaxis]
        end
    end

    ## legend
    if length(p.series)>1
       options["legend"]=Dict("data"=>map(x->x.name,p.series))
    end

    merge!(options,series_d)

    return dict_any(options)
end


text:
import Base: copy
import Base.Broadcast:
    AbstractArrayStyle, Broadcasted, BroadcastStyle, Style, broadcastable, instantiate,
    flatten

# For broadcasting, we think of Lists as being more complex than a 0-dimensional object but
# less complex than a 1-dimensional one (since we lack size). Note that this doesn’t
# actually work for broadcasting with arrays — since Lists don’t implement the array
# interface. This is not an easy problem to resolve.
BroadcastStyle(::Type{<:List}) = Style{List}()
BroadcastStyle(::Style{List}, ::AbstractArrayStyle{0}) = Style{List}()
BroadcastStyle(::Style{List}, s::AbstractArrayStyle{N}) where {N} = s

broadcastable(α::List) = α
instantiate(bc::Broadcasted{Style{List}}) = bc

# This implementation of copy is essentially the same as the implementation of map, except
# that the “only” elements are treated as infinite lists that are nil whenever necessary.
isnil_or_only(::Cons) = false
isnil_or_only(::Any) = true

car_or_only(α::List) = car(α)
car_or_only(x) = x[]

cdr_or_only(α::List) = cdr(α)
cdr_or_only(x) = x

function map_or_only(f, ::Nil, xs...)
    for x in xs
        if !isnil_or_only(xs)
            throw(ArgumentError("not all lists are of equal length"))
        end
    end
    nil
end
map_or_only(f, x::Cons, xs...) =
    Cons(f(car(x), map(car_or_only, xs)...),
         map_or_only(f, cdr(x), map(cdr_or_only, xs)...))

# Pull the scalar out so we may rely on dispatch to determine if the map is empty
map_or_only(f, x, xs...) = map((ys...) -> f(x, ys...), xs...)

function copy(bc::Broadcasted{Style{List}})
    f = flatten(bc)
    map_or_only(f.f, f.args...)
end


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
BoundConstrainedTrustRegionResult

Note that we re-use problem (1) as defined in the documentation of the
solve_bound_constrained_trust_region function.
"""
mutable struct BoundConstrainedTrustRegionResult
  "The approximate solution to (1)."
  solution::Vector{Float64}
  "The value of objective_vector' * (solution - center_point)."
  value::Float64
end

"""
Finds a solution to the problem:

argmin_x objective_vector' * x
s.t. variable_lower_bounds <= x <= variable_upper_bounds                    (1)
     || x - center_point || <= target_radius

where || . || is weighted by norm_weights.

for a positive value of target_radius, by solving the related problem

argmin_t objective_vector' * x
s.t. x := min(max(center_point - t * objective_vector, variable_lower_bounds),
              variable_upper_bounds)
     || x - center_point || <= target_radius

Note that the definition of x is just applying the lower/upper bound constraints
to center_point - t * objective_vector.

This problem is solved by computing the breakpoint at which each component of
x switches from varying with t to being fixed at its bounds.  The radius of the
median breakpoint is evaluated, eliminating half of the components.  The process
is iterated until the argmin is identified.

# Inputs:
- `center_point::Vector{AbstractFloat64}`. See (1). It is assumed that
  center_point satisfies the bounds.
- `objective_vector::Vector{Float64}`. See (1).
- `variable_lower_bounds::Vector{Float64}`. See (1).
- `variable_upper_bounds::Vector{Float64}`. See (1).
- `norm_weights::AbstractVector{Float64}`. The weights used for the norm
  || . ||.
- `target_radius::Float64`. See (1).
- `solve_approximately::Bool`. If true, solve the approximate version of (1)
  which only applies the bound constraints to variables already at their bound.

# Output:
BoundConstrainedTrustRegionResult struct.
"""
function solve_bound_constrained_trust_region(
  center_point::AbstractVector{Float64},
  objective_vector::Vector{Float64},
  variable_lower_bounds::Vector{Float64},
  variable_upper_bounds::Vector{Float64},
  norm_weights::AbstractVector{Float64},
  target_radius::Float64,
  solve_approximately::Bool,
)
  if solve_approximately
    return approximately_solve_bound_constrained_trust_region(
      center_point,
      objective_vector,
      variable_lower_bounds,
      variable_upper_bounds,
      norm_weights,
      target_radius,
    )
  end

  @assert 0.0 <= target_radius < Inf
  if target_radius == 0.0 || norm(objective_vector, 2) == 0.0
    return BoundConstrainedTrustRegionResult(copy(center_point), 0.0)
  end

  direction = zeros(length(center_point))
  threshold = zeros(length(center_point))
  for idx in 1:length(center_point)
    if center_point[idx] >= variable_upper_bounds[idx] &&
       objective_vector[idx] <= 0
      continue
    end
    if center_point[idx] <= variable_lower_bounds[idx] &&
       objective_vector[idx] >= 0
      continue
    end
    direction[idx] = -objective_vector[idx] / norm_weights[idx]
    if direction[idx] > 0
      threshold[idx] =
        (variable_upper_bounds[idx] - center_point[idx]) / direction[idx]
    elseif direction[idx] < 0
      threshold[idx] =
        (variable_lower_bounds[idx] - center_point[idx]) / direction[idx]
    else
      # Variable doesn't move.  Rather than an infinite threshold, or a NaN if
      # the corresponding bound was infinite as well, treat it as fixed, which
      # is effectively equivalent.
      threshold[idx] = 0.0
    end
  end

  # The weighted radius squared of the indices discarded because they are below
  # the threshold.
  low_radius_sq = 0.0
  # The weighted norm squared of the objective coefficients of the indices
  # discarded because they are above the threshold.
  high_radius_sq = 0.0

  indices = collect(1:length(center_point))
  # Infinite thresholds can combine with zeros to create NaNs.  To avoid this,
  # handle indices with infinite thresholds separately.
  infinite_indices = filter(i -> isinf(threshold[i]), indices)
  high_radius_sq +=
    weighted_norm(direction[infinite_indices], norm_weights[infinite_indices])^2
  filter!(i -> isfinite(threshold[i]), indices)

  while length(indices) > 0
    test_threshold = median(threshold[indices])
    test_point =
      clamp.(
        center_point[indices] + test_threshold * direction[indices],
        variable_lower_bounds[indices],
        variable_upper_bounds[indices],
      )
    test_radius =
      weighted_norm(test_point - center_point[indices], norm_weights[indices])
    if low_radius_sq + test_radius^2 + test_threshold^2 * high_radius_sq >=
       target_radius^2
      # test_threshold is too high.  Discard indices greater than it.
      discard_indices = filter(i -> threshold[i] >= test_threshold, indices)
      high_radius_sq +=
        weighted_norm(
          direction[discard_indices],
          norm_weights[discard_indices],
        )^2
      filter!(i -> threshold[i] < test_threshold, indices)
    else
      # test_threshold is too low.  Discard indices less than it.
      discard_indices = filter(i -> threshold[i] <= test_threshold, indices)
      discard_point =
        clamp.(
          center_point[discard_indices] +
          test_threshold * direction[discard_indices],
          variable_lower_bounds[discard_indices],
          variable_upper_bounds[discard_indices],
        )
      low_radius_sq +=
        weighted_norm(
          discard_point - center_point[discard_indices],
          norm_weights[discard_indices],
        )^2
      filter!(i -> threshold[i] > test_threshold, indices)
    end
  end

  # target_threshold is the solution of
  # low_radius_sq + target_threshold^2 * high_radius_sq = target_radius^2.
  if high_radius_sq <= 0.0
    # Special case: high_radius_sq = 0.0, means all bounds hit before reaching
    # target radius.
    target_threshold = maximum(threshold)
  else
    target_threshold = sqrt((target_radius^2 - low_radius_sq) / high_radius_sq)
  end
  candidate_point =
    clamp.(
      center_point + target_threshold * direction,
      variable_lower_bounds,
      variable_upper_bounds,
    )
  return BoundConstrainedTrustRegionResult(
    candidate_point,
    dot(objective_vector, candidate_point - center_point),
  )
end

function approximately_solve_bound_constrained_trust_region(
  center_point::AbstractVector{Float64},
  objective_vector::Vector{Float64},
  variable_lower_bounds::Vector{Float64},
  variable_upper_bounds::Vector{Float64},
  norm_weights::AbstractVector{Float64},
  target_radius::Float64,
)
  direction = zeros(length(center_point))
  for idx in 1:length(center_point)
    if center_point[idx] >= variable_upper_bounds[idx] &&
       objective_vector[idx] <= 0
      continue
    end
    if center_point[idx] <= variable_lower_bounds[idx] &&
       objective_vector[idx] >= 0
      continue
    end
    direction[idx] = -objective_vector[idx] / norm_weights[idx]
  end

  direction_norm = weighted_norm(direction, norm_weights)
  if direction_norm > 0.0
    direction *= target_radius / direction_norm
  end

  return BoundConstrainedTrustRegionResult(
    center_point + direction,
    dot(objective_vector, direction),
  )
end

struct OptimalObjectiveBoundResult
  lagrangian_value::Float64
  lower_bound_value::Float64
  upper_bound_value::Float64
  "The primal solution that minimizes the localized duality gap."
  primal_solution::AbstractVector{Float64}
  "The dual solution that minimizes the localized duality gap."
  dual_solution::AbstractVector{Float64}
end

function get_gap(result::OptimalObjectiveBoundResult)
  return result.upper_bound_value - result.lower_bound_value
end

"""
The norm that defines a ball in joint primal-dual space (x,y). If
norm = MAX_NORM then max{|x|_2, |y|_2} where |.|_2 is the Euclidean norm or
the Euclidean norm in the joint space if norm = EUCLIDEAN_NORM.
"""
@enum LocalizedDualityGapNorm MAX_NORM EUCLIDEAN_NORM

"""
Given a maximum distance to optimality this function returns a tuple with
a lower bound, lagrangian value, and a upper bound on the optimal objective.

These are computed by solving
min grad_x L(primal_solution, dual_solution)^T (x - primal_solution)
- grad_y L(primal_solution, dual_solution)^T (y - dual_solution)
s.t. (x,y) in ball of radius distance_to_optimality centered
at (primal_solution,dual_solution).

See the definition of LocalizedDualityGapNorm for information about
how the ball norm is defined. Note if norm=MAX_NORM then this optimization can
be split into a seperate problems in the primal and the dual.

The lower bound is computed via
L(primal_solution, dual_solution) +
grad_x L(primal_solution, dual_solution)^T (x - primal_solution)
where x is the solution to the trust region problem.

The upper bound is computed via
L(primal_solution, dual_solution) +
grad_y L(primal_solution, dual_solution)^T (y - dual_solution)
where y is the solution to the trust region problem.
"""
function bound_optimal_objective(
  problem::QuadraticProgrammingProblem,
  primal_solution::AbstractVector{Float64},
  dual_solution::AbstractVector{Float64},
  primal_norm_weights::AbstractVector{Float64},
  dual_norm_weights::AbstractVector{Float64},
  distance_to_optimality::Float64,
  norm::LocalizedDualityGapNorm;
  solve_approximately::Bool = false,
)
  # TODO: These are also unnecessarily recomputed.
  primal_gradient =
    compute_primal_gradient(problem, primal_solution, dual_solution)

  lagrangian_value =
    compute_lagrangian_value(problem, primal_solution, dual_solution)

  dual_variable_lower_bounds = -Inf * ones(length(dual_solution))
  dual_variable_upper_bounds = Inf * ones(length(dual_solution))
  dual_variable_lower_bounds[inequality_range(problem)] .= 0.0
  dual_gradient = compute_dual_gradient(problem, primal_solution)

  if norm == MAX_NORM
    # We can split the max norm into two cases.
    # Compute the primal part
    primal_result = solve_bound_constrained_trust_region(
      primal_solution,
      primal_gradient,
      problem.variable_lower_bound,
      problem.variable_upper_bound,
      primal_norm_weights,
      distance_to_optimality,
      solve_approximately,
    )

    # Compute the dual part
    dual_result = solve_bound_constrained_trust_region(
      dual_solution,
      -dual_gradient,
      dual_variable_lower_bounds,
      dual_variable_upper_bounds,
      dual_norm_weights,
      distance_to_optimality,
      solve_approximately,
    )
    OptimalObjectiveBoundResult(
      lagrangian_value,
      # lower bound
      lagrangian_value + primal_result.value,
      # upper bound
      lagrangian_value - dual_result.value,
      # primal solution
      primal_result.solution,
      # dual solution
      dual_result.solution,
    )
  elseif norm == EUCLIDEAN_NORM
    z = [primal_solution; dual_solution]
    z_gradient = [primal_gradient; -dual_gradient]
    z_lower_bound = [problem.variable_lower_bound; dual_variable_lower_bounds]
    z_upper_bound = [problem.variable_upper_bound; dual_variable_upper_bounds]
    norm_weights = [primal_norm_weights; dual_norm_weights]

    # Compute the lower bound
    result = solve_bound_constrained_trust_region(
      z,
      z_gradient,
      z_lower_bound,
      z_upper_bound,
      norm_weights,
      distance_to_optimality,
      solve_approximately,
    )
    primal_tr_solution = result.solution[1:length(primal_solution)]
    dual_tr_solution = result.solution[(length(primal_solution)+1):end]

    return OptimalObjectiveBoundResult(
      lagrangian_value,
      # lower bound
      lagrangian_value +
      dot(primal_tr_solution - primal_solution, primal_gradient),
      # upper bound
      lagrangian_value + dot(dual_tr_solution - dual_solution, dual_gradient),
      primal_tr_solution,
      dual_tr_solution,
    )
  else
    error("unknown norm = $norm, value unknown")
  end
end


text:
# Proximal gradient methods

# Assumption: loss has gradient; penalty has prox e.g.: Lasso
# J(θ) = f(θ) + r(θ) where f is smooth
function _fit(glr::GLR, solver::ProxGrad, X, y, scratch)
    _,p,c = npc(scratch)
    c > 0 && (p *= c)
    # vector caches + eval cache
    θ   = zeros(p)   # θ_k
    Δθ  = zeros(p)   # (θ_k - θ_{k-1})
    θ̄   = zeros(p)   # θ_k + ρ Δθ // extrapolation
    ∇fθ̄ = zeros(p)
    fθ̄  = 0.0        # useful for backtracking function
    θ̂   = zeros(p)   # candidate before becoming θ_k
    # cache for extrapolation constant and stepsizes
    ω   = 0.0   # ω_k
    ω_  = 0.0   # ω_{k-1}
    ω__ = 0.0   # ω_{k-2}
    η   = 1.0   # stepsize (1/L)
    acc = ifelse(solver.accel, 1.0, 0.0) # if 0, no extrapolation (ISTA)
    # functions
    _f      = smooth_objective(glr, X, y; c=c)
    _fg!    = smooth_fg!(glr, X, y, scratch)
    _prox!  = prox!(glr)
    bt_cond = θ̂ ->
                _f(θ̂) > fθ̄ + dot(θ̂ .- θ̄, ∇fθ̄) + sum(abs2.(θ̂ .- θ̄)) / (2η)
    # loop-related
    k, tol = 1, Inf
    while k ≤ solver.max_iter && tol > solver.tol
        # --------------------------------------------------
        # This loop corresponds to the implementation of the
        # FISTA + Backtracking  algorithm in Beck & Teboulle
        # "A Fast Iterative Shrinkage Thresholding Algorithm
        # for Linear Inverse Problems" (page 193)
        # --------------------------------------------------
        # 1. linear extrapolation of past iterates
        ω   = (1.0 + sqrt(1.0 + 4.0 * ω_^2)) / 2.0
        ρ   = acc * ω__ / ω  # ω_{k-2}/ω; note that ρ != 0 only as k > 2
        θ̄  .= θ + ρ * Δθ
        # 2. attempt a prox step, modify the step until verifies condition
        fθ̄  = _fg!(∇fθ̄, θ̄)          # f and ∇f at θ̄
        _prox!(θ̂, η, θ̄ .- η .* ∇fθ̄) # candidate update
        inner = 0
        while bt_cond(θ̂) && inner < solver.max_inner
            η *= solver.beta            # shrink stepsize
            _prox!(θ̂, η, θ̄ .- η .* ∇fθ̄) # try another candidate
            inner += 1
        end
        if inner == solver.max_inner
            @warn "No appropriate stepsize found via backtracking; " *
                  "interrupting. The reason could be input data that is not standardized."
            break
        end
        # update caches
        ω__ = ω_
        ω_  = ω
        Δθ .= θ̂ .- θ
        copyto!(θ, θ̂)
        # update tolerance
        tol = norm(Δθ) / (norm(θ) + eps())
        # update niter
        k += 1
    end
    tol ≤ solver.tol || @warn "Proximal GD did not converge in " *
                              "$(solver.max_iter) iterations."
    return θ
end


text:
"""
Read-only Unix File Descriptor.
"""

@db 2 function raw_transfer(fd::FD{In}, buf, count)
    n = C.read(fd.fd, buf, count)
    @db 2 return n
end

"""
### Reading from a pseudoterminal device.

Reading from a pseudoterminal device is a special case.

On Linux `read(2)` return `EIO` when called before a client has connected.
If a client connects, writes some data and then disconnects before `read(2)`
is called: `read(2)` returns the data and then returns `EIO` if called again.

On macOS `read(2)` returns `0` when called before a client has connected.
If a client connects, writes some data and then disconnects before `read(2)`
is called: `read(2)` simply returns `0` and the data is lost.

To avoid this situation a a duplicate client fd is held open for the lifetime
of the pseudoterminal device (`fd.extra[:pt_clientfd]`). This has the effect that
`read(2)` will always return `EAGAIN` if there is no data available.
(The reader then waits for `poll(2)` to indicate when data is ready as usual.)

The remaining problem is that because `read(2)` will never return `0` there
is no way to detect when the client has closed the terminal.
This `raw_transfer` method handles this by checking if the client process
is still alive and returning `0` if it has terminated.
"""
@db 2 function raw_transfer(fd::FD{In,Pseudoterminal}, buf, count)
    n = C.read(fd.fd, buf, count)
    if n == -1
        err = errno()
        if err == C.EAGAIN && !isalive(fd.extra[:pt_client])
            @db 2 return 0 "Pseudoterminal client process died. -> HUP!"
        end
    end
    @db 2 return n
end

@db function Base.close(fd::FD{In})
    take!(fd.buffer)
    @invoke Base.close(fd::FD)
end

shutdown(fd::FD{In}) = shutdown(fd, C.SHUT_RD)

@static if isdefined(Base, :shutdown)
    Base.shutdown(fd::FD{In}) = UnixIO.shutdown(fd, C.SHUT_RD)
end

Base.isreadable(fd::FD{In}) = Base.isopen(fd)

Base.iswritable(::FD{In}) = false

Base.bytesavailable(fd::FD{In}) = bytesavailable(fd, TransferSizeMechanism(fd))

Base.bytesavailable(fd::FD, ::NoSizeMechanism) = bytesavailable(fd.buffer)

@db 1 function Base.bytesavailable(fd::FD, ::SuppoutsStatSize)
    n = bytesavailable(fd.buffer)
    pos = @cerr allow=C.EBADF C.lseek(fd, 0, C.SEEK_CUR)
    if pos != -1
        n += stat(fd).size - pos
    end
    @db 1 return n
end

@db 1 function Base.bytesavailable(fd::FD, ::SupportsFIONREAD)
    @assert bytesavailable(fd.buffer) == 0
    x = Ref(Cint(0))
    @cerr C.ioctl(fd, C.FIONREAD, x)
    @db 1 return x[]
end

Base.eof(fd::FD{In}; kw...) = eof(fd, TotalSize(fd); kw...)

Base.eof(fd::FD, ::InfiniteTotalSize; kw...) = @assert false
    # FIXME should block

Base.eof(fd::FD, ::KnownTotalSize; kw...) = bytesavailable(fd) == 0

@db 1 function Base.eof(fd::FD, ::UnknownTotalSize; kw...)
    if bytesavailable(fd.buffer) > 0
        @db 1 return false
    end
    event = @dblock fd wait(fd)                                    ;@db 1 event
    if event == C.POLLIN
        @db 1 return false
    end
    if !fd.isclosed
        Base.write(fd.buffer, readavailable(fd; kw...))
    end
    @db 1 return bytesavailable(fd.buffer) == 0
end

@db 1 function Base.unsafe_read(fd::FD{In}, buf::Ptr{UInt8}, nbytes::UInt;
                                timeout=fd.timeout)
    @require !fd.isclosed

    @with_timeout fd timeout begin
        nread = 0
        while nread < nbytes
            n = UnixIO.read(fd, buf + nread, nbytes - nread)
            if n == 0
                throw(EOFError())
            end
            nread += n
        end
        @ensure nread == nbytes
        nothing
    end
end

@db 2 function Base.read(fd::FD{In}, ::Type{UInt8}; kw...)
    @require !fd.isclosed
    eof(fd; kw...) && throw(EOFError())
    if bytesavailable(fd.buffer) == 0
        Base.write(fd.buffer, readavailable(fd; kw...))
    end
    r = Base.read(fd.buffer, UInt8)
    @db 2 return r                                "$(repr(Char(r))) $(repr(r))"
end

Base.readbytes!(fd::FD{In}, buf::Vector{UInt8}, nbytes=length(buf); kw...) =
    readbytes!(fd, buf, UInt(nbytes); kw...)

@db 1 function Base.readbytes!(fd::FD{In}, buf::Vector{UInt8}, nbytes::UInt;
                               all::Bool=true, timeout=fd.timeout)
    @require !fd.isclosed

    @with_timeout fd timeout begin
        lb::Int = length(buf)
        nread = 0
        while nread < nbytes
            @assert nread <= lb
            if (lb - nread) == 0
                lb = lb == 0 ? nbytes : min(lb * 10, nbytes)
                resize!(buf, lb)                         ;@db 1 "resize -> $lb"
            end
            @assert lb > nread
            n = UnixIO.read(fd, buf, nread + 1, lb - nread)
            if n == 0 || !all
                break
            end
            nread += n
        end
        @ensure nread <= nbytes
        @db 1 return nread
    end
end

const BUFFER_SIZE = 65536

@db 5 function Base.readavailable(fd::FD{In}; timeout=0)
    @require !fd.isclosed
    n = bytesavailable(fd)
    if n == 0
        n = BUFFER_SIZE
    end
    buf = Vector{UInt8}(undef, n)
    n = UnixIO.read(fd, buf; timeout)                          ;@db 5 n
    resize!(buf, n)
end

@db 2 function Base.readline(fd::FD{In}; timeout=fd.timeout, kw...)
    @with_timeout fd timeout begin
        @db 2 return readline(fd, ReadFragmentation(fd); kw...)
    end
end

Base.readline(fd::FD, ::ReadsBytes; kw...) = @invoke readline(fd::IO; kw...)

"""
### `readline(::UnixIO.FD{In,S_IFCHR})`

Character or Terminal devices (`S_IFCHR`) are usually used in
"canonical mode" (`ICANON`).

> In canonical mode: Input is made available line by line.
(termios(3))[https://man7.org/linux/man-pages/man3/termios.3.html].

For these devices calling `read(2)` will usually return exactly one line.
It will only ever return an incomplete line if length exceeded `MAX_CANON`.
Note that in canonical mode a line can be terminated by `CEOF` rather than
"\n", but `read(2)` does not return the `CEOF` character (e.g. when the
shell sends a "bash\$ " prompt without a newline).

This `readline` implementation is optimised for the case where `read(2)`
returns exactly one line. However it does not assume that behaviour and
will still work for devices not in canonical mode.

Notes:
 - TIOCSTI could be used to Insert a byte into the input queue
"""
@db 1 function Base.readline(fd::FD, ::ReadsLines; keep::Bool=false)
    @require !fd.isclosed
    @assert bytesavailable(fd.buffer) == 0

    linebuf = Vector{UInt8}(undef, C.MAX_CANON)
    n = UnixIO.transfer(fd, linebuf)
    if n == 0
        @db 1 return "" "EOF or timeout!"
    end

    # Check for '\n'.
    i = find_newline(linebuf, 1, n)
    @assert i == 0 || i == n

    # Trim end of line characters.
    while !keep && n > 0 && (linebuf[n] == UInt8('\r') ||
                             linebuf[n] == UInt8('\n'))
        n -= 1
    end
    resize!(linebuf, n)
    @db 1 return String(linebuf)
end

#=
@db 1 function Base.readline(fd::ReadFD{<:S_IFCHR};
                             timeout=fd.timeout,
                             wait=true,
                             keep::Bool=false)
    @require !fd.isclosed

    linebuf = Vector{UInt8}(undef, C.MAX_CANON)
    fdbuf = fd.buffer

    canonical_mode = iscanon(fd) && !(fd isa ReadFD{Pseudoterminal})

    @with_timeout fd timeout while true

        # If the fd buffer already contains a complete line, return it.
        fdbuf_n = bytesavailable(fdbuf)
        if fdbuf_n > 0
            if canonical_mode || !wait ||
               find_newline(fdbuf.data, fdbuf.ptr, fdbuf.size) != 0
                @db 1 return readline(fdbuf; keep=keep)
            end
        end

        # Read new data from fd into line buffer.
        n = UnixIO.transfer(fd, linebuf, 1, length(linebuf))
        if n == 0
            line = take!(fdbuf)
            @db 1 return String(line) "EOF or timeout!"
        end

        # Prepend old data from the fd buffer.
        if fdbuf_n > 0
            prepend!(linebuf, take!(fdbuf))
            n += fdbuf_n
        end

        # Search for end of line...
        i = find_newline(linebuf, 1, n)

        # In canonical mode, `read(2)` never returns partial lines and
        # may return lines without a newline character (if C.CEOF was sent).
        if i == 0 && (canonical_mode || !wait)
            i = n
        end

        if i != n
            # Copy everything after the newline (or after 0) into fd buffer.
            Base.write(fd.buffer, view(linebuf, i+1:n))
            @db 1 "$(repr(String(take!(copy(fd.buffer))))) -> buffer"
        end

        if i != 0
            # Found end of line.
            while !keep && i > 0 && (linebuf[i] == UInt8('\r') ||
                                     linebuf[i] == UInt8('\n'))
                i -= 1
            end
            resize!(linebuf, i)

            line = String(linebuf)
            @db 1 return line
        end
    end
end
=#

@db 2 function find_newline(buf, i, j)
    checkbounds(buf, i:j)
    while i <= j
        if @inbounds buf[i] == UInt8('\n')
            @db 2 return i
        end
        i += 1
    end
    @db 2 return 0
end

@db 2 function Base.readuntil(fd::FD{In}, d::AbstractChar;
                              timeout=fd.timeout, kw...)
    @with_timeout(fd, timeout,
                  @invoke Base.readuntil(fd::IO, d::AbstractChar; kw...))
end

@db 2 function Base.read(fd::FD{In}, n::Integer=typemax(Int);
                         timeout=fd.timeout)
    @with_timeout(fd, timeout, @invoke Base.read(fd::IO, n::Integer))
end

@db 2 function Base.read(fd::FD{In}, x::Type{String}; kw...)
    String(Base.read(fd; kw...))
end

# End of file: ReadFD.jl


text:
export NaiveMsg, NaiveCF, NaiveController

struct NaiveMsg{X,Z}
    x::X
    z::Z
end

struct NaiveCF{X,Z,U,S} <: ControllerFramework{X,Z,U,NaiveMsg{X,Z},Nothing}
    num_agents::ℕ
    central::CentralControl{U,S}
    msg_queue_len::ℕ
    ΔT::𝕋
end

NaiveCF(X, Z, N, central::CentralControl{U,S}, msg_queue_len, ΔT) where {U,S} =
    NaiveCF{X,Z,U,S}(N, central, msg_queue_len, ΔT)

struct NaiveController{X,Z,U,S} <: Controller{X,Z,U,NaiveMsg{X,Z},Nothing}
    self::ℕ
    cf::NaiveCF{X,Z,U,S}
    "Central control state."
    s_c::Ref{S}
    "Current action"
    u::Ref{U}
    t::Ref{𝕋}
    t0::𝕋
end

# NaiveCF has nothing to log
OneVision.write_logs(::NaiveController) = Dict{𝕋,Nothing}()

function OneVision.control!(
    ctrl::NaiveController{X,Z,U},
    x::X,
    z::Z,
    msgs::Each{NaiveMsg{X,Z}},
)::Tuple{U,Each{NaiveMsg{X,Z}}} where {X,Z,U}
    msgs[ctrl.self] = NaiveMsg(x, z)
    @unpack ΔT, central = ctrl.cf
    xs = [m.x for m in msgs]
    zs = [m.z for m in msgs]
    if mod(ctrl.t[] - ctrl.t0, ΔT) == 0
        ctrl.u[] = control_one(central, ctrl.s_c[], xs, zs, ctrl.t[], ctrl.self)
    end
    msgs′ = fill(NaiveMsg(x, z), length(msgs))
    ctrl.t[] += 1
    ctrl.u[], msgs′
end

function OneVision.make_controllers(
    framework::NaiveCF{X,Z,U,S},
    init_status::Each{Tuple{X,Z,U}},
    init_t::𝕋,
) where {X,Z,U,S}
    central = framework.central
    s_c = init_state(central, init_t)
    ctrls = ntuple(framework.num_agents) do i
        u = init_status[i][3]
        NaiveController(i, framework, Ref(s_c), Ref(u), Ref(init_t), init_t)
    end
    init_msg() = begin
        receives = [NaiveMsg(x, z) for (x, z, u) in init_status]
        constant_queue(receives, framework.msg_queue_len)
    end
    msgs = ntuple(_ -> init_msg(), framework.num_agents)
    ctrls, msgs
end


text:
function brown_res(x, r)
    n = length(x)
    t = x[1] - Float64(n + 1)
    for i = 2:n
        t = t + x[i]
    end
    nm1 = n - 1
    for i = 1:nm1
        r[i] = t + x[i]
    end
    t = x[1]
    for i = 2:n
        t = t * x[i]
    end
    r[n] = t - 1.0e0
    return r
end

# note for these problems m == p
function brown_jac(x, j, n)
    nm1 = n - 1
    for k = 1:n
        for i = 1:nm1
            j[i, k] = 1.0e0
            if i == k
                j[i, k] = 2.0e0
            end
        end
    end
    for k = 1:n
        t = 1.0e0
        for i = 1:n
            t =  i != k ? t * x[i] : t
        end
        j[n, k] = t
    end
    return j
end

const brown5 = let res_init=zeros(5), jac_init=zeros(5,5), x_init=fill(0.5,5)
    jac(x, jac) = brown_jac(x, jac, 5)
    f(;scale=1, verbose=false, print_steps=false) = testone("brown5", brown_res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end

const brown10 = let res_init=zeros(10), jac_init=zeros(10,10), x_init=fill(0.5,10)
    jac(x, jac) = brown_jac(x, jac, 10)
    f(;scale=1, verbose=false, print_steps=false) = testone("brown10", brown_res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end

const brown30 = let res_init=zeros(30), jac_init=zeros(30,30), x_init=fill(0.5, 30)
    jac(x, jac) = brown_jac(x, jac, 30)
    f(;scale=1, verbose=false, print_steps=false) = testone("brown30", brown_res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end

const brown40 = let res_init=zeros(40), jac_init=zeros(40,40), x_init=fill(0.5, 40)
    jac(x, jac) = brown_jac(x, jac, 40)
    f(;scale=1, verbose=false, print_steps=false) = testone("brown40", brown_res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end


text:
import MPI

function main()
    MPI.init()

    comm = MPI.COMM_WORLD

    MPI.barrier(comm)

    rank = MPI.rank(comm)
    size = MPI.size(comm)

    dst = mod(rank+1, size)
    src = mod(rank-1, size)

    N = 4

    send_mesg = Array(Float64, N)
    recv_mesg = Array(Float64, N)

    fill!(send_mesg, float64(rank))

    rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)

    println("$rank: Sending   $rank -> $dst = $send_mesg")
    sreq = MPI.Isend!(send_mesg, dst, rank+32, comm)

    stats = MPI.waitall!([rreq, sreq])

    println("$rank: Receiving $src -> $rank = $recv_mesg")

    MPI.barrier(comm)

    MPI.finalize()
end

main()


text:
# This file was generated, do not modify it. # hide
searchtime = @elapsed sknns, sdists = searchbatch(G, db, k; parallel=true)
nothing # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
"""
Module Internet

Provides functions for navigating in the internet by voice.

# Functions

###### Searching
- [`search`](@ref)

To see a description of a function type `?<functionname>`.
"""
module Internet

import DefaultApplication
import ..JustSayIt: @voiceargs, Keyboard

## CONSTANTS

const SEARCH_ENGINE = "https://google.com"

## Functions

@doc """
    internet `search`

Navigate the internet, performing one of the following actions:
    - `search`
"""
internet
@enum Action search
@voiceargs action=>(valid_input_auto=true, use_max_accuracy=true) function internet(action::Action)
    if (action == search) search_internet()
    else @info "unknown action" #NOTE: this should never happen.
    end
end

@doc "Search in internet: start search engine and automatically change to type command waiting for search keywords to be spoken."
function search_internet()
    DefaultApplication.open(SEARCH_ENGINE)
    Keyboard.type(Keyboard.words)
end

end # module Internet


text:
export GridWorldBase

using MacroTools:@forward
using Random

"""
    GridWorldBase{O} <: AbstractArray{Bool, 3}

A basic representation of grid world.
The first dimension uses multi-hot encoding to encode objects in a tile.
The second and third dimension means the height and width of the grid.
"""
struct GridWorldBase{O} <: AbstractArray{Bool, 3}
    world::BitArray{3}
    objects::O
end

get_object(w::GridWorldBase) = w.objects

function GridWorldBase(objects::Tuple{Vararg{AbstractObject}}, x::Int, y::Int)
    world = BitArray{3}(undef, length(objects), x, y)
    fill!(world, false)
    GridWorldBase(world, objects)
end

@forward GridWorldBase.world Base.size, Base.getindex, Base.setindex!

@generated function Base.to_index(::GridWorldBase{O}, x::X) where {X<:AbstractObject, O}
    i = findfirst(X .=== O.parameters)
    isnothing(i) && error("unknow object $x")
    :($i)
end

Base.setindex!(w::GridWorldBase, v::Bool, o::AbstractObject, x::Int, y::Int) = setindex!(w.world, v, Base.to_index(w, o), x, y)
Base.setindex!(w::GridWorldBase, v::Bool, o::AbstractObject, i::CartesianIndex{2}) = setindex!(w, v, o, i[1], i[2])

Base.getindex(w::GridWorldBase, o::AbstractObject, x::Int, y::Int) = getindex(w.world, Base.to_index(w, o), x, y)
Base.getindex(w::GridWorldBase, o::AbstractObject, i::CartesianIndex{2}) = getindex(w, o, i[1], i[2])
Base.getindex(w::GridWorldBase, o::AbstractObject, x::Colon, y::Colon) = getindex(w.world, Base.to_index(w, o), x, y)

#####
# utils
#####

switch!(world::GridWorldBase, x, src::CartesianIndex{2}, dest::CartesianIndex{2}) = world[x, src], world[x, dest] = world[x, dest], world[x, src]

function switch!(world::GridWorldBase, src::CartesianIndex{2}, dest::CartesianIndex{2})
    for x in axes(world, 1)
        switch!(world, x, src, dest)
    end
end

function Random.rand(f::Function, w::GridWorldBase; max_try=typemax(Int), rng=Random.GLOBAL_RNG)
    inds = CartesianIndices((size(w, 2), size(w, 3)))
    for _ in 1:max_try
        pos = rand(rng, inds)
        f(view(w, :, pos)) && return pos
    end
    @warn "a rare case happened when sampling from GridWorldBase"
    return nothing
end

#####
# get_agent_view
#####

get_agent_view_inds((i, j), (m, n), ::Left) = CartesianIndices((i-(n-1)÷2:i+(n-(n-1)÷2)-1, j-m+1:j))
get_agent_view_inds((i, j), (m, n), ::Right) = CartesianIndices((i-(n-1)÷2:i+(n-(n-1)÷2)-1, j:j+m-1))
get_agent_view_inds((i, j), (m, n), ::Up) = CartesianIndices((i-m+1:i, j-(n-1)÷2:j+(n-(n-1)÷2)-1))
get_agent_view_inds((i, j), (m, n), ::Down) = CartesianIndices((i:i+m-1, j-(n-1)÷2:j+(n-(n-1)÷2)-1))

ind_map((i,j), (m, n), ::Left) = (m-j+1, i)
ind_map((i,j), (m, n), ::Right) = (j, n-i+1)
ind_map((i,j), (m, n), ::Up) = (m-i+1, n-j+1)
ind_map((i,j), (m, n), ::Down) = (i,j)

function get_agent_view!(v::AbstractArray{Bool,3}, a::AbstractArray{Bool,3}, p::CartesianIndex, dir::LRUD)
    view_size = (size(v, 2), size(v, 3))
    grid_size = (size(a,2),size(a,3))
    inds = get_agent_view_inds(p.I, view_size, dir)
    valid_inds = CartesianIndices(grid_size)
    for ind in CartesianIndices(inds)
        if inds[ind] ∈ valid_inds
            v[:, ind_map(ind.I, view_size, dir)...] .= a[:, inds[ind]]
        end
    end
    v
end


text:
######################################################################
# gps.jl
# handles some gps functions
######################################################################

using GPS2XY

export gps_sim, gps_offset
export parse_log
export parse_belief, parse_action

# create simulation and drawing from
# j is the jammer start point (lat, lon) decimal gps coords
# x0 is the uav start point (lat, lon) decimal gps coords
function gps_sim(j::LocTuple, x0::LocTuple, actions::Vector{Pose}, observations::Vector{Float64}, m::SearchDomain, x::Vehicle, f::AbstractFilter; show_mean::Bool=false, show_cov::Bool=false)

	# always start the vehicle in the center
	x.x = m.length / 2.0
	x.y = m.length / 2.0
	dx,dy = gps_offset(x0, j)
	theta!(m, (x.x+dx, x.y+dy))

	# warn user if jammer is not where it should be
	if abs(dx) > x.x || abs(dy) > x.y
		println("WARNING: jammer outside search domain.")
	end

	# loop through all observations...
	for (oi,o) in enumerate(observations)
		update!(f, x, o)
		plot(m, f, x,show_mean=show_mean,show_cov=show_cov)
		savefig("temp_$(oi).png", format="png")
		hold(false)
		act!(m, x, actions[oi])
	end
end

function gps_sim(j::LocTuple, x0::LocTuple, actions::Vector{Pose}, observations::Vector{Float64}, b::Vector{Matrix{Float64}}, m::SearchDomain, x::Vehicle; show_mean=false, show_cov=false)

	# always start the vehicle in the center
	x.x = m.length / 2.0
	x.y = m.length / 2.0
	dx,dy = gps_offset(x0, j)
	theta!(m, (x.x+dx, x.y+dy))

	# warn user if jammer is not where it should be
	if abs(dx) > x.x || abs(dy) > x.y
		println("WARNING: jammer outside search domain.")
	end

	# loop through all observations...
	for (oi,o) in enumerate(observations)
		#update!(f, x, o)
		#plot(m, b[oi], x, show_mean=true, show_cov=true)
		plot(m, b[oi], x, show_mean=show_mean, show_cov=show_cov, obs=o)
		savefig("temp_$(oi).png", format="png")
		hold(false)
		act!(m, x, actions[oi])
	end
end

# This is not really gps related
# Reads my log file
# TODO: handle different types of beliefs
function parse_log(filename::String)

    obs = Float64[]
    actions = Pose[]
    beliefs = Array{Matrix{Float64}}(0)
	#states = Array(Pose, 0)
    states = Pose[]

	logfile = open(filename, "r")

	# create the filter and stuff
	line = readline(logfile)
	while !contains(line, "#####")
		# read header stuff

		line = readline(logfile)
	end

	# now get the actions, observations and beliefs
	while (line = readline(logfile)) != ""
		if contains(line, "observation")
			# observation line
			push!(obs, parse(Float64, split(line, " = ")[2]) )
		elseif contains(line, "action")
			# action line
			a_vec = split(split(line, " = ")[2], ",")
			# the actions we get are north,east,yaw. we want dx,dy,yaw
			dx  = parse(Float64, a_vec[2])
			dy  = parse(Float64, a_vec[1])
			dh  = parse(Float64, a_vec[3])
			push!(actions, (dx, dy, dh))
		elseif contains(line, "belief")
			# now read the next bunch of lines
			bline = readline(logfile)
			arr = split(bline, ",")
			n = length(arr)
            b = zeros(n,n)
			for i = 1:n
				b[i,n] = parse(Float64, arr[i])
			end
			for i = 2:n
				bline = readline(logfile)
				arr = split(bline, ",")
				for j = 1:n
					b[j,n-i+1] = parse(Float64, arr[j])
				end
			end
			push!(beliefs, b)
		elseif contains(line, "state")
			# action line
			s_vec = split(split(line, " = ")[2], ",")
			# the actions we get are north,east,yaw. we want dx,dy,yaw
			sx  = parse(Float64, s_vec[2])
			sy  = parse(Float64, s_vec[1])
			sh  = parse(Float64, s_vec[3])
			push!(states, (sx, sy, sh))
		end
	end
	close(logfile)

	return obs,actions,beliefs,states
end

# I'm not sure this would work...
function parse_belief(logfile)
	# now read the next bunch of lines
	bline = readline(logfile)
	if bline == ""
		error("No belief here.")
	end
	arr = split(bline, ",")
	n = length(arr)
	b = Array(Float64, n, n)
	for i = 1:n
		b[i,n] = parse(Float64, arr[i])
	end
	for i = 2:n
		bline = readline(logfile)
		arr = split(bline, ",")
		for j = 1:n
			b[j,n-i+1] = parse(Float64, arr[j])
		end
	end

	return b
end

"""
`parse_action(logfile)`

Returns tuple `(d_east,d_north,d_yaw)`
"""
function parse_action(logfile)
	line = readline(logfile)
	a_vec = split(split(line, " = ")[2], ",")
	# the actions we get are north,east,yaw. we want dx,dy,yaw
	d_east  = parse(Float64, a_vec[2])
	d_north  = parse(Float64, a_vec[1])
	d_yaw  = parse(Float64, a_vec[3])
	return d_east, d_north, d_yaw
end

# Meh, I don't know if I will use this...
function parse_obs(logfile)
	line = readline(logfile)
	if line == ""
		error("NO OBS HERE")
	end
	d_yaw  = parse(Float64, a_vec[3])
	return 1
end


text:
module CanvasWebIO

using WebIO, JSExpr

export Canvas, addmovable!, addclickable!, addstatic!

mutable struct Canvas
    w::Scope
    size
    movables::Array
    clickables::Array
    static::Array
    getter::Dict
    selected::String #actually just the field for the selection
    handler::Observable
    selection::Observable
end

function Canvas(size)
    w = Scope()
    handler = Observable(w, "handler", ["id", 0, 0])
    selection = Observable(w, "selection", "id")
    getter = Dict()
    selected = WebIO.newid("selected")
    on(selection) do val
        val
    end
    on(handler) do val
        selection[] = val[1]
        try
            getter[val[1]][] = val[2:3]
        catch
            println("Failed to assign value $(val[2:3]) to $(val[1])")
        end
    end
    Canvas(w, size, [], [], [], getter, selected, handler, selection)
end

function Canvas()
    Canvas([800,800])
end

function Base.getindex(canvas::Canvas, i)
    canvas.getter[i]
end
function (canvas::Canvas)()
    canvas_events = Dict()

    canvas_events["dragstart"]  = @js function(event)
        event.stopPropagation()
        event.preventDefault()
    end
    handler = canvas.handler
    canvas_events["drop"]  = @js function(event)
        event.preventDefault()
        event.stopPropagation()
        name = document.getElementById($(canvas.selected)).innerHTML
        #make this saner with selected as attribute, internal data? who knows
        selected_obj = document.getElementById(name)
        selected_obj.style.stroke = "none"
        #We perform the section below several times, how to remove duplicates?
        dim = selected_obj.parentElement.getBoundingClientRect()
        x = event.pageX-dim.x
        y = event.pageY-dim.y
        console.log("dropping", name, "at", x, y)
        if(selected_obj.tagName=="rect")
            xpos = x-selected_obj.getAttribute("width")/2
            ypos = y-selected_obj.getAttribute("height")/2
            selected_obj.setAttribute("x", xpos)
            selected_obj.setAttribute("y", ypos)
        elseif(selected_obj.tagName=="circle")
            xpos = x
            ypos = y
            selected_obj.setAttribute("cx", xpos)
            selected_obj.setAttribute("cy", ypos)
        end

        $handler[] = [name, xpos, ypos]
        document.getElementById($(canvas.selected)).innerHTML = ""
    end
    canvas_events["dragover"]  = @js function(event)
        console.log("dragover")
        event.preventDefault()
        event.stopPropagation()
    end
    canvas_events["click"]  = @js function(event)
        name = document.getElementById($(canvas.selected)).innerHTML
        if(name!="")
            selected_obj = document.getElementById(name)
            if selected_obj.getAttribute("draggable")=="true"
                selected_obj.style.stroke = "none"
                dim = selected_obj.parentElement.getBoundingClientRect()
                x = event.pageX-dim.x
                y = event.pageY-dim.y
                console.log("click (drop)", name, "at", x, y)
                if(selected_obj.tagName=="rect")
                    xpos = x-selected_obj.getAttribute("width")/2
                    ypos = y-selected_obj.getAttribute("height")/2
                    selected_obj.setAttribute("x", xpos)
                    selected_obj.setAttribute("y", ypos)
                elseif(selected_obj.tagName=="circle")
                    xpos = x
                    ypos = y
                    selected_obj.setAttribute("cx", xpos)
                    selected_obj.setAttribute("cy", ypos)
                end
                $handler[] = [name, xpos, ypos]
                document.getElementById($(canvas.selected)).innerHTML = ""
            end
        end

        event.preventDefault()
        event.stopPropagation()
    end

    canvas.w(dom"svg:svg[id = canvas,
        height = $(canvas.size[1]),
        width = $(canvas.size[2])]"(
                                    canvas.static...,
                                    canvas.clickables...,
                                    canvas.movables...,
                                    Node(:div, attributes=Dict("id"=>canvas.selected), ""),
                                    events = canvas_events))
end

function addclickable!(canvas::Canvas, svg::WebIO.Node)
    attr = svg.props[:attributes]
    if "id" in keys(attr)
        id = attr["id"]
    else
        id = WebIO.newid("svg")
    end
    selection = canvas.selection
    clickable_events = Dict()
    clickable_events["click"]  = @js function(event)
        name = document.getElementById($(canvas.selected)).innerHTML
        if name == this.id
            this.style.stroke = "none"
            document.getElementById($(canvas.selected)).innerHTML = ""
        else
            if name != ""
                selected_obj = document.getElementById(name)
                selected_obj.style.stroke = "none"
            end
            this.style.stroke = "green" #Change this later
            this.style.strokeWidth = 2 #Change this later
            document.getElementById($(canvas.selected)).innerHTML = this.id
            $selection[] = this.id
        end
    end
    push!(canvas.clickables,
          Node(svg.instanceof, attributes=attr, events=clickable_events))
end

function addmovable!(canvas::Canvas, svg::WebIO.Node)
    attr = svg.props[:attributes]
    if "id" in keys(attr)
        id = attr["id"]
    else
        id = WebIO.newid("svg")
    end
    if svg.instanceof.tag==:rect
        pos = Observable(canvas.w, id, parse.([attr["x"], attr["y"]]))
    elseif svg.instanceof.tag==:circle
        pos = Observable(canvas.w, id, parse.([attr["cx"], attr["cy"]]))
    end
    on(pos) do val #This ensures it's updated from the js side
        val
    end
    canvas.getter[id] = pos

    handler = canvas.handler
    attr["draggable"] = "true"
    style = Dict(:cursor => "move")
    movable_events = Dict()
    movable_events["dragstart"]  = @js function(event)
        event.stopPropagation()
        console.log("dragging", this.id)
        this.style.stroke = "red" #Change this later
        this.style.strokeWidth = 5 #Change this later
        document.getElementById($(canvas.selected)).innerHTML = this.id
    end
    movable_events["click"]  = @js function(event)
        console.log("clicking", this.id)
        name = document.getElementById($(canvas.selected)).innerHTML
        if name == ""
            this.style.stroke = "red" #Change this later
            this.style.strokeWidth = 5 #Change this later
            document.getElementById($(canvas.selected)).innerHTML = this.id
        else
            selected_obj = document.getElementById(name)
            selected_obj.style.stroke = "none"
            if selected_obj.getAttribute("draggable")=="true"
                dim = selected_obj.parentElement.getBoundingClientRect()
                x = event.pageX-dim.x
                y = event.pageY-dim.y

                console.log("click (drop)", name, "at", x, y)
                if(selected_obj.tagName=="rect")
                    xpos = x-selected_obj.getAttribute("width")/2
                    ypos = y-selected_obj.getAttribute("height")/2
                    selected_obj.setAttribute("x", xpos)
                    selected_obj.setAttribute("y", ypos)
                elseif(selected_obj.tagName=="circle")
                    xpos = x
                    ypos = y
                    selected_obj.setAttribute("cx", xpos)
                    selected_obj.setAttribute("cy", ypos)
                end
                $handler[] = [name, xpos, ypos]
            end
            document.getElementById($(canvas.selected)).innerHTML = ""
        end
    end
    push!(canvas.movables,
          Node(svg.instanceof, attributes=attr, style=style, events=movable_events))
end

function addstatic!(canvas::Canvas, svg::WebIO.Node)
    push!(canvas.static, svg)
end
end


text:
# Common Utils
# ============
#
# Common utilities shared among several algorithms (internal use only).
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioAlignments.jl/blob/master/LICENSE.md

# k: gap length
function affinegap_score(k, gap_open_penalty, gap_extend_penalty)
    return -(gap_open_penalty + gap_extend_penalty * k)
end

# Trace type
# ----------

# trace type for pairwise alignment
const Trace = UInt8

# trace bitmap
const TRACE_NONE   = 0b00000
const TRACE_MATCH  = 0b00001
const TRACE_DELETE = 0b00010
const TRACE_INSERT = 0b00100
const TRACE_EXTDEL = 0b01000
const TRACE_EXTINS = 0b10000

# Traceback
# ---------

macro start_traceback()
    esc(quote
        __alnpos = 0
        anchor_point = (i, j, __alnpos)
        op = OP_INVALID
    end)
end

# reverses the anchors sequence at the end of the traceback
# and offsets the alignment positions, so that (by default) it starts from 0
function reverse_anchors!(v::AbstractVector{AlignmentAnchor},
                          alignment_offset=!isempty(v) ? -v[end].alnpos : 0)
    r = n = length(v)
    @inbounds for i in 1:fld1(n, 2)
        vr = v[r]
        vi = v[i]
        v[i] = AlignmentAnchor(vr.seqpos, vr.refpos, vr.alnpos + alignment_offset, vr.op)
        (i != r) && (v[r] = AlignmentAnchor(vi.seqpos, vi.refpos, vi.alnpos + alignment_offset, vi.op))
        r -= 1
    end
    return v
end

macro finish_traceback()
    esc(quote
        push!(anchors, AlignmentAnchor(anchor_point..., op))
        push!(anchors, AlignmentAnchor(i, j, __alnpos, OP_START))
        reverse_anchors!(anchors)
        pop!(anchors)  # remove OP_INVALID
    end)
end

macro anchor(ex)
    esc(quote
        if op != $ex
            push!(anchors, AlignmentAnchor(anchor_point..., op))
            op = $ex
            anchor_point = (i, j, __alnpos)
        end
        __alnpos -= 1
        if ismatchop(op)
            i -= 1
            j -= 1
        elseif isinsertop(op)
            i -= 1
        elseif isdeleteop(op)
            j -= 1
        else
            @assert false
        end
    end)
end


text:
using StanModels, Distributions

# Define the Stan language model

m8_4s = "
data{
  int N;
  vector[N] y;
}
parameters{
  real sigma;
  real alpha;
  real a1;
  real a2;
}
model{
  real mu;
  mu = a1 + a2;
  y ~ normal( mu , sigma );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m8.4s", m8_4s);

# Input data for cmdstan

m8_4_data = Dict("N" => 100, "y" => rand(Normal(0, 1), 100));

# Sample using cmdstan

(sample_file, log_file) = stan_sample(sm, data=m8_4_data);

rethinking = "
        mean   sd  5.5% 94.5% n_eff Rhat
alpha 0.06 1.90 -2.22  2.49  1321    1
sigma 2.15 2.32  0.70  5.21   461    1
";

# Describe the draws

if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end