text:
using DataStructures
export interactive_evolution, interactive_evolution_timeseries

"""
    interactive_evolution(ds::DynamicalSystem, u0s; kwargs...)
Launch an interactive application that can evolve the initial conditions `u0s`
(vector of vectors) of the given dynamical system.
All initial conditions are evolved in parallel and at exactly the same time.
Two controls allow you to pause/resume the evolution and to adjust the speed.
The application can run forever (trajectories are computed on demand).

The function returns `fig, obs`. `fig` is the overarching fig
(the entire GUI) and can be recorded.
`obs` is a vector of observables, each containing the current state of the trajectory.

## Keywords
* `transform = identity` : Transformation applied to the state of the dynamical system
  before plotting. Can even return a vector that is of higher dimension than `ds`.
* `idxs = 1:min(length(transform(ds.u0)), 3)` : Which variables to plot (up to three can be chosen).
  Variables are selected after `transform` has been applied.
* `colors` : The color for each trajectory. Random colors are chosen by default.
* `lims` : A tuple of tuples (min, max) for the axis limits. If not given, they are
  automatically deduced by evolving each of `u0s` 1000 units and picking most extreme
  values (limits cannot be adjust after application is launched).
* `m = 1.0` : The trajectory endpoints have a marker. A heuristic is done to choose
  appropriate marker size given the trajectory size. `m` is a multiplier that scales
  the marker size.
* `plotkwargs = NamedTuple()` : A named tuple of keyword arguments propagated to
  the plotting function (`lines` for continuous, `scatter` for discrete systems).
  `plotkwargs` can also be a vector of named tuples, in which case each initial condition
  gets different arguments.
* `tail = 1000` : Length of plotted trajectory (in step units).
* `diffeq = NamedTuple()` : Named tuple of keyword arguments propagated to
  the solvers of DifferentialEquations.jl (for continuous systems). Because trajectories
  are not pre-computed and interpolated, it is recommended to use a combination of
  arguments that limit maximum stepsize, to ensure smooth curves. For example:
  ```julia
  using OrdinaryDiffEq
  diffeq = (alg = Tsit5(), dtmax = 0.01)
  ```
"""
function interactive_evolution(
        ds::DynamicalSystems.DynamicalSystem{IIP}, u0s;
        transform = identity, idxs = 1:min(length(transform(ds.u0)), 3),
        colors = [CYCLIC_COLORS[i] for i in 1:length(u0s)],
        tail = 1000, diffeq = NamedTuple(),
        plotkwargs = NamedTuple(), m = 1.0,
        lims = traj_lim_estimator(ds, u0s, diffeq, DynamicalSystems.SVector(idxs...), transform),
    ) where {IIP}

    @assert length(idxs) ≤ 3 "Only up to three variables can be plotted!"
    @assert length(colors) ≥ length(u0s) "You need to provide enough colors!"
    idxs = DynamicalSystems.SVector(idxs...)
    fig = Figure(resolution = (1000, 800), )
    pinteg = DynamicalSystems.parallel_integrator(ds, u0s; diffeq)
    obs, finalpoints = init_trajectory_observables(length(u0s), pinteg, tail, idxs, transform)

    # Initialize main plot with correct dimensionality
    main = fig[1,1] = init_main_trajectory_plot(
        ds, fig, idxs, lims, pinteg, colors, obs, plotkwargs, finalpoints, m
    )

    # here we define the main updating functionality
    run, sleslider = trajectory_plot_controls(fig)

    isrunning = Observable(false)
    on(run) do clicks; isrunning[] = !isrunning[]; end
    on(run) do clicks
        @async while isrunning[]
        # while isrunning[]
            DynamicalSystems.step!(pinteg)
            for i in 1:length(u0s)
                ob = obs[i]
                # last_state = iipcds ? @view(pinteg.u[:, i])[idxs] : pinteg.u[i][idxs]
                last_state = transform(DynamicalSystems.get_state(pinteg, i))[idxs]
                ob[] = push!(ob[], last_state) # push and trigger update with `=`
            end
            finalpoints[] = [x[][end] for x in obs]
            sleslider[] == 0 ? yield() : sleep(sleslider[])
            isopen(fig.scene) || break # crucial, ensures computations stop if closed window
        end
    end
    display(fig)
    fig, obs
end

function init_trajectory_observables(L, pinteg, tail, idxs, transform)
    obs = Observable[]
    T = length(idxs) == 2 ? Point2f0 : Point3f0
    for i in 1:L
        cb = CircularBuffer{T}(tail)
        fill!(cb, T(transform(DynamicalSystems.get_state(pinteg, i))[idxs]))
        push!(obs, Observable(cb))
    end
    finalpoints = Observable([x[][end] for x in obs])
    return obs, finalpoints
end

function init_main_trajectory_plot(
        ds, fig, idxs, lims, pinteg, colors, obs, plotkwargs, finalpoints, m
    )
    is3D = length(idxs) == 3
    mm = maximum(abs(x[2] - x[1]) for x in lims)
    ms = m*(is3D ? 4000 : 15)
    main = !is3D ? Axis(fig) : Axis3(fig)
    # Initialize trajectory plotted element
    for (i, ob) in enumerate(obs)
        pk = plotkwargs isa Vector ? plotkwargs[i] : plotkwargs
        if ds isa DynamicalSystems.ContinuousDynamicalSystem
            Makie.lines!(main, ob;
                color = colors[i], linewidth = 2.0, pk...
            )
        else
            Makie.scatter!(main, ob; color = colors[i],
                markersize = 2ms/3, strokewidth = 0.0, pk...
            )
        end
    end
    finalargs = if ds isa DynamicalSystems.ContinuousDynamicalSystem
        (marker = :circle, )
    else
        (marker = :diamond, )
    end
    Makie.scatter!(main, finalpoints;
        color = colors, markersize = ms, finalargs...
    )
    if !isnothing(lims)
        main.limits = lims
    end
    is3D && (main.protrusions = 50)
    return main
end
function trajectory_plot_controls(fig)
    fig[2, 1] = controllayout = GridLayout(tellwidth = false)
    run = controllayout[1, 1] = Button(fig; label = "run")
    _s, _v = 10.0 .^ (-5:0.1:0), 0.1
    pushfirst!(_s, 0.0)
    slesl = labelslider!(fig, "sleep =", _s;
    sliderkw = Dict(:startvalue => _v), valuekw = Dict(:width => 100),
    format = x -> "$(round(x; digits = 5))")
    controllayout[1, 2] = slesl.layout
    return run.clicks, slesl.slider.value
end

function traj_lim_estimator(ds, u0s, diffeq, idxs, transform)
    largedt_diffeq = merge(diffeq, (;dtmax = Inf))
    _tr = DynamicalSystems.trajectory(ds, 2000.0, u0s[1]; Δt = 1, diffeq = largedt_diffeq)
    tr = DynamicalSystems.Dataset(transform.(_tr.data))
    _mi, _ma = DynamicalSystems.minmaxima(tr)
    mi, ma = _mi[idxs], _ma[idxs]
    for i in 2:length(u0s)
        _tr = DynamicalSystems.trajectory(ds, 2000.0, u0s[i]; Δt = 1, diffeq = largedt_diffeq)
        tr = DynamicalSystems.Dataset(transform.(_tr.data))
        _mii, _maa = DynamicalSystems.minmaxima(tr)
        mii, maa = _mii[idxs], _maa[idxs]
        mi = min.(mii, mi)
        ma = max.(maa, ma)
    end
    # Alright, now we just have to put them into limits and increase by 10%
    mi = mi .- 0.05mi
    ma = ma .+ 0.05ma
    lims = [(mi[i], ma[i]) for i in 1:length(idxs)]
    lims = (lims...,)
end

"""
    interactive_evolution_timeseries(ds::DynamicalSystem, u0s, ps = nothing; kwargs...)
If `ps === nothing`, this function does the same as [`interactive_evolution`](@ref),
but in addition to the state space plot there is
a panel with the variable timeseries plotted and animated in real time.

If `ps` is not nothing, then it must be a dictionary, mapping keys
of the system parameter container (`ds.p`) to possible ranges of values. The app then will
add some additional controls on the left side which allow one to interactively change
system parameters and then click the "update" button to translate the new parameters to
system evolution. This can be done without stopping the live system evolution.
Notice that in this scenario it is recommended to provide the `lims` keyword manually.
An extra argument is returned in this case: a dictionary mapping parameter keys
to _observables_ containing their current values. You can use this to generate additional
plot elements that may depend on system parameters and thus need to be changed
if the sliders are changed.

The following additional keywords also apply:
- `total_span` : How much the x-axis of the timeseries plots should span (in real time units)
- `linekwargs = NamedTuple()` : Extra keywords propagated to the timeseries plots.
- `pnames = Dict(keys(ps) .=> keys(ps))` : Dictionary mapping parameter keys to labels.
  Only valid if `ps` is a dictionary and not `nothing`.
"""
function interactive_evolution_timeseries(
        ds::DynamicalSystems.DynamicalSystem{IIP}, u0s, ps = nothing;
        transform = identity, idxs = 1:min(length(transform(ds.u0)), 3),
        colors = [CYCLIC_COLORS[i] for i in 1:length(u0s)],
        tail = 1000, diffeq = NamedTuple(),
        plotkwargs = NamedTuple(), m = 1.0,
        lims = traj_lim_estimator(ds, u0s, diffeq, DynamicalSystems.SVector(idxs...), transform),
        total_span = ds isa DynamicalSystems.ContinuousDynamicalSystem ? 10 : 50,
        linekwargs = ds isa DynamicalSystems.ContinuousDynamicalSystem ? (linewidth = 4,) : (),
        pnames = isnothing(ps) ? nothing : Dict(keys(ps) .=> keys(ps)),
    ) where {IIP}

    # TODO: Separate into traj_idxs and ts_idxs, in case one wants to plot
    # more timeseries than the ones in the state space plot

    N = length(u0s)
    @assert length(idxs) ≤ 3 "Only up to three variables can be plotted!"
    @assert length(colors) ≥ length(u0s) "You need to provide enough colors!"
    idxs = DynamicalSystems.SVector(idxs...)
    fig = Figure(resolution = (1600, 800), )
    pinteg = DynamicalSystems.parallel_integrator(ds, u0s; diffeq)
    obs, finalpoints = init_trajectory_observables(length(u0s), pinteg, tail, idxs, transform)
    tdt = total_span/50

    # Initialize main plot with correct dimensionality
    main = fig[1,1] = init_main_trajectory_plot(
        ds, fig, idxs, lims, pinteg, colors, obs, plotkwargs, finalpoints, m
    )

    # Initialize timeseries data:
    allts = [] # each entry is one timeseries plot, which contains N series
    for i in 1:length(idxs)
        individual_ts = Observable[]
        for j in 1:N
            cb = CircularBuffer{Point2f0}(tail)
            fill!(cb, Point2f0(
                pinteg.t, transform(DynamicalSystems.get_state(pinteg, j))[idxs][i])
            )
            push!(individual_ts, Observable(cb))
        end
        push!(allts, individual_ts)
    end

    # here we define the main updating functionality
    run, sleslider = trajectory_plot_controls(fig)

    # Initialize timeseries plots:
    ts_axes = []
    for i in 1:length(idxs)
        ax = fig[:, 2][i, 1] = Axis(fig; xticks = LinearTicks(5))
        push!(ts_axes, ax)
        individual_ts = allts[i]
        for j in 1:N
            lines!(ax, individual_ts[j]; color = colors[j], linekwargs...)
            if ds isa DynamicalSystems.DiscreteDynamicalSystem
                scatter!(ax, individual_ts[j]; color = colors[j])
            end
        end
        tight_xticklabel_spacing!(ax)
        ax.ylabel = string(('x':'z')[i])
        ylims!(ax, lims[i])
    end
    for i in 1:length(idxs)-1
        linkxaxes!(ts_axes[i], ts_axes[end])
        hidexdecorations!(ts_axes[i], grid = false)
    end
    xlims!(ts_axes[end], pinteg.t - tdt, total_span+tdt)

    # Functionality of live evolution
    isrunning = Observable(false)
    on(run) do clicks; isrunning[] = !isrunning[]; end
    on(run) do clicks
        @async while isrunning[]
        # while isrunning[]
            DynamicalSystems.step!(pinteg)
            for i in 1:N
                ob = obs[i]
                last_state = transform(DynamicalSystems.get_state(pinteg, i))[idxs]
                pushupdate!(ob, last_state) # push and trigger update with `=`
                for k in 1:length(idxs)
                    pushupdate!(allts[k][i], Point2f0(pinteg.t, last_state[k]))
                end
            end
            finalpoints[] = [x[][end] for x in obs]
            sleslider[] == 0 ? yield() : sleep(sleslider[])
            t_current = pinteg.t
            t_prev = max(0, t_current - total_span)
            xlims!(ts_axes[end], t_prev-tdt, max(t_current, total_span)+tdt)
            isopen(fig.scene) || break # crucial, ensures computations stop if closed window
        end
    end

    # Live parameter changing
    if !isnothing(ps)
        playout = fig[:, 3] = GridLayout(tellheight = false)
        slidervals, returnvals = add_param_controls!(fig, playout, ps, copy(ds.p), pnames)
        update = Button(fig, label = "update", tellwidth = false)
        playout[length(ps)+1, :] = update
        on(update.clicks) do clicks
            update_ds_parameters!(ds, slidervals, returnvals)
        end
    else
        returnvals = nothing
    end

    display(fig)
    return fig, obs, returnvals
end

function add_param_controls!(fig, playout, ps, p0, pnames)
    slidervals = Dict{keytype(ps), Observable}()
    returnvals = Dict{keytype(ps), Observable}()
    for (i, (l, vals)) in enumerate(ps)
        startvalue = p0[l]
        label = string(pnames[l])
        sll = labelslider!(fig, label, vals; sliderkw = Dict(:startvalue => startvalue))
        slidervals[l] = sll.slider.value # directly add the observable
        returnvals[l] = Observable(sll.slider.value[]) # will only get updated on button
        playout[i, :] = sll.layout
    end
    return slidervals, returnvals
end

function update_ds_parameters!(ds, slidervals, returnvals)
    for l in keys(slidervals)
        v = slidervals[l][]
        returnvals[l][] = v
        DynamicalSystems.set_parameter!(ds, l, v)
    end
end


text:
"""
    PRN{name}

Package Registry Name
"""
struct PRN{name} end

"""
    PRN(name::String)

Create a `PRN` (Pacakge Registry Name) object.
"""
PRN(name::String) = PRN{Symbol(name)}()

macro PRN_str(name::String)
    return PRN{Symbol(name)}
end

Base.show(io::IO, ::PRN{registry}) where {registry} = print(io, "Pacakge Registry ", string(registry))

Base.@kwdef struct VersionTokens
    major::String = "major"
    minor::String = "minor"
    patch::String = "patch"
end

const VERSION_TOKENS = VersionTokens()
Base.show(io::IO, vt::VersionTokens) = print(io, "(", vt.major, ", ", vt.minor, ", ", vt.patch, ")")

Base.in(version::String, tokens::VersionTokens) = (version == tokens.major) ||
    (version == tokens.minor) || (version == tokens.patch)

function is_version_number(version)
    occursin(r"[0-9]+.[0-9]+.[0-9]+", version) ||
        occursin(r"v[0-9]+.[0-9]+.[0-9]+", version)
end

struct Project
    path::String
    toml::String
    pkg::Pkg.Types.Project
    git::Cmd
    branch::String
    quiet::Bool
end

function Project(path::String=pwd(); gitconfig=Dict(), branch="master", quiet=false)
    toml = Base.current_project(path)
    toml === nothing && error("cannot find (Julia)Project.toml in $path")
    path = dirname(toml)
    pkg = Pkg.Types.read_project(toml)
    git = RegistryTools.gitcmd(path, gitconfig)
    return Project(path, toml, pkg, git, branch, quiet)
end

Base.show(io::IO, p::Project) = print(io, "Project(", p.path, ")")

function current_branch(p::Project)
    return readchomp(`$(p.git) rev-parse --abbrev-ref HEAD`)
end

function commit_toml(project::Project; push::Bool=false)
    git = project.git
    version_number = project.pkg.version
    run(`$git add $(project.toml)`)
    run(`$git commit -m"bump version to $version_number"`)

    if push
        run(`$git push origin $(project.branch)`)
    end
end

function reset_last_commit(project::Project; push=false)
    git = project.git
    run(`$git revert --no-edit --no-commit HEAD`)
    run(`$git commit -m "revert version bump due to an error occured in IonCLI"`)
    if push
        run(`$git push origin $(project.branch)`)
    end
end

function checkout(f, p::Project)
    old = current_branch(p)

    if old != p.branch
        @info "checking out to $(p.branch)"
        run(`$(p.git) checkout $(p.branch)`)
    end

    f()

    if old != p.branch
        run(`$(p.git) checkout $old`)
    end
end

"""
release a package.

# Arguments

- `version`: version number you want to release. Can be a specific version, "current" or either of $(VERSION_TOKENS)
- `path`: path to the project you want to release.

# Options

- `-r,--registry <registry name>`: registry you want to register the package.
    If the package has not been registered, ion will try to register
    the package in the General registry. Or the user needs to specify
    the registry to register using this option.

- `-b, --branch <branch name>`: branch you want to register.

# Flags

- `-q,--quiet`: do not promote anything.
"""
@cast function release(version::String, path::String=pwd(); registry="", branch="master", quiet::Bool=false)
    project = Project(path; branch=branch)
    # new version needs to be pushed
    # so the JuliaRegistrator can find
    # it later
    checkout(project) do
        if LocalRegistry.is_dirty(project.path, Dict())
            error("package repository is dirty, please commit or stash changes.")
        end

        if version != "current"
            update_version!(project, version)
            commit_toml(project; push=true)
        end

        try
            register(registry, project)
        catch e
            if version != "current"
                reset_last_commit(project; push=true)
            end
            rethrow(e)
        end
    end
    return
end

function register(registry::String, project::Project)
    if isempty(registry) # registered package
        path = LocalRegistry.find_registry_path(nothing, project.pkg)
        return register(PRN(basename(path)), project)
    else
        return register(PRN(registry), project)
    end
end

function register(::PRN"General", project::Project)
    github_token = read_auth()
    auth = GitHub.authenticate(github_token)
    HEAD = read_head(project.git)
    comment_json = Dict{String, Any}(
        "body" => registrator_msg(project),
    )

    repo = github_repo(project.git)
    if repo === nothing
        error("not a GitHub repository")
    end

    comment = GitHub.create_comment(repo, HEAD, :commit; params=comment_json, auth=auth)
    println(comment)
    return
end

function register(registry::PRN, project::Project)
    error("register workflow is not defined for $registry")
end

function registrator_msg(project)
    msg = "Released via [Ion CLI](https://github.com/Roger-luo/IonCLI.jl)\n"
    msg *= "@JuliaRegistrator register"
    if project.branch == "master"
        return msg
    else
        return msg * " branch=$(project.branch)"
    end
end

function read_auth()
    for key in ENV_TOKEN_NAMES
        if haskey(ENV, key)
            return ENV[key]
        end
    end

    buf = Base.getpass("GitHub Access Token (https://github.com/settings/tokens)")
    auth = read(buf, String)
    Base.shred!(buf)
    return auth
end

function read_head(git, branch="master")
    return readchomp(`$git rev-parse --verify HEAD`)
end

function read_remote_push(git, remote="origin")
    return readchomp(`$git config --get remote.$remote.url`)
end

function github_repo(git, remote="origin")
    url = read_remote_push(git, remote)
    github_https = "https://github.com/"
    github_ssh = "git@github.com:"
    if startswith(url, github_https)
        if endswith(url, ".git")
            return url[length(github_https)+1:end-4]
        else
            return url[length(github_https)+1:end]
        end
    elseif startswith(url, github_ssh)
        return url[length(github_ssh)+1:end-4]
    else
        return
    end
end

function update_version!(project::Project, version)
    if is_version_number(version)
        version_number = VersionNumber(version)
    elseif version in VERSION_TOKENS
        version_number = bump_version(project, version)
    else
        error("invalid version $version")
    end

    if !project.quiet
        latest_version = find_max_version(project.pkg.name)

        if latest_version === nothing
            println("package not found in local registries")
        else
            println("latest registered version: ", latest_version)
            if latest_version > version_number
                @warn "input version is smaller than registered version"
            end
        end

        println(" "^10, "current version: ", project.pkg.version)
        println(" "^7, "version to release: ", version_number)
        if !prompt("do you want to update Project.toml?")
            exit(0)
        end
    end

    write_version(project, version_number)
    return project
end

function bump_version(project::Project, token::String)
    if project.pkg.version === nothing
        return bump_version(v"0.0.0", token)
    else
        return bump_version(project.pkg.version, token)
    end
end

function bump_version(version::VersionNumber, token::String)
    if token == VERSION_TOKENS.major
        return VersionNumber(version.major+1, 0, 0)
    elseif token == VERSION_TOKENS.minor
        return VersionNumber(version.major, version.minor+1, 0)
    elseif token == VERSION_TOKENS.patch
        return VersionNumber(version.major, version.minor, version.patch+1)
    else
        error("invalid token $token")
    end
end

function write_version(project::Project, version::VersionNumber)
    project.pkg.version = version
    open(project.toml, "w+") do f
        Pkg.TOML.print(f, to_dict(project))
    end
end

to_dict(p::Project) = to_dict(p.pkg)

function to_dict(project::Pkg.Types.Project)
    project_keys = [:name, :uuid, :authors, :version, :deps, :compat, :extras, :targets]
    t = []
    for key in project_keys
        push_maybe!(t, project, key)
    end

    # write other part back
    for key in keys(project.other)
        if !(key in string.(project_keys))
            push!(t, project.other[key])
        end
    end
    return OrderedDict(t)
end

function push_maybe!(t::Vector, project::Pkg.Types.Project, name::Symbol)
    key = string(name)
    if hasfield(Pkg.Types.Project, name)
        member = getfield(project, name)
        if member !== nothing
            push!(t, key => member)
        end
    else
        member = get(project.other, key, nothing)
        if member !== nothing
            push!(t, key => member)
        end
    end
    return t
end

ext: jl


text:
padding_spec =
[
	(:top, Number, 80),
	(:left, Number, 80),
	(:bottom, Number, 80),
	(:right, Number, 80)
]

primitivefactory(:VegaPadding, padding_spec)

max_stars_repo_name: jpthor/vega


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAppsV1beta1ScaleSpec <: SwaggerModel
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas

    function IoK8sApiAppsV1beta1ScaleSpec(;replicas=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta1ScaleSpec, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        o
    end
end # type IoK8sApiAppsV1beta1ScaleSpec

const _property_map_IoK8sApiAppsV1beta1ScaleSpec = Dict{Symbol,Symbol}(Symbol("replicas")=>Symbol("replicas"))
const _property_types_IoK8sApiAppsV1beta1ScaleSpec = Dict{Symbol,String}(Symbol("replicas")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAppsV1beta1ScaleSpec }) = collect(keys(_property_map_IoK8sApiAppsV1beta1ScaleSpec))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta1ScaleSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAppsV1beta1ScaleSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta1ScaleSpec }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta1ScaleSpec[property_name]

function check_required(o::IoK8sApiAppsV1beta1ScaleSpec)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta1ScaleSpec }, name::Symbol, val)
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# A basic alias analysis with limited scope:
#
# 1. Only objects pointed to by variables, but not elements of arrays or other structs
# 2. We only consider variables that are assigned once
# 3. Very limited inter-precedural (i.e. function call) analysis
#
# The only useful result from this alias analysis is whether a variable definitely
# doesn't alias with anything else, or its aliases are no longer alive.
#
# We are NOT interested in the set of potentially aliased variables.
#
# The algorithm is basically an abstract interpreter of Julia AST.

module AliasAnalysis

import ..DebugMsg
DebugMsg.init()

using CompilerTools.LambdaHandling
using CompilerTools
using CompilerTools.Helper

# state to keep track of variable values
const Unknown  = -1
const NotArray = 0

type State
  linfo  :: LambdaVarInfo
  baseID :: Int
  locals :: Dict{LHSVar, Int}
  revmap :: Dict{Int, Set{LHSVar}}
  nest_level :: Int
  top_level_idx :: Int
  liveness :: CompilerTools.LivenessAnalysis.BlockLiveness
  noReAssign :: Bool
end

init_state(linfo, liveness, noReAssign) = State(linfo, 0, Dict{LHSVar,Int}(), Dict{Int, Set{LHSVar}}(), 0, 0, liveness, noReAssign)

function increaseNestLevel(state)
state.nest_level = state.nest_level + 1
end

function decreaseNestLevel(state)
state.nest_level = state.nest_level -1
end

function next_node(state)
  local n = state.baseID + 1
  state.baseID = n
  return n
end

function update_node(state, v, w)
  if !haskey(state.locals, v)
    # new initialization
    state.locals[v] = w
    if haskey(state.revmap, w)
      push!(state.revmap[w], v)
    else
      state.revmap[w] = push!(Set{LHSVar}(), v)
    end
  else
    # when a variable is initialized more than once, set to Unknown
    state.locals[v] = Unknown
    if haskey(state.revmap, w)
      for u in state.revmap[w]
        state.locals[u] = Unknown
      end
      pop!(state.revmap, w)
    end
  end
end

function update_unknown(state, v)
  state.locals[v] = Unknown
end

function update_notarray(state, v)
  state.locals[v] = NotArray
end

function update(state, v, w)
  if w > 0
    update_node(state, v, w)
  elseif w == NotArray
    update_notarray(state, v)
  else
    update_unknown(state, v)
  end
end

function lookup(state, v)
  if haskey(state.locals, v)
    state.locals[v]
  else
    Unknown
  end
end

# (:lambda, {param, meta@{localvars, types, freevars}, body})
function from_innerlambda(state, expr :: Expr)
  local head = expr.head
  local ast  = expr.args
  local typ  = expr.typ
  assert(length(ast) == 3)
  linfo, body = lambdaToLambdaVarInfo(expr)
  # very conservative handling by setting free variables to Unknown.
  # TODO: may want to simulate function call at call site to get
  #       more accurate information.
  for (v, vd) in linfo.escaping_defs
    update_unknown(state, v)
  end
  return NotArray
end

function from_exprs(state, ast, callback=not_handled, cbdata :: ANY = nothing)
  local len  = length(ast)
  [ from_expr(state, exp, callback, cbdata) for exp in ast ]
end

function from_body(state, expr :: Expr, callback, cbdata :: ANY)
  exprs = expr.args
  local ret = NotArray       # default return
  for i = 1:length(expr.args)
    if state.nest_level == 0
      state.top_level_idx = i
    end
    ret = from_expr(state, exprs[i], callback, cbdata)
  end
  return ret
end

function from_assignment(state, expr :: Expr, callback, cbdata :: ANY)
  local head = expr.head
  local ast  = expr.args
  local typ  = expr.typ
  assert(length(ast) == 2)
  local lhs = ast[1]
  local rhs = ast[2]
  @dprintln(2, "AA ", lhs, " = ", rhs)
  lhs = toLHSVar(lhs)
  assert(isa(lhs, LHSVar))
  if lookup(state, lhs) != NotArray
    rhs = from_expr(state, rhs, callback, cbdata)
    # if all vars that have rhs are not alive afterwards
    # then we can safely give v a fresh ID.
    @dprintln(3, "from_assignment rhs after from_expr = ", rhs)
    # Optionally allow re-assignment to be considered as non-aliasing when RHS is dead.
    if state.noReAssign == false && state.nest_level == 0
      tls = CompilerTools.LivenessAnalysis.find_top_number(state.top_level_idx, state.liveness)
      if tls == nothing
        @dprintln(0, "state.liveness = ", state.liveness)
        @dprintln(0, "state.top_level_idx = ", state.top_level_idx)
      else
        @dprintln(3, "tls = ", tls)
      end
      assert(tls != nothing)
      assert(CompilerTools.LivenessAnalysis.isDef(lhs, tls))
      if (haskey(state.revmap, rhs))
        @dprintln(3, "rhs in state.revmap")
        dead = true
        for v in state.revmap[rhs]
          @dprintln(4, "dead = ", dead, " v = ", v, " live_out = ", tls.live_out)
          dead = dead && !in(v, tls.live_out)
        end
        if dead
          @dprintln(3, "dead is true so calling next_node")
          rhs = next_node(state)
        end
      else
        @dprintln(3, "rhs not in state.revmap = ", state.revmap)
      end
    end
    @dprintln(2, "AA update ", lhs, "::", typeof(lhs), " <- ", rhs)
    update(state, lhs, rhs)
  end
end

function from_call(state, expr :: Expr, callback, cbdata)
  # The assumption here is that the program has already been translated
  # by DomainIR, and all args are either SymbolNode or Constant.
  local head = expr.head
  local ast = expr.args
  local typ = expr.typ
  assert(length(ast) >= 1)
  if head == :invoke
    # skip the first argument that is LambdaInfo
    ast = ast[2:end]
  end
  local fun  = ast[1]
  local args = ast[2:end]
  @dprintln(2, "AA from_call: fun=", fun, " typeof(fun)=", typeof(fun), " args=",args, " typ=", typ)
  #fun = from_expr(state, fun)
  #@dprintln(2, "AA from_call: new fun=", fun)
  if isa(fun, RHSVar)
    ftyp = getType(fun, state.linfo)
    @dprintln(2, "AA from_call: ftyp=", ftyp)
    if ftyp <: Function
       fun = ftyp.name.name
    end
  end
  @dprintln(2, "AA from_call: fun=", fun)
  fun = isa(fun, TopNode) ? fun.name : fun
  fun = isa(fun, GlobalRef) ? fun.name : fun
  if is(fun, :reshape)
    # assume reshape's return result to always alias its first argument
    return from_expr(state, args[1], callback, cbdata)
  elseif is(fun, :arrayref) || is(fun, :arrayset) || is(fun, :getindex) || is(fun, :setindex!)
    # This is actually a conservative answer since arrayref might return
    # an array too, but we don't consider it as a case to handle.
    return Unknown
  elseif fun == :ccall && (args[1] == QuoteNode(:jl_new_array) || args[1] == QuoteNode(:jl_alloc_array_1d) ||
         args[1] == QuoteNode(:jl_alloc_array_2d) || args[1] == QuoteNode(:jl_alloc_array_3d))
    return next_node(state)
  elseif isFromBase(state, ast[1]) && endswith(string(fun), '!')
    # if the function is from base, and follow julia convention of ending with !, and
    # only one input is an array, then we consider it aliases the output
    alias = nothing
    for exp in args
      if (isa(exp, Expr) && (exp.typ == nothing || !iselementarytype(exp.typ)))
        # non-flattened input AST is not handled
        alias = Unknown
      elseif isa(exp, RHSVar)
        # A local var
        typ = getType(exp, state.linfo)
        if !iselementarytype(typ)
          # complicated type
          if isArrayType(typ)
            # an array type!
            if !is(alias, nothing) # more than one array as input
              alias = Unknown
            else
              # remember possible match
              alias = lookup(state, toLHSVar(exp))
            end
          else
            # non array, non elementary type
            alias = Unknown
          end
        end
      end
    end
    alias = is(alias, nothing) ? Unknown : alias
    if alias == Unknown
      # if we don't know the result, it could alias any input
      for exp in args
        if isa(exp, RHSVar)
            update_unknown(state, toLHSVar(exp))
        end
      end
    end
    return alias
  else
    @dprintln(2, "AA: unknown call ", fun)
    # For unknown calls, conservative assumption is that the result
    # might alias one of (or an element of) the non-leaf-type inputs.
    for exp in args
        if isa(exp, RHSVar)
            update_unknown(state, toLHSVar(exp))
        end
    end
    return Unknown
  end
end

function from_return(state, expr :: Expr, callback, cbdata :: ANY)
  local head = expr.head
  local typ  = expr.typ
  local args = from_exprs(state, expr.args, callback, cbdata)
  if length(args) == 1
    return args[1]
  else
    return Unknown
  end
end

function process_callback(handled::Array, state, ast, callback, cbdata)
    @dprintln(3,"Processing expression from callback for ", ast)
    @dprintln(3,handled)
    return from_exprs(state, handled, callback, cbdata)
end

# AST node replaced
function process_callback(handled::Expr, state, ast, callback, cbdata)
    return from_expr_inner(state, handled, callback, cbdata)
end

function process_callback(handled::Integer, state, ast, callback, cbdata)
    return handled
end

function process_callback(handled::Any, state, ast, callback, cbdata)
    return from_expr_inner(state, ast, callback, cbdata)
end

function from_expr_inner(state, ast::Expr, callback, cbdata)

    local head = ast.head
    local args = ast.args
    local typ  = ast.typ
    @dprintln(2, " --> ", head)
    if is(head, :lambda)
        return from_innerlambda(state, ast)
    elseif is(head, :body)
        return from_body(state, ast.args, callback, cbdata)
    elseif is(head, :(=))
        return from_assignment(state, ast, callback, cbdata)
    elseif is(head, :return)
        return from_return(state, ast, callback, cbdata)
    elseif is(head, :invoke)
        return from_call(state, ast, callback, cbdata)
    elseif is(head, :call)
        return from_call(state, ast, callback, cbdata)
        # TODO: catch domain IR result here
    elseif is(head, :call1)
      return from_call(state, ast, callback, cbdata)
    elseif is(head, :method)
        # skip
    elseif is(head, :line)
        # skip
    elseif is(head, :new)
        # skip
    elseif is(head, :boundscheck) || is(head, :inbounds)
        # skip?
    elseif is(head, :type_goto)
        # skip?
    elseif is(head, :gotoifnot)
        # skip
    elseif is(head, :loophead)
        # skip
    elseif is(head, :loopend)
        # skip
    elseif is(head, :meta)
        # skip
    elseif is(head, :static_parameter)
        # skip
    else
        throw(string("from_expr: unknown Expr head :", head))
    end

  return Unknown
end

function from_expr_inner(state, ast::RHSVar, callback, cbdata)
    @dprintln(2, "RHSVar", ast, " :: ", typeof(ast))
    return lookup(state, toLHSVar(ast))
end

function from_expr_inner(state, ast::ANY, callback, cbdata)
  @dprintln(2, " not handled ", ast)
  return Unknown
end

function not_handled(a,b,c)
  nothing
end

function from_expr(state, ast::LambdaInfo, callback=not_handled, cbdata :: ANY = nothing)
    return from_expr(state, getBody(ast), callback, cbdata)
end

function from_expr(state, ast :: ANY, callback=not_handled, cbdata :: ANY = nothing)
  # "nothing" output means couldn't be handled
  handled = callback(ast, state, cbdata)
  return process_callback(handled, state, ast, callback, cbdata);
end

function isFromBase(state, x::GlobalRef)
  startswith(string(Base.resolve(x, force=true).mod), "Base")
end

function isFromBase(state, x::RHSVar)
  tname = getType(x, state.linfo).name
  isFromBase(state, GlobalRef(tname.module, tname.name))
end

function isFromBase(state, x::ANY)
  false
end

function iselementarytype(typ::DataType)
  (typ.name == Type.name || eltype(typ) == typ)
end

function iselementarytype(typ::Any)
  return false
end

function from_lambda(LambdaVarInfo :: LambdaVarInfo, body, liveness, callback=not_handled, cbdata :: ANY = nothing; noReAssign = false)
  local state = init_state(LambdaVarInfo, liveness, noReAssign)
  #@dprintln(2, "AA ", isa(body, Expr), " ", is(body.head, :body))
  for v in getLocalVariables(LambdaVarInfo)
    vtyp = getType(v, LambdaVarInfo)
    @dprintln(2, "v = ", v, " vtyp = ", vtyp)
    if !isArrayType(vtyp)
      update_notarray(state, v)
    elseif isInputParameter(v, LambdaVarInfo)
      update_unknown(state, v)
    end
  end
  @dprintln(2, "AA locals=", state.locals)
  from_body(state, body, callback, cbdata)
  @dprintln(2, "AA locals=", state.locals)
  local revmap = Dict{Int, LHSVar}()
  local unique = Set{LHSVar}()
  # keep only variables that have unique object IDs.
  # TODO: should consider liveness either here or during analysis,
  #       since its ok to alias dead vars.
  for (v, w) in state.locals
    if w > 0
      if haskey(revmap, w)
        delete!(unique, revmap[w])
      else
        push!(unique, v)
        revmap[w] = v
      end
    end
  end
  @dprintln(2, "AA after alias analysis: ", unique)
  # return the set of variables that are confirmed to have no aliasing
  return unique
end

function from_lambda(expr, liveness, callback=not_handled, cbdata :: ANY = nothing; noReAssign = false)
  LambdaVarInfo, body = lambdaToLambdaVarInfo(expr)
  from_lambda(LambdaVarInfo, body, liveness, callback, cbdata; noReAssign = noReAssign)
end

end


text:
using Dates, HTTP
"""
    latest_beagle()
---
# Description
- Return the latest beagle URL.
################################################################################
Note Beagle is not version controlled by git.  It seems that the names for new
versions are following the pattern like `beagle.ddMmmyy.xxx.jar`.  So I just
grab all the file names in https://faculty.washington.edu/browning/beagle/.
There must also be better ways to get the URL.
"""
function latest_beagle()
    later(a, b) = begin
        mon = Dict("Jan" => 1, "Feb" => 2, "Mar" => 3, "Apr" => 4,
                   "May" => 5, "Jun" => 6, "Jul" => 7, "Aug" => 8,
                   "Sep" => 9, "Oct" => 10, "Nov" => 11, "Dec" => 12, "mmm" => 0)
        da, db, ya, yb = parse.(Int, [a[8:9], b[8:9], a[13:14], b[13:14]])
        ma, mb = mon[a[10:12]], mon[b[10:12]]
        x = ya * 10000 + ma * 100 + da
        y = yb * 10000 + mb * 100 + db
        return x > y ? a : b
    end

    BeagleURL = "https://faculty.washington.edu/browning/beagle/"
    beagle = "beagle.00mmm00.xxx.jar"
    web = HTTP.request("GET", BeagleURL)
    txt = split(String(web.body), "\n")
    for line in txt
        r = match(r"beagle\.\d{2}\w{3}\d{2}.{5}jar", line)
        r == nothing || (beagle = later(beagle, r.match))
    end
    message("The current Beagle is $beagle.")
    return joinpath(BeagleURL, beagle)
end

"""
    latest_plink()
---
Return URL of the latest plink version 1 for Linux x86_64.
"""
function update_plink()
    message("Updating plink to the latest")
    plinkURL = "http://s3.amazonaws.com/plink1-assets"
    latest = "plink_linux_x86_64_latest.zip"
    odir = pwd()
    cd(abg_bin)
    download(joinpath(plinkURL, latest), latest)
    run(`unzip -u $latest`)
    cd(odir)
end

"""
    make_bins()
---
This is to make binaries from my C++ files.
It is supposed to run at the startup.
But this can also be run after new C++ codes were written.
"""
function make_bins()
    subtitle("C++ binaries")
    isdir(abg_bin) || mkdir(abg_bin)

    message("Compiling ABG C++ codes in `cpp` into `bin`")
    # free-space is removed, as it uses some new c++ features.
    width = 70
    for cpp in readdir(abg_cpp)
        bin = cpp[1:end-4]
        if (!isfile("$abg_bin/$bin")) || (stat("$abg_bin/$bin").mtime < stat("$abg_cpp/$bin.cpp").mtime)
            print(lpad("g++ -O3 -Wall -std=c++17 cpp/$cpp -o bin/$bin", width))
            run(`g++ -O3 -Wall -std=c++17 -o $abg_bin/$bin $abg_cpp/$cpp`)
            done()
        else
            print(lpad("bin/$bin", width))
            done("OK")
        end
    end
end

"""
    Update()
---
Every time one starts this package, the package will automatically check if `plink` and `beagle.jar` exist.
If not, the package will download the latest version of these two files.
The package will update, or download them on 17th every month anyway.
"""
function Update(force = false)
    isdir(abg_bin) || mkdir(abg_bin)
    if Dates.day(now()) == 17 || force
        rm(beagle, force=true)
        rm(plink,  force=true)
    end

    subtitle("Updating the binaries")

    item("Beagle related")
    beagle2vcfURL = "https://faculty.washington.edu/browning/beagle_utilities/"
    beagle2vcf = "beagle2vcf.jar"

    isfile(beagle) || download(latest_beagle(), beagle)
    isfile(joinpath(abg_bin, beagle2vcf)) || download(joinpath(beagle2vcfURL, beagle2vcf), joinpath(abg_bin, beagle2vcf))
    done()

    item("plink")
    isfile(plink) || update_plink()
    done()

    item("C++ binaries")
    make_bins()
end


text:
"""
    ReversedEasing

Will be used to reverse an easing inside [`easing_to_animation`](@ref).
Can be constructed from an easing function using [`rev`](@ref).
"""
struct ReversedEasing
    easing::Easing
end

"""
    rev(e::Easing)

Reverse an easing function such that `easing_to_animation` maps it to `[1.0, 0.0]` instead of `[0.0, 1.0]`.
An example can be seen in [`rotate`](@ref)
"""
rev(e::Easing) = ReversedEasing(e)

"""
    easing_to_animation(easing)

Converts an easing to an Animation with time goes from `0.0` to `1.0` and value from `0` to `1`.
"""
easing_to_animation(easing) = Animation(0.0, 0.0, easing, 1.0, 1.0)

"""
    easing_to_animation(rev_easing::ReversedEasing)

Converts an easing to an Animation with time goes from `0.0` to `1.0` and value from `1` to `0`.
"""
easing_to_animation(rev_easing::ReversedEasing) =
    Animation(0.0, 1.0, rev_easing.easing, 1.0, 0.0)


text:
@testset "Sanity Matrix/Array" begin
    @test distance_matrix([ ], [ ]) == 0

    @testset "Sanity Matrix/Array/Left" begin
        @test distance_matrix([ 1 ], [ ])     == 1
        @test distance_matrix([ 1 2 ], [ ])   == 2
        @test distance_matrix([ 1 2 3 ], [ ]) == 3
    end

    @testset "Sanity Matrix/Array/Right" begin
        @test distance_matrix([ ], [ 1 ])     == 1
        @test distance_matrix([ ], [ 1 2 ])   == 2
        @test distance_matrix([ ], [ 1 2 3 ]) == 3
    end

    @testset "Sanity Matrix/Array/Both" begin
        @test distance_matrix([ 1 ], [ 1 ])         == 0
        @test distance_matrix([ 1 ], [ 2 ])         == 1
        @test distance_matrix([ 1 2 ], [ 2 2 ])     == 1
        @test distance_matrix([ 1 2 3 ], [ 2 2 3 ]) == 1
        @test distance_matrix([ 1 2 3 ], [ 1 2 2 ]) == 1
    end
end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
include("test_functional.jl")
include("test_issue17.jl")

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
import Base.Random: rand, zeros

export MetaPopulation, popsize, popct, popfits

@doc """
Meta populations can be used to simulation spatially structured populations,
such as those that are geographically isolated from one another.

A meta population can be used in all instances where an ordinary population is
valid.
"""
type MetaPopulation
  populations::Vector{Population}
end

MetaPopulation(mp::MetaPopulation) = MetaPopulation([Population(p) for p = mp.populations])

@doc """rand(::Type{MetaPopulation}, ls::Landscape, popsize::Int64, popct::Int64)

Create a random meta population based on the given landscape. Each population
will consist of `popsize` individuals and there will be `popct` individual
populations.
"""
rand(::Type{MetaPopulation}, ls::Landscape, popsize::Int64, popct::Int64) =
    MetaPopulation([rand(Population, ls, popsize) for _ = 1:popct])

@doc """zeros(::Type{MetaPopulation}, ls::Landscape, popsize::Int64, popct::Int64)

Create a meta population in which all loci "zeroed". Each population will
consist of `popsize` individuals and there will be `popct` individual
populations.
"""
zeros(::Type{MetaPopulation}, ls::Landscape, popsize::Int64, popct::Int64) =
    MetaPopulation([zeros(Population, ls, popsize) for _ = 1:popct])

@doc """popsize(p::MetaPopulation)

Return the number of individuals in each population within the meta population.
"""
popsize(p::MetaPopulation) = length(p.populations[1].genotypes)

@doc """popct(p::MetaPopulation)

Return the number of populations within the meta population.
"""
popct(p::MetaPopulation) = length(p.populations)

@doc """popfits(p::MetaPopulation, ls::Landscape)

Return a matrix of fitness values where the ith row represents the fitness
values of the ith member of each population. Each column represents a
population.
"""
function popfits(p::MetaPopulation)
  fs = zeros(Float64, popsize(p), popct(p))
  for ip = 1:popct(p)
    for ig = 1:popsize(p)
      fs[ig, ip] = fitness(p.populations[ip].genotypes[ig])
    end
  end
  return fs
end


text:
function null_error_handler(error_code::Cint, mod::Char, func::Char, eh_data::Ptr{Cvoid})
    nothing
end


text:
"""
    vandermonde(elem::AbstractElemShape, N, rst...)

Computes the generalized Vandermonde matrix V of degree N at points (r,s,t).
"""
vandermonde(elem::AbstractElemShape, N, rst...) = first(basis(elem,N,rst...))
vandermonde(elem::Line, N, r) = first(basis(elem,N,r)) # specialize for 1D

"""
    grad_vandermonde(elem::AbstractElemShape, N, rst...)

Computes the generalized Vandermonde derivative matrix V of degree N at points (r,s,t).
"""
grad_vandermonde(elem::AbstractElemShape, N, rst...) = basis(elem, N, rst...)[2:end]
grad_vandermonde(elem::Line, N, r) = last(basis(elem,N,r)) # specialize for 1D

"""
    nodes(elem::AbstractElemShape,N)

Computes interpolation nodes of degree N. Edge nodes coincide with (N+1)-point Lobatto points.
Default routine for elem = Tet(), Pyr(), Tri().

For Quad(), Hex(), Wedge() elements, nodes(...) returns interpolation points constructed
using a tensor product of lower-dimensional nodes.
"""
nodes(elem::AbstractElemShape,N) = build_warped_nodes(elem,N,nodes(Line(),N))

"""
    basis(elem::AbstractElemShape, N, rst...)

Computes orthonormal basis of degree N at tuple of coordinate arrays (r,s,t).
"""
basis(elem::AbstractElemShape,N,rst...)

"""
    equi_nodes(elem::AbstractElemShape, N)

Compute equispaced nodes of degree N.
"""
equi_nodes(elem::AbstractElemShape, N)

"""
    quad_nodes(elem::AbstractElemShape, N)

Compute quadrature nodes and weights exact for (at least) degree 2N polynomials.
"""
quad_nodes(elem::AbstractElemShape, N)

"""
    stroud_quad_nodes(elem::AbstractElemShape,N)

Returns Stroud-type quadrature nodes and weights constructed from the tensor product
of (N+1)-point Gauss-Jacobi rules. Exact for degree 2N polynomials
"""
stroud_quad_nodes(elem::AbstractElemShape,N)


text:
##############################################################################################
module TestClassicSmootherExps
##############################################################################################
##############################################################################################
# imports and exports
using DataAssimilationBenchmarks.SmootherExps
using JLD2, Statistics

##############################################################################################
##############################################################################################
# run and analyze the ETKS for state estimation with the Lorenz-96 model
# arguments are
# time_series, method, seed, nanl, lag, shift, obs_un, obs_dim, γ, N_ens, infl = args

function run_smoother_state_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "etks", 0, 3500, 10, 1, 1.0, 40, 1.00, 21, 1.02)
        classic_state(args)
        true
    catch
        false
    end
end

function analyze_smoother_state_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/etks-classic/")
        data = load(path * "etks-classic_L96_state_seed_0000_diff_0.000_sysD_40_obsD_40_" *
                    "obsU_1.00_gamma_001.0_nanl_03500_tanl_0.05_h_0.05_lag_010_shift_001_" *
                    "mda_false_nens_021_stateInfl_1.02.jld2")
        filt_rmse = data["filt_rmse"]
        post_rmse = data["post_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if mean(filt_rmse[501:end]) < 0.2 && mean(post_rmse[501:end]) < 0.15
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# run and analyze the ETKF for joint state-parameter estimation with the Lorenz-96 model
# time_series, method, seed, nanl, lag, shift, obs_un, obs_dim, γ,
# param_err, param_wlk, N_ens, state_infl, param_infl = args

function run_smoother_param_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "etks", 0, 3500, 10, 1, 1.0, 40, 1.0,
                0.10, 0.0010, 21, 1.02, 1.0)
        classic_param(args)
        true
    catch
        false
    end
end

function analyze_smoother_param_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/etks-classic/")
        data = load(path * "etks-classic_L96_param_seed_0000_diff_0.000_sysD_41_obsD_40_" *
                    "obsU_1.00_gamma_001.0_paramE_0.10_paramW_0.0010_nanl_03500_tanl_0.05_" *
                    "h_0.05_lag_010_shift_001_mda_false_nens_021_stateInfl_1.02_" *
                    "paramInfl_1.00.jld2")
        filt_rmse = data["filt_rmse"]
        post_rmse = data["post_rmse"]
        para_rmse = data["param_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if (mean(filt_rmse[501:end]) < 0.2) && (mean(post_rmse[501:end]) < 0.15) &&
            (mean(para_rmse[501:end]) < 0.01)
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# end module

end

max_stars_repo_name: cgrudz/DataAssimilationBenchmarks

max_issues_repo_name: cgrudz/DataAssimilationBenchmarks


text:
export add_gate, gate, gate_help

# These are tables of the gates and their descriptions
gate_table = Dict{Symbol, Matrix{CC}}()
gate_descr = Dict{Symbol,String}()

"""
`add_gate(S::Symbol, M::Matrix, descr::String)` adds a new gate to the
library of gates.
* `S` is the symbol associated with this gate.
* `M` is the matrix for this gate.
* `descr` is a description of the action of the gate.
"""
function add_gate(S::Symbol, M::Matrix, descr::String="No description available")
    gate_table[S] = M
    gate_descr[S] = descr
end

"""
`gate(S::Symbol)` returns the gate associated with the symbol `S`.
"""
function gate(S::Symbol)
    if haskey(gate_table,S)
        return gate_table[S]
    end
    throw(error("No such gate, $S"))
end

"""
`gate_help(S::Symbol)` prints a short help message describing gate `S`.
"""
function gate_help(S::Symbol)
    if haskey(gate_descr,S)
        println(":$S\t",gate_descr[S])
        return
    end
    println("Unknown gate, $S")
end

add_gate(:I, [1 0; 0 1], "Identity (nop) gate")
add_gate(:X, [0 1; 1 0], "Qubit swap (Pauli X)")
add_gate(:H, [1 1; 1 -1]/sqrt(2), "Hadammard gate")


text:
# workaround #60:
@osx_only ENV["PATH"] = JULIA_HOME*":"*ENV["PATH"]

include("IJulia.jl")
include("inline.jl")

using IPythonDisplay
pushdisplay(InlineDisplay())

using IJulia

ccall(:jl_install_sigint_handler, Void, ())

# the size of truncated output to show should not depend on the terminal
# where the kernel is launched, since the display is elsewhere
ENV["LINES"] = 30
ENV["COLUMNS"] = 80

println(IJulia.orig_STDOUT, "Starting kernel event loops.")
IJulia.watch_stdio()

# workaround JuliaLang/julia#4259
delete!(task_local_storage(),:SOURCE_PATH)

IJulia.waitloop()


text:
using MbedTLS, Test, Random, Distributed

# by adding a single worker, we load the julia-shipped mbedtls binaries and thus ensure
# the rest of the MbedTLS.jl tests run propertly even in the presence of alternative
# versioned mbedtls libraries
addprocs(1)

import Sockets: connect

# Message digests
@test digest(MD_SHA256, "test", "secret") ==
   UInt8[0x03, 0x29, 0xa0, 0x6b, 0x62, 0xcd, 0x16, 0xb3, 0x3e, 0xb6, 0x79, 0x2b, 0xe8, 0xc6, 0x0b, 0x15, 0x8d, 0x89, 0xa2, 0xee, 0x3a, 0x87, 0x6f, 0xce, 0x9a, 0x88, 0x1e, 0xbb, 0x48, 0x8c, 0x09, 0x14]

# Symmetric encryption
let
    secret_key = rand(UInt8, 32)
    message = "Testing symmetric encryption"
    iv  = rand(UInt8, 16)

    cipher_text = encrypt(CIPHER_AES, secret_key, message, iv)
    plain_text = decrypt(CIPHER_AES, secret_key, cipher_text, iv)

    @test message == String(plain_text)
end

# openssl rsa -inform PEM -text -noout < key.pem
# keyPEMN =
# """
#     00:bc:42:e0:e8:31:52:fd:af:e3:32:42:c4:28:5b:
#     e9:cf:50:af:e9:20:c7:96:37:5c:cb:09:e6:90:5c:
#     ac:fc:8c:03:22:7c:68:70:32:f0:7b:e4:52:85:a2:
#     12:61:ac:82:16:a6:66:e8:e9:a2:f4:dc:b7:05:0f:
#     e0:21:ea:0d:6d:0d:3e:fe:5e:c2:6a:59:36:52:c9:
#     66:5a:54:b9:f7:80:43:41:48:56:b5:b3:50:bd:2a:
#     29:8c:81:96:29:64:0a:80:82:a4:dc:93:28:8f:43:
#     d5:0b:4b:3a:4f:c4:be:a0:ff:3f:c6:68:16:f0:c5:
#     be:04:16:fc:b6:52:c4:d8:f9:4d:66:bc:d6:b8:ce:
#     4f:a4:5a:97:66:fb:5d:db:1c:e9:e6:89:ab:f6:82:
#     9e:de:93:f7:3d:b4:35:77:5e:ae:1f:67:14:29:a4:
#     00:df:7c:2c:3d:76:42:d5:76:08:ff:11:09:20:fc:
#     bd:8a:8d:d4:a9:ce:2e:35:2b:c6:d2:de:dc:ad:1b:
#     8d:01:7c:c5:32:9a:8e:c4:f7:a6:94:55:d3:4b:96:
#     1b:ee:a0:94:95:5f:b2:a0:b0:f8:bb:02:b0:a5:a9:
#     0e:62:f1:a2:8a:3a:6f:ec:c2:e9:5e:b5:73:45:cb:
#     35:86:24:07:e8:11:28:25:7d:ce:35:c0:48:ff:a0:
#     f6:13
# """
keyPEMN = 23765780372307576820243785798710539106220019827391701276528023465129289386491638200309724085345164197157471148949804917921135961746756129955755123246417065664000080083079385610461361129162124204308601470870363093163972105641380345473969785909726960802402051498557030325897388318234260665809255594699414895859065471822450027547355233222204156733689900551990147367094140650142549959621935662238290254703866686898285439063843022440258882238926771918029490099100463478790421686785230680564068945032213969236010272566088532861337469543246408184082656216674316765444657485587779489134661155841677954252447171077733538133523
keyPEMe = BigInt(65537)

function verify_key_pem(data, signature)
    pubkey = RSA(MbedTLS.MBEDTLS_RSA_PKCS_V15, MD_SHA1)
    @test_throws MbedTLS.MbedException MbedTLS.pubkey_from_vals!(pubkey, keyPEMN, keyPEMe)
    MbedTLS.pubkey_from_vals!(pubkey, keyPEMe, keyPEMN)
    MbedTLS.verify(pubkey, MD_SHA1, MbedTLS.digest(MD_SHA1, data), signature)
end

# RSA
let
    key = MbedTLS.gen_key(MersenneTwister(0))
    # todo: test encryption/decryption

    # Test signature verification.
    # Pre-generated signature of "MbedTLS.jl" using key.pem
    signature =
    UInt8[0x35,0x6d,0xa0,0xdd,0x02,0x06,0xc2,0x8e,0x21,0xcf,0x34,0x48,0x59,0x39,
          0x61,0x2c,0x37,0x50,0x74,0x2f,0x1c,0x25,0x49,0x52,0x1c,0xaf,0xb0,0xb5,
          0x2d,0x52,0x4d,0xa0,0x45,0x86,0x23,0xc8,0xf4,0xc6,0xb1,0xa2,0x25,0x86,
          0x5b,0xe6,0xcc,0x34,0x21,0xd6,0x31,0x36,0x0c,0xaf,0x97,0xa4,0x7e,0xf4,
          0xe3,0x1a,0x5f,0x58,0x78,0x57,0xbf,0x1f,0xfd,0x5e,0xe0,0xc2,0x9a,0x33,
          0x66,0xb2,0x4b,0x41,0xb2,0x37,0xac,0x63,0xd2,0x7d,0x9b,0x76,0x2f,0xe8,
          0x66,0x1e,0x83,0x4f,0x45,0x47,0x76,0x3f,0xcc,0x91,0xbb,0x8b,0xc2,0x12,
          0x8a,0x69,0x83,0x39,0xca,0x87,0xbd,0x36,0xa0,0x5e,0x21,0x60,0x39,0x9b,
          0xf0,0x3c,0x84,0x19,0x92,0xfc,0x76,0xaa,0x71,0xd8,0x39,0x0c,0x5a,0x2c,
          0x82,0xe3,0x77,0xc1,0xcc,0x7b,0xb0,0xfb,0x3f,0xb4,0x4a,0xf1,0xaa,0xa8,
          0xae,0x1e,0x46,0xde,0x42,0xc4,0x1d,0x9b,0xe9,0x58,0xaf,0x36,0xb5,0xbf,
          0x35,0x92,0x94,0x5e,0x73,0xca,0x83,0x34,0x05,0xe3,0xf5,0x8c,0x18,0x3a,
          0xb6,0xea,0xe0,0xc9,0x54,0x83,0x13,0x1f,0xd0,0x4e,0x4e,0x26,0xdb,0xec,
          0x8a,0x98,0x99,0x03,0x09,0x49,0xd5,0x10,0x2f,0xc4,0x61,0xc7,0x9d,0x5f,
          0x37,0x10,0x49,0x28,0x72,0xb5,0x43,0xba,0x96,0x98,0x2c,0xb1,0x86,0x6e,
          0x04,0x80,0x61,0xdb,0x00,0x77,0x31,0x14,0xf3,0x31,0x0e,0xc6,0xda,0x40,
          0x73,0xe7,0x3e,0x1f,0xa6,0x45,0x29,0x43,0xcb,0xe2,0xff,0x1e,0xe8,0xaa,
          0x4f,0x0b,0x01,0xcb,0x4f,0xf2,0xe4,0xc9,0xfd,0x1d,0xb8,0x5c,0x40,0xf2,
          0xeb,0x76,0x2b,0x5f]
    verify_key_pem("MbedTLS.jl", signature)
end

# Basic TLS client functionality
let
    testhost = "httpbin.org"
    sock = connect(testhost, 443)
    entropy = MbedTLS.Entropy()

    function entropy_func(buf)
        buf[:] = rand(RandomDevice(), UInt8, length(buf))
        return Cint(length(buf))
    end

    MbedTLS.add_source!(entropy, entropy_func, 0, true)
    rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)

    ctx = MbedTLS.SSLContext()
    conf = MbedTLS.SSLConfig()

    MbedTLS.config_defaults!(conf)
    MbedTLS.authmode!(conf, MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED)
    MbedTLS.rng!(conf, rng)

    function show_debug(level, filename, number, msg)
        @show level, filename, number, msg
    end

    MbedTLS.dbg!(conf, show_debug)
    MbedTLS.set_dbg_level(MbedTLS.DebugThreshold(1))

    MbedTLS.ca_chain!(conf)

    MbedTLS.setup!(ctx, conf)
    MbedTLS.set_bio!(ctx, sock)
    MbedTLS.hostname!(ctx, testhost)
    MbedTLS.handshake(ctx)

    write(ctx, "GET / HTTP/1.1\r\nHost: $testhost\r\n\r\n")
    buf = String(read(ctx, 100))
    @test occursin(r"^HTTP/1.1 200 OK", buf)
end

# Test ALPN
let
    testhost = "google.com"
    sock = connect(testhost, 443)
    entropy = MbedTLS.Entropy()

    function entropy_func(buf)
        buf[:] = rand(RandomDevice(), UInt8, length(buf))
        return Cint(length(buf))
    end

    MbedTLS.add_source!(entropy, entropy_func, 0, true)
    rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)

    ctx = MbedTLS.SSLContext()
    conf = MbedTLS.SSLConfig()

    MbedTLS.config_defaults!(conf)
    MbedTLS.set_alpn!(conf, ["h2"])
    MbedTLS.authmode!(conf, MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED)
    MbedTLS.rng!(conf, rng)

    function show_debug(level, filename, number, msg)
        @show level, filename, number, msg
    end

    MbedTLS.dbg!(conf, show_debug)
    MbedTLS.set_dbg_level(MbedTLS.DebugThreshold(1))

    MbedTLS.ca_chain!(conf)

    MbedTLS.setup!(ctx, conf)
    MbedTLS.set_bio!(ctx, sock)
    MbedTLS.hostname!(ctx, testhost)
    MbedTLS.handshake(ctx)

    write(ctx, "PRI * HTTP/2.0\r\n\r\nSM\r\n\r\n")
    @test MbedTLS.alpn_proto(ctx) == "h2"
end

# Test pk.jl methods
let
    key = MbedTLS.parse_keyfile(joinpath(@__DIR__, "key.pem"))
    @test MbedTLS.bitlength(key) == 2048
    @test MbedTLS.get_name(key) == "RSA"
    output = fill(0x00, 256)
    @test sizeof(output) == MbedTLS.sign!(key, MD_SHA1,
        MbedTLS.digest(MD_SHA1, "MbedTLS.jl"), output, MersenneTwister(0))
    verify_key_pem("MbedTLS.jl", output)

    signature = MbedTLS.sign(key, MD_SHA1,
        MbedTLS.digest(MD_SHA1, "MbedTLS.jl"), MersenneTwister(0))
    @test signature == output[1:length(signature)]

    key_string = read(open(joinpath(@__DIR__, "key.pem"), "r"))
    key = MbedTLS.PKContext()
    MbedTLS.parse_key!(key, key_string)
    @test MbedTLS.bitlength(key) == 2048
    @test MbedTLS.get_name(key) == "RSA"

    pubkey_string = read(open(joinpath(@__DIR__, "public_key.pem"), "r"))
    pubkey = MbedTLS.PKContext()
    MbedTLS.parse_public_key!(pubkey, pubkey_string)
    @test MbedTLS.bitlength(pubkey) == 2048
    @test MbedTLS.get_name(pubkey) == "RSA"

    key = MbedTLS.parse_keyfile(joinpath(@__DIR__, "key.pem"))
    @test MbedTLS.bitlength(key) == 2048
    @test MbedTLS.get_name(key) == "RSA"

    pubkey = MbedTLS.parse_public_keyfile(joinpath(@__DIR__, "public_key.pem"))
    @test MbedTLS.bitlength(pubkey) == 2048
    @test MbedTLS.get_name(pubkey) == "RSA"

    for md in instances(MbedTLS.MDKind)
        if in(md, (MD_NONE, MD_MD2, MD_MD4))
            continue
        end
        hash = MbedTLS.digest(md, "MbedTLS.jl")
        signature = MbedTLS.sign(key, md, hash, MersenneTwister(0))
        @test MbedTLS.verify(pubkey, md, hash, signature) == 0
    end
end

# Test md.jl
let
    md = MbedTLS.MD(MD_SHA1)
    write(md, UInt8['M', 'b', 'e', 'd'])
    a32 = fill(UInt32(0), 1)
    reinterpret(UInt8, a32) .= UInt8['T', 'L', 'S', '.']
    write(md, a32[])
    a16 = fill(UInt16(0), 1)
    reinterpret(UInt8, a16) .= UInt8['j', 'l']
    write(md, a16)
    @test MbedTLS.finish!(md) == MbedTLS.digest(MD_SHA1,b"MbedTLS.jl")

    # Test reset functionality
    md = MbedTLS.MD(MbedTLS.MD_SHA256, "passcode")
    write(md, "msg")
    digest1 = MbedTLS.finish!(md)
    MbedTLS.reset!(md)
    write(md, "msg")
    digest2 = MbedTLS.finish!(md)
    @test digest1 == digest2
    write(md, "msg")
    digest3 = MbedTLS.finish!(md)
    @test digest1 ≠ digest3
end

# log_secrets
mktempdir() do d

    f = joinpath(d, "secrets.log")

    testhost = "httpbin.org"
    sock = connect(testhost, 443)

    ctx = MbedTLS.SSLContext()
    conf = MbedTLS.SSLConfig(true; log_secrets=f)

    MbedTLS.setup!(ctx, conf)
    MbedTLS.set_bio!(ctx, sock)
    MbedTLS.hostname!(ctx, testhost)
    MbedTLS.handshake(ctx)

    write(ctx, "GET / HTTP/1.1\r\nHost: $testhost\r\n\r\n")
    buf = String(read(ctx, 100))
    @test occursin(r"^HTTP/1.1 200 OK", buf)
    @test occursin(r"^CLIENT_RANDOM [0-9a-f]{64} [0-9a-f]{96}$", read(f, String))
end

let
    include("clntsrvr/clntsrvr.jl")
end


text:
using MultivectorAlgebras
using Base.Test: @testset

@testset "1-dimensional" begin include("one-dimensional_tests.jl") end
@testset "2-dimensional" begin include("two-dimensional_tests.jl") end


text:
using SchoolBusRouting, JLD, DataFrames, CSV
include("SchoolBusRouting.jl")

# experiment parameters
experiments = vcat(["RSRB0$i" for i = 1:8], ["CSCB0$i" for i = 1:8]);
maxtimes = [2700.0, 5400.0];
conditions = ["one", "many", "Chen", "combined", "LBH"]
maxiter = 1

# output arrays
exp_id = Int[]; benchmark = AbstractString[]; exp_name = AbstractString[];
max_riding_time = Float64[]; nbuses = Int[]; iter = Int[];

counter = 0

experiment = "RSRB01"
maxtime = 2700.0
condition = "many"
for maxtime in (maxtimes), experiment in (experiments)
	# load data from ASCII files
	PATH = "$(Pkg.dir("SchoolBusRouting"))/data/input"
	data = SchoolBusRouting.loadSyntheticBenchmark("$PATH/$experiment/Schools.txt",
									  "$PATH/$experiment/Stops.txt")
	data = SchoolBusRouting.loadSyntheticBenchmark("data/input/CSCB01/Schools.txt",
						 "data/input/CSCB01/Stops.txt")
	data.params.max_time_on_bus = maxtime
	println("***********************************************")
	println("$experiment - $maxtime")
	for condition in conditions, i=1:maxiter
		counter += 1
		scenariolist = []
		params = []
		if condition in ["many", "combined"]
			# many scenarios
			# λs = [1e2, 5e2, 5e3, 1e4, 5e4, 1e5, 5e5, 1e6];
			λs = [1e2, 5e2, 1e3, 2e3, 5e3];
			nGreedy = 20; nIterations = 80;
			if experiment == "RSRB02" && maxtime < 4000
				# λs = [1e2, 5e2, 1e3, 2e3, 5e3, 7e3, 1e4, 5e4, 1e5, 5e5, 1e6, 5e6, 1e7, 2e7];
				nGreedy = 40; nIterations = 150;
			end
			params = [SchoolBusRouting.ScenarioParameters(maxRouteTimeLower=maxtime-2000,
											 maxRouteTimeUpper=maxtime+2000,
											 nGreedy=nGreedy, λ=λ,
											 nIterations=nIterations) for λ in λs]
            params = [ScenarioParameters(maxRouteTimeLower=maxtime-2000,
											 maxRouteTimeUpper=maxtime+2000,
											 nGreedy=nGreedy, λ=λ,
											 nIterations=nIterations) for λ in λs]
			# srand(i)
			Random.seed!(i)
			scenariolist = computescenarios(data, params;
			                               ## Gurobi arguments
			 									# OutputFlag=0,
			                                    # MIPGap=ifelse(maxtime < 4000, 1e-4, 0.05),
			                                    # TimeLimit=ifelse(maxtime < 4000, 90, 30)
                                           ## Cbc arguments
												logLevel=1,
			                                    allowableGap=ifelse(maxtime < 4000, 1e-4, 0.05),
			                                    seconds=ifelse(maxtime < 4000, 90, 30),
												threads = 8
												);
		end
		if condition == "one"
			params = [SchoolBusRouting.ScenarioParameters(maxRouteTimeLower=maxtime-2000,
											 maxRouteTimeUpper=maxtime+2000,
											 nGreedy=ifelse(maxtime < 4000, 30, 20), λ=λ,
											 nIterations=ifelse(maxtime < 4000, 120, 100)
											 ) for λ in [1e8]]
			srand(i)
			scenariolist = SchoolBusRouting.computescenarios(data, params;
			                               ## Gurobi arguments
			 									# OutputFlag=0,
			                                    # MIPGap=ifelse(maxtime < 4000, 1e-4, 0.05),
			                                    # TimeLimit=ifelse(maxtime < 4000, 90, 30)
                                           ## Cbc arguments
												logLevel=1,
			                                    allowableGap=ifelse(maxtime < 4000, 1e-4, 0.05),
			                                    seconds=ifelse(maxtime < 4000, 90, 30),
												threads = 8
												);
		end
		if condition in ["Chen", "combined"]
			mapping = Dict(stop.originalId => (stop.schoolId,
			                                   stop.id) for stop in vcat(data.stops...))
			routestops = [Vector{Int}[] for school in data.schools]
			servicetimes = [Float64[] for school in data.schools]
			found = false
			for line in readlines("$(Pkg.dir("SchoolBusRouting"))/data/routes/$(experiment[1:4])$(Int(round(maxtime)))/SBRP_$experiment.txt")
				if ((line != "NODE_COORD_SECTION" && !found) ||
					(split(line, ":")[1] == "DIMENSION" && found))
					found = false
					continue
				elseif !found
					found = true
					continue
				end
				elements = split(line, "\t")
				stoplist = [mapping[parse(Int64, a)][2] for a in elements[11:end]]
				servicetime = parse(Float64, elements[4])
				destination = parse(Int64, elements[8]) - 200000
				@assert all([mapping[parse(Int64,a)][1] for a in elements[11:end]] .== destination)
				if length(stoplist) > 0
					push!(routestops[destination], stoplist)
					push!(servicetimes[destination], servicetime)
				end
			end
			routes = [[SchoolBusRouting.Route(j, stops) for (j, stops) in enumerate(routestops[i])] for i = eachindex(routestops)];
			scenariolist = vcat(scenariolist,
			                    [(SchoolBusRouting.Scenario(i, length(params)+1, collect(eachindex(routes[i]))),
			                     routes[i]) for i=eachindex(routes)])
		end
		if condition == "LBH"
			scenarios, routes, usedscenario, buses = SchoolBusRouting.lbh(data)
			SchoolBusRouting.loadsolution!(data, scenarios, routes, usedscenario, buses)
		else
			SchoolBusRouting.loadroutingscenarios!(data, scenariolist);
			SchoolBusRouting.routeBuses!(data, OutputFlag=0, TimeLimit=3600);
		end
		SchoolBusRouting.test(data)
		println("$counter - $condition: $(length(data.buses)) buses")
		CSV.write("$(Pkg.dir("SchoolBusRouting"))/output/tableS1/solution-$counter.csv", SchoolBusRouting.output(data))
		push!(nbuses, length(data.buses))
		push!(exp_id, counter)
		push!(benchmark, experiment)
		push!(max_riding_time, maxtime)
		push!(exp_name, condition)
		push!(iter, i)
		JLD.save("$(Pkg.dir("SchoolBusRouting"))/output/tableS1/solution-$counter.jld", "data", data)
	end
end
df = DataFrame()
df[:exp_id] = exp_id; df[:benchmark] = benchmark; df[:exp_name] = exp_name;
df[:max_riding_time] = max_riding_time; df[:buses] = nbuses; df[:iter] = iter;
CSV.write("$(Pkg.dir("SchoolBusRouting"))/output/pnas-tableS1.csv", df)


text:
# Test of the A* algorithm for finding the shortest path between two vertices

using Graphs
using Base.Test

# The "touring Romania" example from Russell and Norvig
g1 = simple_inclist(20, is_directed=false)

g1_wedges = [
    (1, 20, 75.),
    (1, 16, 140.), #2, Arad -- Sibiu
    (1, 17, 118.),
    (20, 13, 71.),
    (13, 16, 151.),
    (16, 6, 99.),
    (16, 15, 80.), #7, Sibiu -- Rimnicu Vilcea
    (6, 2, 211.),
    (17, 10, 111.),
    (10, 11, 70.),
    (11, 4, 75.),
    (4, 3, 120.),
    (3, 15, 146.),
    (3, 14, 138.),
    (15, 14, 97.), #15, Rimnicu Vilcea -- Pitesti
    (14, 2, 101.), #16, Pitesti -- Bucharest
    (2, 7, 90.),
    (2, 18, 85.),
    (18, 8, 98.),
    (8, 5, 86.),
    (18, 19, 142.),
    (19, 9, 92.),
    (9, 12, 87.) ]

g1_heuristics = [
    366, 0, 160, 242, 161, 176, 77, 151, 226, 244, 241, 234, 380, 100, 193,
    253, 329, 80, 199, 374 ]

ne = length(g1_wedges)
eweights1 = zeros(ne)
for i = 1 : ne
    we = g1_wedges[i]
    add_edge!(g1, we[1], we[2])
    eweights1[i] = we[3]
end

sp = shortest_path(g1, eweights1, 1, 2, n -> g1_heuristics[n])
edge_numbers = map(e -> edge_index(e, g1), sp)
@test edge_numbers == [2, 7, 15, 16]


text:
"""
This function does nothing if project is not in queue.
Returns nothing.
"""
function start_construction!(projects::Vector{<: Project{<: BuildPhase}},
                             index::Int64,
                             project::P,
                             iteration_year::Int64) where P <: Project{<: BuildPhase}
    return
end

"""
This function converts a project from queue to planned if
queue time has been completed.
Returns nothing.
"""
function start_construction!(projects::Vector{<: Project{<: BuildPhase}},
                             index::Int64,
                             project::P,
                             iteration_year::Int64) where P <: Project{Queue}

    queue_time = length(get_queue_cost(get_finance_data(project)))

     if get_decision_year(project) + queue_time == iteration_year
        projects[index] = convert(Project{Planned}, project)
     end
end

"""
This function does nothing if project is not a planned project.
Returns nothing.
"""
function finish_construction!(projects::Vector{<: Project{<: BuildPhase}},
                             index::Int64,
                             project::P,
                             sys_UC::Union{Nothing, PSY.System},
                             services::Vector{PSY.Service},
                             simulation_dir::String,
                             iteration_year::Int64,
                             start_year::Int64) where P <: Project{<: BuildPhase}
    return
end

"""
This function converts a project from planned to existing if
construction time has been completed.
Returns nothing.
"""
function finish_construction!(projects::Vector{<: Project{<: BuildPhase}},
                             index::Int64,
                             project::P,
                             sys_UC::Nothing,
                             services::Vector{PSY.Service},
                             simulation_dir::String,
                             iteration_year::Int64,
                             start_year::Int64) where P <: Project{Planned}

     if get_construction_year(project) == iteration_year
        projects[index] = convert(Project{Existing}, project)
     end

     type = get_type(get_tech(project))
     zone = get_zone(get_tech(project))

     availability_df = read_data(joinpath(simulation_dir, "timeseries_data_files", "Availability", "DAY_AHEAD_availability.csv"))
     if in(get_name(project), names(availability_df))
         availability_raw = availability_df[:, Symbol(get_name(project))]
     elseif in("$(type)_$(zone)", names(availability_df))
         availability_raw = availability_df[:, Symbol("$(type)_$(zone)")]
     end

    load_n_vg_df =  read_data(joinpath(simulation_dir, "timeseries_data_files", "Net Load Data", "load_n_vg_data.csv"))
    load_n_vg_df[:, get_name(project)] = availability_raw * get_maxcap(project)

    write_data(joinpath(simulation_dir, "timeseries_data_files", "Net Load Data"), "load_n_vg_data.csv", load_n_vg_df)

    return

end

"""
This function converts a project from planned to existing if
construction time has been completed.
Returns nothing.
"""
function finish_construction!(projects::Vector{<: Project{<: BuildPhase}},
                             index::Int64,
                             project::P,
                             sys_UC:: PSY.System,
                             services::Vector{PSY.Service},
                             simulation_dir::String,
                             iteration_year::Int64,
                             start_year::Int64) where P <: Project{Planned}

     if get_construction_year(project) == iteration_year
        projects[index] = convert(Project{Existing}, project)
        PSY_project = create_PSY_generator(project, sys_UC)

        PSY.add_component!(sys_UC, PSY_project)

        for product in get_products(project)
            add_device_services!(services, PSY_project, product)
        end

        type = get_type(get_tech(project))
        zone = get_zone(get_tech(project))

        availability_df = read_data(joinpath(simulation_dir, "timeseries_data_files", "Availability", "DAY_AHEAD_availability.csv"))
        if in(get_name(project), names(availability_df))
            availability_raw = availability_df[:, Symbol(get_name(project))]
        elseif in("$(type)_$(zone)", names(availability_df))
            availability_raw = availability_df[:, Symbol("$(type)_$(zone)")]
        end

        add_device_forecast!(simulation_dir, sys_UC, PSY_project, availability_raw, start_year)

     end

     return
end

"""
This function does nothing if project is not an existing project.
Returns nothing.
"""
function retire_old!(projects::Vector{<: Project{<: BuildPhase}},
                     index::Int64,
                     project::P,
                     sys::Union{Nothing, PSY. System},
                     simulation_dir::String,
                     iteration_year::Int64) where P <: Project{<: BuildPhase}
    return false
end

"""
This function converts a project from existing to retired if
lifetime has ended.
Returns nothing.
"""
function retire_old!(projects::Vector{<: Project{<: BuildPhase}},
                     index::Int64,
                     project::P,
                     sys::Union{Nothing, PSY. System},
                     simulation_dir::String,
                     iteration_year::Int64) where P <: Project{Existing}
    if get_end_life_year(project) == iteration_year
        set_retirement_year!(projects[index], iteration_year + 1)
        projects[index] = convert(Project{Retired}, project)
        remove_system_component!(sys, project)
        remove_renewable_gen_data!(project, simulation_dir)
    end
end


text:
using TopOpt.MMA
using Test, LinearAlgebra

using TopOpt.MMA: eval_constraint, eval_objective, dim

function f(x::AbstractVector, grad::AbstractVector)
    if length(grad) != 0
        grad[1] = 0.0
        grad[2] = 0.5/sqrt(x[2])
    end
    sqrt(x[2])
end

function g(x::AbstractVector, grad::AbstractVector, a, b)
    if length(grad) != 0
        grad[1] = 3a * (a*x[1] + b)^2
        grad[2] = -1
    end
    (a*x[1] + b)^3 - x[2]
end

ndim = 2
m = MMA.Model(ndim, f)

box!(m, 1, 0.0, 10.0)
box!(m, 2, 0.0, 10.0)

ineq_constraint!(m, (x,grad) -> g(x,grad,2,0))
ineq_constraint!(m, (x,grad) -> g(x,grad,-1,1))

@test dim(m) == 2

# Objective
let
    grad1 = zeros(2)
    grad2 = zeros(2)
    p = [1.234, 2.345]
    @test eval_objective(m, p, grad2) ≈ f(p, grad1)
    @test grad1 ≈ grad2
end

# Box
@test min(m, 1) == 0.0
@test max(m, 1) == 10.0
@test min(m, 2) == 0.0
@test max(m, 2) == 10.0

# Inequalities
let
    grad1 = zeros(2)
    grad2 = zeros(2)
    p = [1.234, 2.345]
    @test eval_constraint(m, 1, p, grad1) ≈ g(p,grad2 ,2,0)
    @test grad1 ≈ grad2

    @test eval_constraint(m, 2, p, grad1) ≈ g(p,grad2,-1,1)
    @test grad1 ≈ grad2
end

r = optimize(m, [1.234, 2.345], MMA.MMA02(), options = MMA.Options(tol = MMA.Tolerances(xtol = 1e-6, ftol = 0.0, grtol = 0.0)))
@test abs(r.minimum - sqrt(8/27)) < 1e-6
@test norm(r.minimizer - [1/3, 8/27]) < 1e-6

r = optimize(m, [1.234, 2.345], MMA.MMA87())
@test abs(r.minimum - sqrt(8/27)) < 1e-6
@test norm(r.minimizer - [1/3, 8/27]) < 1e-6


text:
function Ψ_search(mdl::Model)

    tm = mdl.DR
    mdl.Ψ = zeros( Float64, size(mdl.VF) )
    # Probabilities = ones(t_V,2)
    mdl.Λ = vec( (mdl.λ^200)[1,:] )

    not_converging = true
    check = 0
    iter = 0

    while not_converging

        Λᴼˡᵈ, Ψᴼˡᵈ = copy(mdl.Λ), copy(mdl.Ψ)

        # firstly idx = 0
        mdl.Ψ, mdl.Λ = mdl.alternate_return(mdl)

        # secondly idx ≥ 1
        indices = unique(mdl.DR)
        indices = indices[indices .> 0]

        # Threads.@threads
        for idx = indices
            # for idx = LinearIndices(mdl.Ψ[:,1]) # 4 debug only

            mdl.Ψ[idx,:] += vec(sum((tm .== idx) .* Ψᴼˡᵈ * mdl.λ, dims = 1))

        end

        #aftercscalling

        scale = 1/(sum(mdl.Ψ)+sum(mdl.Λ))
        mdl.Ψ *= scale
        mdl.Λ *= scale

        #CONVERGEANCE CRITERION
        con_val = (sum(abs.(mdl.Λ - Λᴼˡᵈ) ) + sum(abs.(mdl.Ψ - Ψᴼˡᵈ)))
        iter +=1

        if (con_val < Ψ_precission)

            not_converging = check < 7
            check += 1
            not_converging || printstyled("\e[s\e[6G✓\e[u",bold = true, color = :green)

        elseif (iter == Ψ_maxiter)

            not_converging = false
            printstyled("\e[s\e[6G!\e[u",bold = true, color = :red)
            global Ψ_punish = true

        end

    end

end


text:
using Mamba

## Data
equiv = Dict{Symbol, Any}(
  :group => [1, 1, 2, 2, 2, 1, 1, 1, 2, 2],
  :y =>
    [1.40 1.65
     1.64 1.57
     1.44 1.58
     1.36 1.68
     1.65 1.69
     1.08 1.31
     1.09 1.43
     1.25 1.44
     1.25 1.39
     1.30 1.52]
)
equiv[:N] = size(equiv[:y], 1)
equiv[:P] = size(equiv[:y], 2)

equiv[:T] = [equiv[:group] 3 - equiv[:group]]

## Model Specification
model = Model(

  y = Stochastic(2,
    (delta, mu, phi, pi, s2_1, T) ->
      begin
        sigma = sqrt(s2_1)
        UnivariateDistribution[
          begin
            m = mu + (-1)^(T[i, j] - 1) * phi / 2 + (-1)^(j - 1) * pi / 2 +
                delta[i, j]
            Normal(m, sigma)
          end
          for i in 1:10, j in 1:2
        ]
      end,
    false
  ),

  delta = Stochastic(2,
    s2_2 -> Normal(0, sqrt(s2_2)),
    false
  ),

  mu = Stochastic(
    () -> Normal(0, 1000)
  ),

  phi = Stochastic(
    () -> Normal(0, 1000)
  ),

  theta = Logical(
    phi -> exp(phi)
  ),

  pi = Stochastic(
    () -> Normal(0, 1000)
  ),

  s2_1 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  ),

  s2_2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  ),

  equiv = Logical(
    theta -> Int(0.8 < theta < 1.2)
  )

)

## Initial Values
inits = [
  Dict(:y => equiv[:y], :delta => zeros(10, 2), :mu => 0, :phi => 0,
       :pi => 0, :s2_1 => 1, :s2_2 => 1),
  Dict(:y => equiv[:y], :delta => zeros(10, 2), :mu => 10, :phi => 10,
       :pi => 10, :s2_1 => 10, :s2_2 => 10)
]

## Sampling Scheme
scheme = [NUTS(:delta),
          Slice([:mu, :phi, :pi], 1.0),
          Slice([:s2_1, :s2_2], 1.0, Univariate)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, equiv, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)


text:
using Revise
import MitosisStochasticDiffEq as MSDE
using Test, Random
using LinearAlgebra
using BenchmarkTools
using MitosisStochasticDiffEq: NoiseGrid, EM
using Mitosis: WGaussian
using StaticArrays

#=
function MSDE.saveit!(uu::Vector, u::Tuple{Int64, Float64, Vector{Float64}}, P)
  push!(uu, (u[1], u[2], copy(u[3])))
end
=#

# define SDE function
f(u,p,t) = p[1]*u + p[2]
g(u,p,t) = p[3]

# time span
tstart = 0.0
tend = 1.0
dt = 0.0001
trange = tstart:dt:tend

B, β, σ̃ = @SMatrix([-0.1]), @SVector([0.2]), 1.3
p = (B, β, σ̃)

u0 = @SVector [1.0]

k_oop = MSDE.SDEKernel(f, g, trange, p)

# pass noise process and compare with EM()
w = zero(u0)
Ws = [w]
for (i,ti) in enumerate(trange[1:end-1])
    push!(Ws, w + sqrt(trange[i+1]-ti)*@SVector(randn(1)))
end

function MSDE.dZ!(u, dz::Tuple{Int64, Float64, <:SVector}, Z, P::MSDE.SDEKernel)
    i = u[1]
    (Z[i+1][1] - Z[i][1], Z[i+1][2] - Z[i][2],  Z[i+1][3] -  Z[i][3])
end

function MSDE.tangent!(du::Tuple{Int64, Float64, <:SVector}, u, dz, P::MSDE.SDEKernel)
    MSDE.@unpack f, g, p, noise_rate_prototype = P
    k1 = f(u[3],p,u[2])
    g1 = g(u[3],p,u[2])
    du3 = k1*dz[2] + g1*dz[3]
    (dz[1], dz[2], du3)
end

function MSDE.exponential_map!(u::Tuple{Int64, Float64, <:SVector}, du, P::MSDE.SDEKernel)
    (u[1] + du[1], u[2] + du[2], u[3] + du[3])
end

sol2, solend2 = @btime MSDE.sample(k_oop, u0, MSDE.EulerMaruyama!(), Ws)
v = solend2[end]

Z = collect(zip(1:length(trange), trange, Ws))
u = (1, 0.0, u0)

sol7, solend7 = @time MSDE.solve!(MSDE.EulerMaruyama!(), typeof(u)[], u, Z, k_oop)
@test solend2[3] ≈ solend7[3] rtol=1e-8

@btime MSDE.solve!(MSDE.EulerMaruyama!(), nothing, u, Z, k_oop);

message1_, backward1 = MSDE.backwardfilter(k_oop, WGaussian{(:μ, :Σ, :c)}(v, @SMatrix([1.0]), 0.0))
U = reinterpret(Tuple{typeof(message1_.sol.u[1].x[1]), typeof(message1_.sol.u[1].x[2]), typeof(message1_.sol.u[1].x[3])}, message1_.soldis)
message1 = MSDE.Message(message1_.ktilde, message1_.sol, U, message1_.ts, message1_.filter)
gp1 = MSDE.GuidedSDE(k_oop, message1)

ξ = (1, 0.0, u0, 0.0)
sol, solend = MSDE.solve!(MSDE.EulerMaruyama!(), typeof(ξ)[], ξ, Z, gp1)

@btime MSDE.solve!(MSDE.EulerMaruyama!(), nothing, ξ, Z, gp1);


text:
@testset "object" begin
    @testset "pyis" begin
        x = pylist()
        y = Py(x)
        z = pylist()
        @test pyis(x, x)
        @test pyis(x, y)
        @test !pyis(x, z)
    end
    @testset "pyrepr" begin
        @test pyrepr(String, pybuiltins.None) == "None"
        @test pyrepr(String, pyint(123)) == "123"
        @test pyrepr(String, "hello") == "'hello'"
    end
    @testset "pyascii" begin
        @test pyascii(String, pybuiltins.None) == "None"
        @test pyascii(String, pyint(234)) == "234"
        @test pyascii(String, "hello") == "'hello'"
    end
    @testset "pyhasattr" begin
        @test pyhasattr(pybuiltins.None, "__class__")
        @test !pyhasattr(pybuiltins.None, "not_an_attr")
    end
    @testset "pygetattr" begin
        @test pyisinstance(pygetattr(pybuiltins.None, "__class__"), pybuiltins.type)
        @test_throws PyException pygetattr(pybuiltins.None, "not_an_attr")
        @test pyisinstance(pygetattr(pybuiltins.None, "__class__", pybuiltins.None), pybuiltins.type)
        @test pyis(pygetattr(pybuiltins.None, "not_an_attr", pybuiltins.None), pybuiltins.None)
    end
    @testset "pysetattr" begin
        x = pytype("Foo", (), ())()
        @test_throws PyException pysetattr(pybuiltins.None, "not_an_attr", 123)
        pysetattr(x, "foo", 123)
        @test pyeq(Bool, pygetattr(x, "foo"), 123)
    end
    @testset "pydelattr" begin
        x = pytype("Foo", (), ())()
        @test !pyhasattr(x, "foo")
        pysetattr(x, "foo", 123)
        @test pyhasattr(x, "foo")
        pydelattr(x, "foo")
        @test !pyhasattr(x, "foo")
        @test_throws PyException pydelattr(pybuiltins.None, "__class__")
        @test_throws PyException pydelattr(pybuiltins.None, "not_an_attr")
    end
    @testset "pyissubclass" begin
        @test pyissubclass(pybuiltins.object, pybuiltins.object)
        @test pyissubclass(pybuiltins.int, pybuiltins.object)
        @test pyissubclass(pybuiltins.bool, pybuiltins.object)
        @test pyissubclass(pybuiltins.bool, pybuiltins.int)
        @test !pyissubclass(pybuiltins.int, pybuiltins.bool)
        @test !pyissubclass(pybuiltins.object, pybuiltins.int)
        @test !pyissubclass(pybuiltins.object, pybuiltins.bool)
        @test_throws PyException pyissubclass(pybuiltins.None, pybuiltins.None)
        @test_throws PyException pyissubclass(pybuiltins.object, pybuiltins.None)
        @test_throws PyException pyissubclass(pybuiltins.None, pybuiltins.object)
    end
    @testset "pyisinstance" begin
        o = pybuiltins.object()
        i = pybuiltins.int()
        b = pybuiltins.bool()
        @test pyisinstance(o, pybuiltins.object)
        @test pyisinstance(i, pybuiltins.object)
        @test pyisinstance(b, pybuiltins.object)
        @test pyisinstance(b, pybuiltins.int)
        @test !pyisinstance(i, pybuiltins.bool)
        @test !pyisinstance(o, pybuiltins.int)
        @test !pyisinstance(o, pybuiltins.bool)
        @test_throws PyException pyisinstance(o, pybuiltins.None)
    end
    @testset "pyhash" begin
        @test pyhash(pybuiltins.None) isa Integer
        @test pyhash(pybuiltins.None) == pyhash(pybuiltins.None)
        @test pyhash(pytuple((1, 2, 3))) == pyhash(pytuple((1, 2, 3)))
    end
    @testset "pytruth" begin
        @test !pytruth(pybuiltins.None)
        @test pytruth(pybuiltins.object)
        for n in -2:2
            @test pytruth(pyint(n)) == (n != 0)
        end
        for n in 0:3
            @test pytruth(pylist(zeros(n))) == (n != 0)
        end
    end
    @testset "pynot" begin
        @test pynot(pybuiltins.None)
        @test !pynot(pybuiltins.int)
        for n in -2:2
            @test pynot(pyfloat(n)) == (n == 0)
        end
        for n in 0:3
            @test pynot(pylist(zeros(n))) == (n == 0)
        end
    end
    @testset "pylen" begin
        for n in 0:3
            @test pylen(pytuple(zeros(n))) == n
            @test pylen(pylist(zeros(n))) == n
        end
        @test_throws PyException pylen(pybuiltins.None)
        @test_throws PyException pylen(pyint(0))
    end
    @testset "pygetitem" begin
        x = pyrange(3)
        @test_throws PyException pygetitem(x, pybuiltins.None)
        @test_throws PyException pygetitem(x, 3)
        @test_throws PyException pygetitem(x, -4)
        @test pyeq(Bool, pygetitem(x, 0), 0)
        @test pyeq(Bool, pygetitem(x, 1), 1)
        @test pyeq(Bool, pygetitem(x, 2), 2)
        @test pyeq(Bool, pygetitem(x, -1), 2)
        @test pyeq(Bool, pygetitem(x, -2), 1)
        @test pyeq(Bool, pygetitem(x, -3), 0)
    end
    @testset "pysetitem" begin
        x = pydict()
        @test_throws PyException pygetitem(x, "foo")
        pysetitem(x, "foo", 123)
        @test pyeq(Bool, pygetitem(x, "foo"), 123)
        pysetitem(x, "foo", 0)
        @test pyeq(Bool, pygetitem(x, "foo"), 0)
    end
    @testset "pydelitem" begin
        x = pydict()
        pysetitem(x, "foo", 123)
        @test pylen(x) == 1
        pydelitem(x, "foo")
        @test pylen(x) == 0
    end
    @testset "pydir" begin
        x = pytype("Foo", (), (foo=1, bar=2))()
        d = pydir(x)
        @test pycontains(d, "__class__")
        @test pycontains(d, "foo")
        @test pycontains(d, "bar")
        @test !pycontains(d, "baz")
    end
    @testset "pycall" begin
        @test pyeq(Bool, pycall(pybuiltins.int), 0)
        @test_throws PyException pycall(pybuiltins.None)
        @test pyeq(Bool, pycall(pybuiltins.int, "10"), 10)
        @test_throws PyException pycall(pybuiltins.int, "10", "20")
        @test pyeq(Bool, pycall(pybuiltins.int, "10", base=16), 16)
        @test_throws PyException pycall(pybuiltins.int, "10", base="16")
        @test_throws PyException pycall(pybuiltins.int, "10", bad_argument=0)
    end
    @testset "pyeq" begin
        @test pyis(pyeq(pybuiltins.None, pybuiltins.None), pybuiltins.True)
        @test pyeq(Bool, pybuiltins.None, pybuiltins.None)
        @test pyis(pyeq(pybuiltins.None, pybuiltins.int), pybuiltins.False)
        @test !pyeq(Bool, pybuiltins.None, pybuiltins.int)
    end
    @testset "pyne" begin
        @test pyis(pyne(pybuiltins.None, pybuiltins.None), pybuiltins.False)
        @test !pyne(Bool, pybuiltins.None, pybuiltins.None)
        @test pyis(pyne(pybuiltins.None, pybuiltins.int), pybuiltins.True)
        @test pyne(Bool, pybuiltins.None, pybuiltins.int)
    end
    @testset "pylt" begin
        for a in -1:1
            for b in -1:1
                @test pyis(pylt(pyint(a), pyint(b)), pybool(a < b))
                @test pylt(Bool, pyint(a), pyint(b)) == (a < b)
            end
        end
    end
    @testset "pyle" begin
        for a in -1:1
            for b in -1:1
                @test pyis(pyle(pyint(a), pyint(b)), pybool(a ≤ b))
                @test pyle(Bool, pyint(a), pyint(b)) == (a ≤ b)
            end
        end
    end
    @testset "pygt" begin
        for a in -1:1
            for b in -1:1
                @test pyis(pygt(pyint(a), pyint(b)), pybool(a > b))
                @test pygt(Bool, pyint(a), pyint(b)) == (a > b)
            end
        end
    end
    @testset "pyge" begin
        for a in -1:1
            for b in -1:1
                @test pyis(pyge(pyint(a), pyint(b)), pybool(a ≥ b))
                @test pyge(Bool, pyint(a), pyint(b)) == (a ≥ b)
            end
        end
    end
    @testset "pycontains" begin
        x = pyrange(3)
        @test pycontains(x, 0)
        @test pycontains(x, 1)
        @test pycontains(x, 2)
        @test !pycontains(x, 3)
        @test !pycontains(x, -1)
        @test !pycontains(x, pybuiltins.None)
    end
    @testset "pyin" begin
        x = pyrange(3)
        @test pyin(0, x)
        @test pyin(1, x)
        @test pyin(2, x)
        @test !pyin(3, x)
        @test !pyin(-1, x)
        @test !pyin(pybuiltins.None, x)
    end
end

@testset "iter" begin
    @test_throws PyException pyiter(pybuiltins.None)
    @test_throws PyException pyiter(pybuiltins.True)
    it = pyiter(pyrange(2))
    x = PythonCall.pynext(it)
    @test !PythonCall.ispynull(x)
    @test pyeq(Bool, x, 0)
    x = PythonCall.pynext(it)
    @test !PythonCall.ispynull(x)
    @test pyeq(Bool, x, 1)
    x = PythonCall.pynext(it)
    @test PythonCall.ispynull(x)
end

@testset "number" begin
    @testset "pyneg" begin
        for n in -2:2
            @test pyeq(Bool, pyneg(pyint(n)), pyint(-n))
        end
    end
    @testset "pypos" begin
        for n in -2:2
            @test pyeq(Bool, pypos(pyint(n)), pyint(n))
        end
    end
    @testset "pyabs" begin
        for n in -2:2
            @test pyeq(Bool, pyabs(pyint(n)), pyint(abs(n)))
        end
    end
    @testset "pyinv" begin
        for n in -2:2
            @test pyeq(Bool, pyinv(pyint(n)), pyint(-n-1))
        end
    end
    @testset "pyindex" begin
        for n in -2:2
            @test pyeq(Bool, pyindex(pyint(n)), pyint(n))
        end
    end
    @testset "pyadd" begin
        for x in -2:2
            for y in -2:2
                @test pyeq(Bool, pyadd(pyint(x), pyint(y)), pyint(x+y))
            end
        end
    end
    @testset "pysub" begin
        for x in -2:2
            for y in -2:2
                @test pyeq(Bool, pysub(pyint(x), pyint(y)), pyint(x-y))
            end
        end
    end
    @testset "pymul" begin
        for x in -2:2
            for y in -2:2
                @test pyeq(Bool, pymul(pyint(x), pyint(y)), pyint(x*y))
            end
        end
    end
    # TODO
    # @testset "pymatmul" begin
    #     for x in -2:2
    #         for y in -2:2
    #             @test pyeq(Bool, pymul(pyint(x), pyint(y)), pyint(x*y))
    #         end
    #     end
    # end
    @testset "pytruediv" begin
        for x in -2:2
            for y in -2:2
                if y == 0
                    @test_throws PyException pytruediv(pyint(x), pyint(y))
                else
                    @test pyeq(Bool, pytruediv(pyint(x), pyint(y)), pyfloat(x/y))
                end
            end
        end
    end
    @testset "pyfloordiv" begin
        for x in -2:2
            for y in -2:2
                if y == 0
                    @test_throws PyException pyfloordiv(pyint(x), pyint(y))
                else
                    @test pyeq(Bool, pyfloordiv(pyint(x), pyint(y)), pyfloat(fld(x, y)))
                end
            end
        end
    end
    @testset "pymod" begin
        for x in -2:2
            for y in -2:2
                if y == 0
                    @test_throws PyException pymod(pyint(x), pyint(y))
                else
                    @test pyeq(Bool, pymod(pyint(x), pyint(y)), pyint(mod(x, y)))
                end
            end
        end
    end
    @testset "pydivmod" begin
        for x in -2:2
            for y in -2:2
                if y == 0
                    @test_throws PyException pydivmod(pyint(x), pyint(y))
                else
                    @test pyeq(Bool, pydivmod(pyint(x), pyint(y)), pytuple(fldmod(x, y)))
                end
            end
        end
    end
    @testset "pylshift" begin
        for n in 0:3
            @test pyeq(Bool, pylshift(pyint(123), pyint(n)), pyint(123 << n))
        end
    end
    @testset "pyrshift" begin
        for n in 0:3
            @test pyeq(Bool, pyrshift(pyint(123), pyint(n)), pyint(123 >> n))
        end
    end
    @testset "pyand" begin
        for x in 0:3
            for y in 0:3
                @test pyeq(Bool, pyand(pyint(x), pyint(y)), pyint(x & y))
            end
        end
    end
    @testset "pyxor" begin
        for x in 0:3
            for y in 0:3
                @test pyeq(Bool, pyxor(pyint(x), pyint(y)), pyint(x ⊻ y))
            end
        end
    end
    @testset "pyor" begin
        for x in 0:3
            for y in 0:3
                @test pyeq(Bool, pyor(pyint(x), pyint(y)), pyint(x | y))
            end
        end
    end
    # TODO: in-place operators
end


text:
export saveMagneticFieldAsHDF5, MagneticFieldSweepCurrentsMeas

#TODO: Unit handling should be put into GaussMeter

# define measObj
@compat struct MagneticFieldSweepCurrentsMeas <: MeasObj
  su::SurveillanceUnit
  daq::AbstractDAQ
  gauss::GaussMeter
  positions::Positions
  currents::Matrix{Float64}
  gaussRanges::Vector{Int}
  waitTime::Float64
#  voltToCurrent::Float64
  pos::Matrix{typeof(1.0Unitful.m)}
  magneticField::Array{typeof(1.0Unitful.T),3}
  magneticFieldError::Array{typeof(1.0Unitful.T),3}

  MagneticFieldSweepCurrentsMeas(su, daq, gauss, positions, currents, gausMeterRanges, waitTime) =
                 new(su, daq, gauss, positions, currents, gausMeterRanges, waitTime, #voltToCurrent,
                      zeros(typeof(1.0Unitful.m),3,length(positions)),
                      zeros(typeof(1.0Unitful.T),3,length(positions),size(currents,2)),
                      zeros(typeof(1.0Unitful.T),3,length(positions),size(currents,2)))
end

# define preMoveAction
function preMoveAction(measObj::MagneticFieldSweepCurrentsMeas,
                       pos::Vector{typeof(1.0Unitful.mm)}, index)
  @info "moving to position" pos
end

# define postMoveAction
function postMoveAction(measObj::MagneticFieldSweepCurrentsMeas,
                        pos::Vector{typeof(1.0Unitful.mm)}, index)
  @info "post action" pos index length(measObj.positions)
  #sleep(0.05)
  measObj.pos[:,index] = pos

  for l=1:size(measObj.currents,2)
    # set current at DC sources
    #value(measObj.rp,"AOUT0",measObj.currents[1,l]*measObj.voltToCurrent)
    #value(measObj.rp,"AOUT1",measObj.currents[2,l]*measObj.voltToCurrent)
    setSlowDAC(measObj.daq, measObj.currents[1,l], 0)
    setSlowDAC(measObj.daq, measObj.currents[2,l], 1)

    @debug "Set DC source $(measObj.currents[1,l]*Unitful.A)  $(measObj.currents[2,l]*Unitful.A)"
    # set measurement range of gauss meter
    range = measObj.gaussRanges[l]
    setAllRange(measObj.gauss, Char("$range"[1]))
    # wait until magnet is on field
    sleep(0.6)
    # perform field measurment
    magneticField = getXYZValues(measObj.gauss)
    measObj.magneticField[:,index,l] = magneticField
    # perform error estimation based on gauss meter specification
    magneticFieldError = zeros(typeof(1.0Unitful.T),3,2)
    magneticFieldError[:,1] = abs.(magneticField)*1e-3
    magneticFieldError[:,2] = getFieldError(range).*ones(3)
    measObj.magneticFieldError[:,index,l] = maximum(magneticFieldError,dims=2)

    @info "Field $(uconvert.(Unitful.mT,measObj.magneticField[:,index,l]))"
  end
  setSlowDAC(measObj.daq, 0.0, 0)
  setSlowDAC(measObj.daq, 0.0, 1)

  @showprogress for i=1:measObj.waitTime
    sleep(1)
  end

end

function saveMagneticFieldAsHDF5(measObj::MagneticFieldSweepCurrentsMeas,
       filename::String, params=Dict{String,Any}())
  h5open(filename, "w") do file
    write(file, measObj.positions)
    write(file, "/unitCoords", "m")
    write(file, "/unitFields", "T")
    write(file, "/positions", ustrip.(measObj.pos))
    write(file, "/fields", ustrip.(measObj.magneticField))
    write(file, "/fieldsError", ustrip.(measObj.magneticFieldError))
    write(file, "/currents", measObj.currents)
    for (key,value) in params
      write(file, key, value)
    end
  end
end