text:
using BinDeps

@BinDeps.setup

glib = library_dependency("glib", aliases = ["libglib-2.0-0", "libglib-2.0", "libglib-2.0.so.0"])
lcm = library_dependency("lcm", aliases=["liblcm", "liblcm.1"], depends=[glib])
lcmgl_client = library_dependency("bot2-lcmgl-client", aliases=["libbot2-lcmgl-client", "libbot2-lcmgl-client.1"], depends=[lcm])

if Sys.isapple()
    using Homebrew
    provides(Homebrew.HB, "glib", glib, os=:Darwin)
    ENV["PKG_CONFIG_PATH"] = get(ENV, "PKG_CONFIG_PATH", "") * ":" * joinpath(Homebrew.prefix(), "lib", "pkgconfig")
    ENV["INCLUDE_PATH"] = get(ENV, "INCLUDE_PATH", "") * ":" * joinpath(Homebrew.prefix(), "include")
end

provides(AptGet, Dict("libglib2.0-dev" => glib))
provides(Yum,
    Dict("glib" => glib))

libbot_dirname = "libbot-cc8d228b50847c4c55e6963b8ee95c237287547f"
provides(Sources,
    URI("https://github.com/RobotLocomotion/libbot/archive/cc8d228b50847c4c55e6963b8ee95c237287547f.zip"),
    lcmgl_client,
    unpacked_dir=libbot_dirname)
provides(Sources,
    URI("https://github.com/lcm-proj/lcm/releases/download/v1.3.1/lcm-1.3.1.zip"),
    lcm)

prefix = joinpath(BinDeps.depsdir(lcmgl_client), "usr")
provides(SimpleBuild,
    (@build_steps begin
        GetSources(lcm)
        @build_steps begin
            ChangeDirectory(joinpath(BinDeps.depsdir(lcm), "src", "lcm-1.3.1"))
            `./configure --prefix=$(prefix) --with-java=no --with-python=no --with-lua=no`
            MakeTargets(".", [])
            MakeTargets(".", ["install"])
        end
    end), lcm)

provides(SimpleBuild,
    (@build_steps begin
        GetSources(lcmgl_client)
        @build_steps begin
            MakeTargets(joinpath(BinDeps.depsdir(lcmgl_client), "src", libbot_dirname), ["BUILD_PREFIX=$(prefix)"])
            @static if Sys.isapple()
                `install_name_tool -change $(joinpath(prefix, "lib", "liblcm.1.dylib")) "@loader_path/liblcm.1.dylib" $(joinpath(prefix, "lib", "libbot2-lcmgl-client.1.dylib"))`
            end
        end
    end), lcmgl_client)

@BinDeps.install Dict(:lcm => :liblcm, "bot2-lcmgl-client" => :libbot2_lcmgl_client)


text:
using HorizonSideRobots
include("lib_for_practic_10.jl")

function cherez_2(r::Robot, side::HorizonSide)
    while !(isborder(r, side))
        putmarker!(r)
        try_move!(r, side)
        try_move!(r,side)
    end
end

function step_verh(r::Robot, side::HorizonSide)
    if ismarker(r)==false
        move!(r, Nord)
    end

    if ismarker(r)==true
        move!(r, Nord)
        move!(r, side)
    end
end

function mark_2(r::Robot)
    while (isborder(r, Nord)==false && isborder(r, Ost)==false)
        cherez_2(r, Ost)
        step_verh(r, West)
        cherez_2(r, West)
        step_verh(r, Ost)
    end
end

mark_2(r)

ext: jl


text:
using PAINTER
using Base.Test

FOV = 0.01
indwvl = 1:30
nx = 64
eps1 = 1e-4
eps2 = 1e-4
rho_y = 10.
alpha = 1e4
beta = 1e5
rho_spat = 4.
rho_ps = rho_spat
lambda_spat = 1e-5
rho_spec = 0.5
lambda_spec = 1e-5
nbitermax = 4
dptype = "sliding"
dpprm = 5
aff=false

OIDATA, PDATA = painter(nbitermax = nbitermax, nx = nx, lambda_spat = lambda_spat,
                            lambda_spec = lambda_spec, rho_y = rho_y, rho_spat = rho_spat,
                            rho_spec = rho_spec, rho_ps = rho_ps, alpha = alpha, beta = beta,
                            eps1 = eps1, eps2 = eps2, FOV = FOV, indwvl = indwvl,
                            dptype = dptype, dpprm = dpprm, aff=aff)

@test_approx_eq_eps PDATA.crit1[1] 6.25337 1E-4


text:
include("test_sti.jl")
include("test_solver.jl")


text:
struct Check
    checkfn
    name
    throw_error::Bool
    message
end

struct SanityCheckException <: Exception end

Base.show(io::IO, check::Check) = print(io, "Check(\"", check.name, "\")")

function runchecks(checks, learner)
    failed = trues(length(checks))
    for (i, check) in enumerate(checks)
        failed[i] = !check.checkfn(learner)
    end

    failedchecks = checks[failed]

    if isempty(failedchecks)
        return
    end

    println("$(length(failedchecks))/$(length(checks)) sanity checks failed:")

    for (i, check) in enumerate(failedchecks)
        println("---")
        println(i, ": ", check.name, " (", check.throw_error ? "ERROR" : "WARNING", ")")
        println()
        println(check.message)
    end

    if any(getfield.(failedchecks, :throw_error))
        throw(SanityCheckException())
    end
end

"""
    SanityCheck([checks; usedefault = true])

Callback that runs sanity [`Check`](#)s when the `Learner` is initialized.
If `usedefault` is `true`, it will run all checks in FluxTraining.CHECKS
in addition to the ones you pass in.
"""
struct SanityCheck <: Callback
    checks::Vector{Check}
    function SanityCheck(checks = []; usedefault = isempty(checks))
        if usedefault
            checks = vcat(checks, CHECKS)
        end
        return new(checks)
    end
end

stateaccess(::SanityCheck) = (
    data = Read(),
    model = Read(),
    lossfn = Read(),
    optimizer = Read(),
    callbacks = Read(),
)

function on(::Init, phase::AbstractTrainingPhase, cb::SanityCheck, learner)
    runchecks(cb.checks, learner)
end

# Checks

CheckDataIteratorTrain() =
    Check(
        "Has training data iterator",
        true,
"""
`learner` training data iterator is `nothing` (not set). You can pass it to
`Learner` in one of the following ways:

- `Learner(model, traindataiter, opt, lossfn)`
- `Learner(model, (traindataiter, valdataiter), opt, lossfn)`
- `Learner(model, (training = traindataiter,), opt, lossfn)`
"""
        ) do learner
        !isnothing(learner.data.training)
    end

CheckDataIteratorValid() =
    Check(
        "Has validation data iterator",
        false,
"""
`learner` validation data iterator is `nothing` (not set).
It is not mandatory, but you must set it to fit `ValidationPhase`s.

You can pass it to `Learner` like this:

- `Learner(model, (traindataiter, valdataiter), opt, lossfn)`

Or if you want to use training data as validation data:

- `Learner(model, (traindataiter, traindataiter), opt, lossfn)`
"""
        ) do learner
        !isnothing(learner.data.validation)
    end

CheckIteratesTuples() =
    Check(
        "Data iterators iterate over tuples",
        true,
"""
Data iterators need to be iterable and return tuples.
This means that `for (x, y) in dataiter end` works where
`(x, y)` is a pair of model inputs and outputs.
"""
        ) do learner
        try
            batch, _ = Base.iterate(learner.data.training)
            x, y = batch
        catch
            return false
        end
        return true
    end

CheckModelLossStep() =
    Check(
        "Model and loss function compatible with data",
        true,
"""
To perform the optimization step, the model and loss function need
to be compatible with the data. This means the following must work:

- `(x, y), _ = iterate(learner.data.training)`
- `ŷ = learner.model(x)`
- `loss = learner.lossfunction(ŷ, y)`
"""
        ) do learner
        try
            dev = ToGPU() in learner.callbacks.cbs ? gpu : identity
            x, y = dev(iterate(learner.data.training)[1])
            ŷ = dev(learner.model)(x)
            @assert learner.lossfn(ŷ, y) isa Number
        catch
            return false
        end
        return true
    end

const CHECKS = [
    CheckDataIteratorTrain(),
    CheckDataIteratorValid(),
    CheckIteratesTuples(),
    CheckModelLossStep(),
]


text:
module RVPosterior

using AR1
using Ensemble
using Kepler

immutable PhysicalParams
    mu::Float64
    sigma::Float64
    tau::Float64
    nu::Float64
    betas::Array{Float64, 1}
    pers::Array{Float64, 1}
    vels::Array{Float64, 1}
    ecws::Array{Float64, 1}
    esws::Array{Float64, 1}
    chis::Array{Float64, 1}
end

nplanets(p::PhysicalParams) = size(p.pers,1)
npred(p::PhysicalParams) = size(p.betas,1)

ecc(p::PhysicalParams) = sqrt(p.ecws.*p.ecws + p.esws.*p.esws)
omega(p::PhysicalParams) = atan2(p.esws, p.ecws)

function rvs(p::PhysicalParams, ts::Array{Float64, 1})
    es = ecc(p)
    omegas = omega(p)

    rvs = Kepler.rv_model(ts, p.vels, es, omegas, p.chis, p.pers)
    sum(rvs, 2)
end

function cresid(p, ts, vs, preds)
    npl = size(p.pers, 1)
    rs = zeros(ts)

    if npl > 0
        rv = rvs(p, ts)
    else
        rv = zeros(ts)
    end

    for i in eachindex(rs)
        rs[i] = vs[i] - rv[i]
        for j in 1:size(preds,2)
            rs[i] -= p.betas[j]*preds[i,j]
        end
    end

    rs
end

function uncorr_resid(p, ts, vs, dvs, preds)
    crs = cresid(p, ts, vs, preds)

    AR1.residuals(AR1.AR1Model(p.mu, p.sigma, p.tau, p.nu), ts, crs, dvs)
end

function loglike(p, ts, vs, dvs, preds)
    rs, srs = uncorr_resid(p, ts, vs, dvs, preds)

    lp = 0.0
    for i in 1:size(ts, 1)
        lp += -0.91893853320467274180 - log(srs[i]) - 0.5*rs[i]*rs[i]/(srs[i]*srs[i])
    end
    lp
end

immutable Posterior
    ts::Array{Float64, 1}
    vs::Array{Float64, 1}
    dvs::Array{Float64, 1}
    preds::Array{Float64, 2}
    npl::Int
end

function toparams(post, v)
    nplp = post.npl*5
    nbet = size(v,1) - 4 - nplp

    mu = v[1]
    sigma = exp(v[2])
    tau = exp(v[3])
    nu = exp(v[4])

    betas = v[5:5+nbet-1]

    i = 5+nbet
    Ps = zeros(post.npl)
    Ks = zeros(post.npl)
    ecws = zeros(post.npl)
    esws = zeros(post.npl)
    chis = zeros(post.npl)
    for j in 1:post.npl
        Ps[j] = exp(v[i])
        i += 1

        Ks[j] = exp(v[i])
        i += 1

        x = v[i:i+1]
        z = Ensemble.Parameterizations.unit_disk_value(x)
        ecws[j] = z[1]
        esws[j] = z[2]
        i += 2

        chis[j] = Ensemble.Parameterizations.bounded_value(v[i], 0.0, 1.0)
        i += 1
    end

    PhysicalParams(mu, sigma, tau, nu, betas, Ps, Ks, ecws, esws, chis)
end

function tov(post, p)
    nplp = post.npl*5
    nbet = size(p.betas, 1)

    n = 4 + nplp + nbet

    v = zeros(n)
    i = 1

    v[i] = p.mu
    i += 1

    v[i] = log(p.sigma)
    i += 1

    v[i] = log(p.tau)
    i += 1

    v[i] = log(p.nu)
    i += 1

    v[i:i+nbet-1] = p.betas
    i += nbet

    for j in 1:post.npl
        v[i] = log(p.pers[j])
        i += 1

        v[i] = log(p.vels[j])
        i += 1

        z = Float64[p.ecws[j], p.esws[j]]
        x = Ensemble.Parameterizations.unit_disk_param(z)
        v[i:i+1] = x
        i += 2

        v[i] = Ensemble.Parameterizations.bounded_param(p.chis[j], 0.0, 1.0)
        i += 1
    end

    v
end

function logprior(post, p, v)
    nplp = post.npl*5
    nbet = size(p.betas, 1)

    T = post.ts[end] - post.ts[1]
    vmin = minimum(post.vs)
    vmax = maximum(post.vs)

    V = vmax - vmin

    lp = 0.0

    # Jumping in log(sigma), log(tau), log(nu), but want flat priors
    # p(logsigma) = sigma*p(sigma)
    lp += log(p.sigma)
    lp += log(p.tau)
    lp += log(p.nu)

    if post.npl > 0
        for i in 1:post.npl
            # Similarly, want flat priors on P and velocity
            lp += log(p.pers[i])
            lp += log(p.vels[i])

            z = Float64[p.ecws[i], p.esws[i]]
            j = 5 + nbet + 2 + (i-1)*5
            ve = v[j:j+1]
            lp += Ensemble.Parameterizations.unit_disk_logjac(z, ve)

            chiv = v[5 + nbet + 4 + (i-1)*5]
            lp += Ensemble.Parameterizations.bounded_logjac(p.chis[i], chiv, 0.0, 1.0)
        end
    end

    lp
end

function lnprob(post, v)
    p = toparams(post, v)

    lp = logprior(post, p, v)

    if lp == -Inf
        return lp
    else
        return loglike(p, post.ts, post.vs, post.dvs, post.preds) + lp
    end
end

function addplanet(post, ps::Array{Float64, 2}, P0, K0)
    new_post = Posterior(post.ts, post.vs, post.dvs, post.preds, post.npl+1)

    nwalk = size(ps, 2)

    new_ps = zeros(size(ps,1)+5, nwalk)
    for i in 1:nwalk
        p = toparams(new_post, cat(1, ps[:,i], zeros(5)))
        p.pers[end] = P0*(1.0 + 1e-5*randn())
        p.vels[end] = K0*(1.0 + 1e-5*randn())
        p.ecws[end] = 0.1*rand()
        p.esws[end] = 0.1*rand()
        p.chis[end] = rand()

        new_ps[:,i] = tov(new_post, p)
    end

    lnps = Float64[lnprob(new_post, new_ps[:,i]) for i in 1:nwalk]

    ibest = indmax(lnps)
    pbest = new_ps[:,ibest]

    for i in 1:nwalk
        new_ps[:,i] = pbest + 1e-5*randn(size(pbest,1))
    end

    new_post, new_ps, lnps
end

function save_samples(io, ps)
    writedlm(io, reshape(permutedims(ps, [1, 3, 2]), (size(ps,1), size(ps,2)*size(ps,3))))
end

function load_samples(io, nwalk)
    ps = readdlm(io)

    ps = reshape(ps, size(ps,1), div(size(ps,2),nwalk), nwalk)

    permutedims(ps, [1, 3, 2])
end

end

max_stars_repo_name: farr/RVChallenge

max_issues_repo_name: farr/RVChallenge

max_forks_repo_name: farr/RVChallenge


text:
abstract AbstractChain

chain{N<:Number,T<:AbstractFloat}(s::Symbol,nparas::Int,nsamples::Int,::Type{N},::Type{T},args...;keyargs...) = _chain(Val{s},nparas,nsamples,N,T,args...;keyargs...)

numparas(c::AbstractChain) = size(c.values,1)
numsamples(c::AbstractChain) = size(c.values,2)

accepted!(c::AbstractChain,indicator::AbstractIndicatorMatrix) = c.accepted += accepted(indicator)

samples(c::AbstractChain) = c.values
loglikelihood(c::AbstractChain) = c.loglikelihood
logprior(c::AbstractChain,m::AbstractModel) = logprior(parameters(m),c.values,eltype(c.loglikelihood))
logposterior(c::AbstractChain,m::AbstractModel) = c.loglikelihood + logprior(c,m)


text:
using PyPlot

using SubmodularMaximization

discretization = 1000

num_agents = 20
num_sensors = 6
nominal_area = 1.0

max_success_probability = 0.9
sensor_radius = sqrt(nominal_area / (num_agents * pi))
station_radius = 10 * sensor_radius

sensor = ProbabilisticSensor([0.5, 0.5], sensor_radius, max_success_probability)

density = [detection_probability(sensor, [x,y])
           for x in range(0, stop=1, length=discretization),
           y in range(0, stop=1, length=discretization)]

figure()
imshow(density', cmap="viridis", vmin=minimum(density), vmax=maximum(density[:]),
       extent=[0, 1, 0, 1], interpolation="nearest", origin="lower")
colorbar()

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
module Mex

using MxArrays

export jl_mex, input, call_matlab, is_interrupt_pending, check_for_interrupt

# *** ccall stuff ***

const libmex = open_matlab_library("libmex")
const _call_matlab_with_trap = Libdl.dlsym(libmex, :mexCallMATLABWithTrap)
const libut = open_matlab_library("libut")
const _is_interrupt_pending = Libdl.dlsym(libut, :utIsInterruptPending)

# *** interrupt handling ***

is_interrupt_pending() = ccall(_is_interrupt_pending, UInt8, ()) != 0

function check_for_interrupt()
    if is_interrupt_pending()
        throw(InterruptException())
    end
end

# *** exception handling ***

type MatlabException <: Exception
    ptr::Ptr{Void}
    MatlabException(ptr::Ptr{Void}) = new(ptr)
end

function MatlabException(id::String, msg::String)
    mxid = "jl:"*replace(id, r"[^\w\s]", "_")
    mx = call_matlab(1, "MException", MxArray[MxArray(mxid), MxArray(msg)])[1]
    return MatlabException(mx.ptr)
end

MatlabException(mexn::MatlabException) = mexn

function MatlabException(exn)
    buf = IOBuffer()
    showerror(buf, exn)
    seek(buf, 0)
    msg = replace(readstring(buf), "\\", "\\\\")
    return MatlabException(string(typeof(exn)), msg)
end

function add_cause(mexn::MatlabException, cause::MatlabException)
    args = MxArray[MxArray(x.ptr) for x in [mexn, cause]]
    mx = call_matlab(1, "addCause", args)[1]
    return MatlabException(mx.ptr)
end

function add_backtrace(exn, bt)
    buf = IOBuffer()
    print(buf, "Julia backtrace:")
    Base.show_backtrace(buf, bt)
    seek(buf, 0)
    msg = replace(readstring(buf), "\\", "\\\\")
    cause = MatlabException("backtrace", msg)
    return add_cause(MatlabException(exn), cause)
end

# *** calling into MATLAB ***

# Call a matlab function specified by name
# This version allows full control over data marshaling.
function call_matlab(nout::Integer, fn::String, args::Vector{MxArray})
    ins = Ptr{Void}[arg.ptr for arg in args]
    nin = length(ins)
    outs = Vector{Ptr{Void}}(nout)

    ptr = ccall(_call_matlab_with_trap, Ptr{Void},
        (Int32, Ptr{Ptr{Void}}, Int32, Ptr{Ptr{Void}}, Ptr{UInt8}),
        nout, outs, nin, ins, fn)
    if ptr != C_NULL
        throw(MatlabException(ptr))
    end

    return MxArray[MxArray(o) for o in outs]
end

# This version uses default data marshaling.
function call_matlab(nout, fn, args...)
    mxin = MxArray[MxArray(arg) for arg in args]
    return map(jvalue, call_matlab(nout, fn, mxin))
end

# Make MxArray callable. Works for strings or function handles.
# This version allows full control over data marshaling.
function (mx::MxArray)(nout::Integer, args::Vector{MxArray})
    _args = MxArray[mx]
    append!(_args, args)
    return call_matlab(nout, "feval", _args)
end

# This version uses default data marshaling and assumes exactly one return value.
function (mx::MxArray)(args...)
    mxs = MxArray[MxArray(arg) for arg in args]
    return jvalue(mx(1, mxs)[1])
end

# *** stdout/stderr redirection ***

function fwrite(fid, msg)
    call_matlab(0, "fwrite", convert(Float64, fid), msg, "char")
    return
end

function readloop(stream, fid)
    try
        while isopen(stream)
            fwrite(fid, readavailable(stream))
        end
    end
end

const mexstdout = redirect_stdout()[1]
const mexstderr = redirect_stderr()[1]
@schedule readloop(mexstdout, 1)
@schedule readloop(mexstderr, 2)

# *** stdin ***

input(prompt::String="julia> ") = call_matlab(1, "input", prompt, "s")[1]

# the entry point for calling into julia from matlab
global jl_mex_call_depth = 0
function jl_mex(outs::Vector{Ptr{Void}}, ins::Vector{Ptr{Void}})
    global jl_mex_call_depth += 1
    try
        if jl_mex_call_depth == 1
            jl_mex_outer(outs, ins)
        else
            jl_mex_inner(outs, ins)
        end
    finally
        jl_mex_call_depth -= 1
    end
end

function jl_mex_outer(outs::Vector{Ptr{Void}}, ins::Vector{Ptr{Void}})
    @sync begin
        # do the desired computation
        main_task = @async jl_mex_inner(outs, ins)

        # check for interrupt
        @async while(!istaskdone(main_task))
            if is_interrupt_pending()
                Base.throwto(main_task, InterruptException())
            end
            yield()
        end
    end
end

function jl_mex_inner(outs::Vector{Ptr{Void}}, ins::Vector{Ptr{Void}})
    nouts = length(outs)
    none = MxArray(false)
    for ix in 1:nouts
        outs[ix] = none.ptr
    end
    try
        args = [MxArray(arg) for arg in ins]
        vals = eval(Main, parse(jvalue(args[1])))(args[2:end])
        outix = 2
        for val in vals
            if outix > nouts
                break
            end
            mx = MxArray(val)
            outs[outix] = mx.ptr
            outix += 1
        end
    catch exn
        outs[1] = add_backtrace(exn, catch_backtrace()).ptr
    end
end

# evaluate Julia expressions
jl_eval(exprs::Vector{MxArray}) = [eval(Main, parse(jvalue(e))) for e in exprs]

# Call a julia function, possibly with keyword arguments.
#
# The values in the args array are interpreted as follows:
#      Index = Meaning
# -----------------------------------------------------------
#          1 = the function to call
#          2 = an integer, npos, the number of positional arguments
#   3:2+npos = positional arguments
# 3+npos:end = keyword arguments, in keyword/value pairs
#
# If npos < 0, all arguments are assumed positional.
function jl_call_kw(args::Vector{MxArray})
    vals = map(jvalue, args)
    nvals = length(vals)
    expr = Expr(:call, parse(vals[1]))
    npos = vals[2]
    if npos < 0
        npos = nvals - 2
    end
    for ix in 3:(2+npos)
        push!(expr.args, vals[ix])
    end
    for ix in (3+npos):2:nvals
        push!(expr.args, Expr(:kw, parse(vals[ix]), vals[ix+1]))
    end
    return [eval(Main, expr)]
end

end # module


text:
module Datasets

import BinDeps: unpack_cmd

# UCR datasets #
export ucr_traindata, ucr_testdata, download_ucr

# synthetic datasets #
export fakedata

# place to store the data
const DATAPATH = joinpath(@__DIR__, "../../data") |> normpath

include("ucr_datasets.jl")
include("fake_datasets.jl")

end # DATASETS module


text:
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: February, 2020
# Copyright: Georgia Institute of Technology, 2020
#
# Test gradient

using LinearAlgebra, Statistics, PyPlot, Test
using Random; Random.seed!(1)
push!(LOAD_PATH, "./src/"); using WaveOpsAVP

# Model
nz = 201
dz = R(5)
m_geom = ModelGeom(nz, dz)
nb = 32
mtrue = R(1/2000^2)*ones(R, nz, nb)
mtrue[101-50:101+50, :] .= R(1/2200^2)
Mtrue = Model(m_geom, mtrue)

# Data parameters
nt = 1001
dt = R(0.001)
θ = range(R(0), stop = R(sqrt(3)/2), length = 11)
p = sin.(θ)*sqrt(max(mtrue...))
dat_pars = DataParams(nt, dt, p)

# Source parameters
wav = cfreq_wavelet(50, nt, dt, (R(5), R(10), R(30), R(60)))
src_pars = SrcParams(wav, dt)

# ABC parameters
abc_size = (50, 50)
abc_fact = R(1e-4)
abc_geom = ABC_Geom(nz, abc_size, abc_fact)

# Modeling operator
acc_mod = 4
F = seismod_fun(dat_pars, src_pars; abc_geom = abc_geom, acc = acc_mod)

# Generate synthetics
dat = F.eval(Mtrue)

## Log-posterior

# Setting objective functional
σ2_dat = mean(norm(dat).^R(2))
mpr = randn(R, size(mtrue)); mpr = norm(mtrue, Inf)*mpr/norm(mpr, Inf)
Mpr = Model(m_geom, mpr)
σ2_pr = mean(norm(Mpr).^R(2))
D = derivative_linop(m_geom)
fun(m::Array{R, 2}; compute_grad::Bool = false) = mlogpost_gauss(dat, Model(m_geom, m), σ2_dat, σ2_pr, F; m0 = Mpr, D = D, compute_grad = compute_grad)

# Computing gradient
m0 = R(1/2000^2)*ones(R, nz, nb)
f0, g0 = fun(m0; compute_grad = true)

# Finite-difference approx
dm = randn(R, size(m0))
dm = norm(m0, Inf)*dm/norm(dm, Inf)
t = R(1e-4)
a = (-R(1/12)*fun(m0+2*t*dm)+R(2/3)*fun(m0+t*dm)-R(2/3)*fun(m0-t*dm)+R(1/12)*fun(m0-2*t*dm))/t

# Analytic
b = dot(g0, dm)

# Relative error
@test isapprox(a, b, rtol = R(1e-6))


text:
Test.@testset "assert.jl" begin
    Test.@test_nowarn CompatHelper.always_assert(true)
    Test.@test CompatHelper.always_assert(true) isa Nothing
    Test.@test CompatHelper.always_assert(true) == nothing
    Test.@test Test.@test_nowarn CompatHelper.always_assert(true) isa Nothing
    Test.@test Test.@test_nowarn CompatHelper.always_assert(true) == nothing
    Test.@test_throws CompatHelper.AlwaysAssertionError CompatHelper.always_assert(false)
end

Test.@testset "ci_service.jl" begin
    withenv("GITHUB_REPOSITORY" => "foo/bar") do
        Test.@test CompatHelper.auto_detect_ci_service() isa CompatHelper.CIService
        Test.@test CompatHelper.auto_detect_ci_service() isa CompatHelper.GitHubActions
        Test.@test_throws ArgumentError CompatHelper.main(; keep_existing_compat = false, drop_existing_compat = false)
    end
    withenv("GITHUB_REPOSITORY" => nothing) do
        Test.@test_throws ErrorException CompatHelper.auto_detect_ci_service()
    end
end

Test.@testset "git.jl" begin
    Test.@test CompatHelper.git_decide_master_branch(CompatHelper.DefaultBranch(), "  abcdefg  ") == "abcdefg"
    Test.@test CompatHelper.git_decide_master_branch(" abc ", " xyz ") == "abc"
end

Test.@testset "new_versions.jl" begin
    Test.@test_throws ArgumentError CompatHelper.old_compat_to_new_compat("", "", :abc)
end

Test.@testset "pull_requests.jl" begin
    a = GitHub.Repo(Dict("name" => "Foo", "full_name" => "Foo", "owner" => "Foo", "id" => 1, "url" => "Foo", "html_url" => "Foo", "fork" => false))
    b = GitHub.Repo(Dict("name" => "Foo", "full_name" => "Foo", "owner" => "Foo", "id" => 1, "url" => "Foo", "html_url" => "Foo", "fork" => false))
    c = GitHub.Repo(Dict("name" => "Foo", "full_name" => "Foo", "owner" => "Foo", "id" => 1, "url" => "Foo", "html_url" => "Foo", "fork" => true))
    Test.@test !CompatHelper._repos_are_the_same(a, nothing)
    Test.@test !CompatHelper._repos_are_the_same(nothing, a)
    Test.@test !CompatHelper._repos_are_the_same(nothing, nothing)
    Test.@test CompatHelper._repos_are_the_same(a, a)
    Test.@test CompatHelper._repos_are_the_same(a, b)
    Test.@test !CompatHelper._repos_are_the_same(a, c)
    Test.@test CompatHelper._repos_are_the_same(b, a)
    Test.@test CompatHelper._repos_are_the_same(b, b)
    Test.@test !CompatHelper._repos_are_the_same(b, c)
    Test.@test !CompatHelper._repos_are_the_same(c, a)
    Test.@test !CompatHelper._repos_are_the_same(c, b)
    Test.@test CompatHelper._repos_are_the_same(c, c)
end

Test.@testset "utils.jl" begin
    Test.@test CompatHelper.generate_pr_title_parenthetical(:keep, true) == " (keep existing compat)"
    Test.@test CompatHelper.generate_pr_title_parenthetical(:drop, true) == " (drop existing compat)"
    Test.@test CompatHelper.generate_pr_title_parenthetical(:brandnewentry, true) == ""
    Test.@test CompatHelper.generate_pr_title_parenthetical(:keep, false) == ""
    Test.@test CompatHelper.generate_pr_title_parenthetical(:drop, false) == ""
    Test.@test CompatHelper.generate_pr_title_parenthetical(:brandnewentry, false) == ""
    project = Dict{Any, Any}()
    Test.@test !haskey(project, "compat")
    CompatHelper.add_compat_section!(project)
    Test.@test haskey(project, "compat")
    CompatHelper.add_compat_section!(project)
    Test.@test haskey(project, "compat")
end

Test.@testset "version_numbers.jl" begin
    Test.@test CompatHelper.generate_compat_entry(v"1.2.3") == "1.2"
    Test.@test CompatHelper.generate_compat_entry(v"1.2.0") == "1.2"
    Test.@test CompatHelper.generate_compat_entry(v"1.0.3") == "1.0"
    Test.@test CompatHelper.generate_compat_entry(v"1.0.0") == "1.0"
    Test.@test CompatHelper.generate_compat_entry(v"0.2.3") == "0.2"
    Test.@test CompatHelper.generate_compat_entry(v"0.2.0") == "0.2"
    Test.@test CompatHelper.generate_compat_entry(v"0.0.3") == "0.0.3"
    Test.@test CompatHelper.generate_compat_entry(v"0.0.0") == "0.0.0"

    Test.@test CompatHelper._remove_trailing_zeros(v"11.22.33") == "11.22.33"
    Test.@test CompatHelper._remove_trailing_zeros(v"11.22.0") == "11.22"
    Test.@test CompatHelper._remove_trailing_zeros(v"11.0.33") == "11.0.33"
    Test.@test CompatHelper._remove_trailing_zeros(v"11.0.0") == "11"
    Test.@test CompatHelper._remove_trailing_zeros(v"0.22.33") == "0.22.33"
    Test.@test CompatHelper._remove_trailing_zeros(v"0.22.0") == "0.22"
    Test.@test CompatHelper._remove_trailing_zeros(v"0.0.33") == "0.0.33"
    Test.@test_throws DomainError CompatHelper._remove_trailing_zeros(v"0.0.0")
end


text:
function Tide(arg0::jint)
    return Tide((jint,), arg0)
end

function Tide(arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::jint, arg5::jint)
    return Tide((jint, jint, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5)
end

function equals(obj::Tide, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_delaunay_multipliers(obj::Tide)
    return jcall(obj, "getDelaunayMultipliers", Vector{jint}, ())
end

function get_doodson_multipliers(obj::Tide)
    return jcall(obj, "getDoodsonMultipliers", Vector{jint}, ())
end

function get_doodson_number(obj::Tide)
    return jcall(obj, "getDoodsonNumber", jint, ())
end

function get_phase(obj::Tide, arg0::BodiesElements)
    return jcall(obj, "getPhase", jdouble, (BodiesElements,), arg0)
end

function get_rate(obj::Tide, arg0::BodiesElements)
    return jcall(obj, "getRate", jdouble, (BodiesElements,), arg0)
end

function get_tau_multiplier(obj::Tide)
    return jcall(obj, "getTauMultiplier", jint, ())
end

function hash_code(obj::Tide)
    return jcall(obj, "hashCode", jint, ())
end


text:
using PlanktonIndividuals.Fields
using PlanktonIndividuals.Biogeochemistry

function test_fields()
    grid = RectilinearGrid(size = (4,6,2), x = (0,12), y = (0,12), z = (0,-8))

    nut = nutrients_init(CPU(), grid)

    @test Tuple(collect(keys(nut))) == nut_names
    @test nut.DIC.data == zeros(8,10,6)
    @test interior(nut.DIC.data, grid) == zeros(4,6,2)

    nuts = generate_nutrients(CPU(), grid, default_nut_init())
    @test maximum(nuts.DIC.data) < 23.0
    @test minimum(nuts.DIC.data) > 17.0

    zero_fields!(nuts)
    @test nuts.DIC.data == zeros(8,10,6)

    return nothing
end
function test_fill_halos()
    grid = RectilinearGrid(size = (4,6,2), x = (0,12), y = (0,12), z = (0,-8))
    nuts = generate_nutrients(CPU(), grid, default_nut_init())
    Nx,Ny,Nz = grid.Nx, grid.Ny, grid.Nz
    Hx,Hy,Hz = grid.Hx, grid.Hy, grid.Hz

    PlanktonIndividuals.Fields.fill_halo_east!(nuts.DIC.data, Hx, Nx, Periodic())
    @test nuts.DIC.data[Nx+Hx+1:Nx+2*Hx, :, :] == nuts.DIC.data[1+Hx:2*Hx, :, :]

    PlanktonIndividuals.Fields.fill_halo_west!(nuts.DIC.data, Hx, Nx, Periodic())
    @test nuts.DIC.data[1:Hx, :, :] == nuts.DIC.data[Nx+1:Nx+Hx, :, :]

    PlanktonIndividuals.Fields.fill_halo_north!(nuts.DIC.data, Hy, Ny, Periodic())
    @test nuts.DIC.data[:, Ny+Hy+1:Ny+2*Hy, :] == nuts.DIC.data[:, 1+Hy:2*Hy, :]

    PlanktonIndividuals.Fields.fill_halo_south!(nuts.DIC.data, Hy, Ny, Periodic())
    @test nuts.DIC.data[:, 1:Hy, :] == nuts.DIC.data[:, Ny+1:Ny+Hy, :]

    PlanktonIndividuals.Fields.fill_halo_bottom!(nuts.DIC.data, Hz, Nz, Periodic())
    @test nuts.DIC.data[:, :, Nz+Hz+1:Nz+2*Hz] == nuts.DIC.data[:, :, 1+Hz:2*Hz]

    PlanktonIndividuals.Fields.fill_halo_top!(nuts.DIC.data, Hz, Nz, Periodic())
    @test nuts.DIC.data[:, :, 1:Hz] == nuts.DIC.data[:, :, Nz+1:Nz+Hz]

    PlanktonIndividuals.Fields.fill_halo_east!(nuts.DIC.data, Hx, Nx, Bounded())
    for i in 1:Hx
        @test nuts.DIC.data[Nx+Hx+i, :, :] == nuts.DIC.data[Nx+Hx, :, :]
    end

    PlanktonIndividuals.Fields.fill_halo_west!(nuts.DIC.data, Hx, Nx, Bounded())
    for i in 1:Hx
        @test nuts.DIC.data[i, :, :] == nuts.DIC.data[Hx+1, :, :]
    end

    PlanktonIndividuals.Fields.fill_halo_north!(nuts.DIC.data, Hy, Ny, Bounded())
    for i in 1:Hy
        @test nuts.DIC.data[:, Ny+Hy+i, :] == nuts.DIC.data[:, Ny+Hy, :]
    end

    PlanktonIndividuals.Fields.fill_halo_south!(nuts.DIC.data, Hy, Ny, Bounded())
    for i in 1:Hy
        @test nuts.DIC.data[:, i, :] == nuts.DIC.data[:, Hy+1, :]
    end

    PlanktonIndividuals.Fields.fill_halo_bottom!(nuts.DIC.data, Hz, Nz, Bounded())
    for i in 1:Hz
        @test nuts.DIC.data[:, :, Nz+Hz+i] == nuts.DIC.data[:, :, Nz+Hz]
    end

    PlanktonIndividuals.Fields.fill_halo_top!(nuts.DIC.data, Hz, Nz, Bounded())
    for i in 1:Hz
        @test nuts.DIC.data[:, :, i] == nuts.DIC.data[:, :, Hz+1]
    end

    PlanktonIndividuals.Fields.fill_halo_east_vel!(nuts.DIC.data, Hx, Nx, Bounded())
    for i in 1:Hx-1
        @test nuts.DIC.data[Nx+Hx+1+i, :, :] == nuts.DIC.data[Nx+Hx+1, :, :]
    end

    PlanktonIndividuals.Fields.fill_halo_north_vel!(nuts.DIC.data, Hy, Ny, Bounded())
    for i in 1:Hy-1
        @test nuts.DIC.data[:, Ny+Hx+1+i, :] == nuts.DIC.data[:, Ny+Hy+1, :]
    end

    PlanktonIndividuals.Fields.fill_halo_bottom_vel!(nuts.DIC.data, Hz, Nz, Bounded())
    for i in 1:Hz-1
        @test nuts.DIC.data[:, :, Nz+Hx+1+i] == nuts.DIC.data[:, :, Nz+Hz+1]
    end

    PlanktonIndividuals.Fields.fill_halo_east_Gc!(nuts.DIC.data, Hx, Nx, Bounded())
    for i in 1:Hx
        @test nuts.DIC.data[Nx+Hx+i, :, :] == zeros(10,6)
    end

    PlanktonIndividuals.Fields.fill_halo_north_Gc!(nuts.DIC.data, Hy, Ny, Bounded())
    for i in 1:Hy
        @test nuts.DIC.data[:, Ny+Hx+i, :] == zeros(8,6)
    end

    PlanktonIndividuals.Fields.fill_halo_bottom_Gc!(nuts.DIC.data, Hz, Nz, Bounded())
    for i in 1:Hz
        @test nuts.DIC.data[:, :, Nz+Hx+i] == zeros(8,10)
    end

    return nothing
end

function test_boundary_conditions()
    grid = RectilinearGrid(size = (4, 4, 4), x = (0,32), y = (0,32), z = (0,-32))
    model = PlanktonModel(CPU(), grid; mode = CarbonMode())
    set_bc!(model, :DIC, :west, 0.1)
    @test model.nutrients.DIC.bc.west == 0.1
    set_bc!(model, :DIC, :west, ones(4,4))
    @test model.nutrients.DIC.bc.west == ones(4,4)
    set_bc!(model, :DIC, :west, ones(4,4,10))
    @test model.nutrients.DIC.bc.west == ones(4,4,10)

    Gcs = nutrients_init(CPU(), grid)
    apply_bcs!(Gcs, model.nutrients, model.grid, 10, 1, CPU())
    @test Gcs.DIC.data[3,3:6,3:6] == ones(4,4) ./ 8.0

end

@testset "Fields" begin
    @testset "Fields" begin
        test_fields()
    end
    @testset "Fill halos" begin
        test_fill_halos()
    end
    @testset "Boundary conditions" begin
        test_boundary_conditions()
    end
end


text:
using Weave

filenames = ["01-basics", "02-collections", "03-types" ,
              "04-project-1DCA", "05-project-errorprop"]
notesDir = "notes/"
notebookDir = "notebooks/"

function cleanTemps(filename,notesDir;exts=[".out", ".log", ".aux"])
  for ext in exts
    rm(notesDir * filename * ext)
  end
end

for filename in filenames
  # Basics
  println("Building $(filename)...")
  file = filename * ".jmd"                              # file = filename + extension

  weave(file; doctype="md2pdf",
        out_path = notesDir)                            # weaving notes
  cleanTemps(filename,notesDir)                         # removing ugly .tex perifery
  convert_doc(file, notebookDir * filename * ".ipynb")  # weaving notebooks
end


text:
function aggregate_all_experiments(data_folder::String; write_data::Bool)
    experiment_names = readdir(data_folder)
    aggregated_dataframe = init_aggregated_experiment_dataframe()
    for current_experiment in experiment_names
        aggregated_dataframe = vcat(
            aggregated_dataframe,
            aggregate_experiment(data_folder, current_experiment)
        )
    end
    if write_data
        Feather.write(
            joinpath(data_folder, "aggregated_data.feather"),
            aggregated_dataframe
        )
    end
    return aggregated_dataframe
end

function aggregate_experiment(data_folder::String, experiment_name::String; write_data::Bool=false)
    data = read_simulation_dataframe(data_folder, experiment_name)
    graphs = read_simulation_graphs(data_folder, experiment_name)
    aggregated_dataframe = build_aggregated_dataframe(data, graphs)
    if write_data
        Feather.write(
            joinpath(data_folder, experiment_name, experiment_name * ".feather"),
            aggregated_dataframe
        )
    end
    return aggregated_dataframe
end

function read_simulation_dataframe(data_folder::String, experiment_name::String)
    data = Feather.read(
        joinpath(data_folder, experiment_name, "agents", "adata.feather")
    )
    return data
end

function read_simulation_graphs(data_folder::String, experiment_name::String)
    graph_file_names = readdir(joinpath(data_folder, experiment_name, "graphs"))
    graphs = [
        LightGraphs.loadgraph(
            joinpath(data_folder, experiment_name, "graphs", file_name),
            GraphIO.GML.GMLFormat()
        ) for file_name in graph_file_names
    ]
    return graphs
end

function build_aggregated_dataframe(data::DataFrames.DataFrame, graphs::AbstractArray)
    aggregated_dataframe = init_aggregated_experiment_dataframe()
    for rep in 1:length(graphs)
        filtered_data = filter_by_rep(data, rep)
        row = build_row(filtered_data, graphs[rep])
        push!(aggregated_dataframe, row)
    end
    return aggregated_dataframe
end

function init_aggregated_experiment_dataframe()
    aggregated_experiment_dataframe = DataFrames.DataFrame(
        NetworkType = String[],
        StubbornCount = Int64[],
        Replicate = Int64[],
        SeperateComponents = Int64[],
        SizeBiggestComponent = Int64[],
        UniqueCultures = Int64[],
        FractionZeros = Float64[]
    )
    return aggregated_experiment_dataframe
end

function filter_by_rep(dataframe::DataFrames.DataFrame, replicate)
    final_tick = maximum(dataframe.TickNr)
    filtered_dataframe = dataframe |>
        @filter(_.Replicate == replicate && _.TickNr == final_tick) |>
        DataFrames.DataFrame
    return filtered_dataframe
end

function build_row(data::DataFrames.DataFrame, graph::LightGraphs.AbstractGraph)
    network_type = get_network_type(data)
    stubborn_count = data.StubbornCount[1]
    replicate = data.Replicate[1]
    connected_component_sizes = get_connected_component_sizes(graph)
    seperate_components = length(connected_component_sizes)
    size_biggest_component = maximum(connected_component_sizes)
    unique_cultures = get_unique_cultures(data)
    fraction_zeros = get_fraction_zeros(data)
    return (
        network_type,
        stubborn_count,
        replicate,
        seperate_components,
        size_biggest_component,
        unique_cultures,
        fraction_zeros
    )
end

function get_network_type(data::DataFrames.DataFrame)
    network_type = data.NetworkType[1]
    if network_type == "watts_strogatz"
        network_type = network_type * "_" * string(data.Beta[1])
    end
    return network_type
end

function get_connected_component_sizes(graph::LightGraphs.AbstractGraph)
    connected_component_sizes = [
        length(component) for component in LightGraphs.connected_components(graph)
    ]
    return connected_component_sizes
end

function get_unique_cultures(data::DataFrames.DataFrame)
    unique_cultures = data |>
        @map(_.Culture) |>
        @unique() |>
        collect |>
        length
    return unique_cultures
end

function get_fraction_zeros(data)
    culture_dimensions = length(data.Culture[1])
    agent_count = data |> nrow
    cultures_concat = data |>
        @map(_.Culture) |>
        collect |>
        join
    fraction_zeros = length(collect(eachmatch(r"0", cultures_concat))) / (culture_dimensions * agent_count)
    return fraction_zeros
end


text:
path = "docs/src/assets/"

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->100,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exSHC1.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->100,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exSHC2.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->100,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.001:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exSHC3.png")

#------------------------------------
using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->10π,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exSHC4.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->10π,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exSHC5.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->10π,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.001:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exSHC6.png")

#-----------------------------------
using ISA, Plots
𝐶₀ = AMFMtriplet(t->10,t->exp(π),1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exSHC7.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10,t->exp(π),1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exSHC8.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10,t->exp(π),1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.001:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exSHC9.png")


text:
using IRTools, Test
using IRTools: CFG, dominators, domtree

relu(x) = (y = x > 0 ? x : 0)
ir = @code_ir relu(1)

@test dominators(CFG(ir)) ==
  Dict(1 => Set([1]), 2 => Set([1, 2]), 3 => Set([1, 3]), 4 => Set([1, 4]))

@test domtree(CFG(ir)) == (1 => [2 => [], 3 => [], 4 => []])

@test domtree(CFG(ir)', entry = 4) == (4 => [1 => [], 2 => [], 3 => []])


text:
module TrackRunning

using DataFrames: readtable, isna, deleterows!
using FileIO: load
using Images: channelview, rawview, permuteddimsview
using SmoothingSplines: predict, fit, SmoothingSpline
using Reexport
@reexport using PyPlot #: xticks, yticks, imshow, scatter, ColorMap

export readttwatch, readrun
export dist_m, time_sec, speed_kmh, pace_minkm
export mapzoom, mapprojection, downloadmap, showmap, plotrun

include("run.jl")
include("io.jl")
include("metric.jl")
include("map.jl")
include("plot.jl")

end # module TrackRunning


text:
# Add fallbacks for predict_* which require mean, mode, median.

const BadModeTypes   = Union{AbstractArray{Continuous},Table(Continuous)}
const BadMeanTypes   = Union{AbstractArray{<:Finite},Table(Finite)}
const BadMedianTypes = Union{AbstractArray{<:Finite},Table(Finite)}

# mode:
predict_mode(m::Probabilistic, fitresult, Xnew) =
    predict_mode(m, fitresult, Xnew, Val(target_scitype(m)))
predict_mode(m, fitresult, Xnew, ::Any) =
    mode.(predict(m, fitresult, Xnew))
predict_mode(m, fitresult, Xnew, ::Val{<:BadModeTypes}) =
    throw(ArgumentError("Attempting to compute mode of predictions made "*
                        "by a model expecting `Continuous` targets. "))

# mean:
predict_mean(m::Probabilistic, fitresult, Xnew) =
    predict_mean(m, fitresult, Xnew, Val(target_scitype(m)))
predict_mean(m, fitresult, Xnew, ::Any) =
    mean.(predict(m, fitresult, Xnew))
predict_mean(m, fitresult, Xnew, ::Val{<:BadMeanTypes}) =
    throw(ArgumentError("Attempting to compute mean of predictions made "*
                        "by a model expecting `Finite` targets. "))

# median:
predict_median(m::Probabilistic, fitresult, Xnew) =
    predict_median(m, fitresult, Xnew, Val(target_scitype(m)))
predict_median(m, fitresult, Xnew, ::Any) =
    median.(predict(m, fitresult, Xnew))
predict_median(m, fitresult, Xnew, ::Val{<:BadMedianTypes}) =
    throw(ArgumentError("Attempting to compute median of predictions made "*
                        "by a model expecting `Finite` targets. "))

# not in MLJModelInterface as methodswith requires InteractiveUtils
MLJModelInterface.implemented_methods(::FI, M::Type{<:MLJType}) =
    getfield.(methodswith(M), :name) |> unique

# serialization fallbacks:
MLJModelInterface.save(filename, model, fitresult; kwargs...) = fitresult
MLJModelInterface.restore(filename, model, serializable_fitresult) =
                          serializable_fitresult

# to suppress inclusion of abstract types in the model registry.
for T in (:Supervised, :Unsupervised,
          :Interval, :Static, :Deterministic, :Probabilistic)
    ex = quote
        MLJModelInterface.is_wrapper(::Type{$T}) = true
    end
    eval(ex)
end


text:
using FIX
using GDAXClient
using DataStructures

import FIX: onFIXMessage

#keep two quotes with fixed distance from average price.
#move them if price goes away further than `price_tolerance`

gadx_keys = JSON.parsefile(joinpath(homedir(),".gdax_keys"))
api_key = gadx_keys["Key"]
api_secret = gadx_keys["Secret"]
passphrase = gadx_keys["Passphrase"]

user = GDAXUser("https://api.gdax.com", "wss://ws-feed.gdax.com", api_key, api_secret, passphrase)
subscription = Dict("type" => "subscribe",
                    "product_ids" => ["BTC-EUR"],
                    "channels" => ["heartbeat", "level2", "full"])

bid_quote = FIXQuote("buy",  "BTC-EUR", 0.0001001, NaN, price_tolerance = 0.1)
ask_quote = FIXQuote("sell", "BTC-EUR", 0.0001001, NaN, price_tolerance = 0.1)

struct SpreadQuotes <: AbstractGDAXMessageHandler
    bid::FIXQuote
    ask::FIXQuote
end

function onFIXMessage(this::SpreadQuotes, msg::OrderedDict{Int64, String})
    onFIXMessage(this.bid, msg)
    onFIXMessage(this.ask, msg)
end

quoter = SpreadQuotes(bid_quote, ask_quote)

condition = Condition()
books     = GDAXOrderBooks(subscription["product_ids"], condition)

client = fixconnect(user, quoter)

#this is ugly
quoter.bid.client = client
quoter.ask.client = client

wsclient = wsconnect(user, subscription, books)

spread = NaN
exec_task  = @async begin
    while true
        wait(condition)
        avg_bid = getBidAveragePrice(books["BTC-EUR"], 0.5)
        avg_ask = getAskAveragePrice(books["BTC-EUR"], 0.5)

        mid_price = (avg_bid + avg_ask) / 2.0
        quoter.bid.price = mid_price - spread
        quoter.ask.price = mid_price + spread

        GDAXClient.maintain(quoter.bid)
        GDAXClient.maintain(quoter.ask)
    end
    println("exec task exiting")
end

#now adjust spread and tolerances as you wish
spread = 100.0
quoter.bid.price_tolerance = 10.0
quoter.ask.price_tolerance = 10.0

max_stars_repo_name: DatName/gdax-julia

max_issues_repo_name: DatName/gdax-julia

max_forks_repo_name: DatName/gdax-julia


text:
typealias OptionalFunction Union{Base.Callable,Void}

immutable FunctionMap{T}<:AbstractLinearMap{T}
    f::Function
    M::Int
    N::Int
    _ismutating::Bool
    _issym::Bool
    _ishermitian::Bool
    _isposdef::Bool
    _fT::OptionalFunction
    _fC::OptionalFunction
    function FunctionMap(f::Function,M::Int,N::Int=M;ismutating::Bool=false,issym::Bool=false,ishermitian::Bool=(T<:Real && issym),isposdef::Bool=false,ftranspose::OptionalFunction=nothing,fctranspose::OptionalFunction=nothing)
        # sanity checks
        (issym || ishermitian || isposdef) && (M!=N) && error("a symmetric or hermitian map should be square")
        T<:Real && (issym!=ishermitian) && error("a real symmetric map is also hermitian")
        isposdef && !ishermitian && error("a positive definite map should be hermitian")
        # construct
        new(f,M,N,ismutating,issym,ishermitian,isposdef,ftranspose,fctranspose)
    end
end
# additional constructor
function FunctionMap(f::Function,M::Int,N::Int=M;ismutating::Bool=false,isreal::Bool=true,issym::Bool=false,ishermitian::Bool=(isreal && issym),isposdef::Bool=false,ftranspose::OptionalFunction=nothing,fctranspose::OptionalFunction=nothing)
    T=(isreal ? Float64 : Complex128) # default assumption
    FunctionMap{T}(f,M,N;ismutating=ismutating,issym=issym,ishermitian=ishermitian,isposdef=isposdef,ftranspose=ftranspose,fctranspose=fctranspose)
end

# show
function Base.show{T}(io::IO,A::FunctionMap{T})
    print(io,"FunctionMap{$T}($(A.f),$(A.M),$(A.N);ismutating=$(A._ismutating),issym=$(A._issym),ishermitian=$(A._ishermitian),isposdef=$(A._isposdef),transpose=$(A._fT),ctranspose=$(A._fC))")
end

# properties
Base.size(A::FunctionMap,n)=(n==1 ? A.M : (n==2 ? A.N : error("AbstractLinearMap objects have only 2 dimensions")))
Base.size(A::FunctionMap)=(A.M,A.N)
Base.issym(A::FunctionMap)=A._issym
Base.ishermitian(A::FunctionMap)=A._ishermitian
Base.isposdef(A::FunctionMap)=A._isposdef

# multiplication with vector
Base.A_mul_B!(y::AbstractVector,A::FunctionMap,x::AbstractVector)=begin
    (length(x)==A.N && length(y)==A.M) || throw(DimensionMismatch())
    A._ismutating ? A.f(y,x) : copy!(y,A.f(x))
    return y
end
*(A::FunctionMap,x::AbstractVector)=begin
    length(x)==A.N || throw(DimensionMismatch())
    A._ismutating ? A.f(similar(x,promote_type(eltype(A),eltype(x)),A.M),x) : A.f(x)
end

Base.At_mul_B!(y::AbstractVector,A::FunctionMap,x::AbstractVector)=begin
    A._issym && return Base.A_mul_B!(y,A,x)
    (length(x)==A.M && length(y)==A.N) || throw(DimensionMismatch())
    if A._fT!=nothing
        return (A._ismutating ? A._fT(y,x) : copy!(y,A._fT(x)))
    elseif A._fC!=nothing
        return (A._ismutating ? conj!(A._fC(y,conj(x))) : copy!(y,conj(A._fC(conj(x)))))
    else
        error("transpose not implemented for $A")
    end
end
Base.At_mul_B(A::FunctionMap,x::AbstractVector)=begin
    A._issym && return A*x
    length(x)==A.M || throw(DimensionMismatch())
    if A._fT!=nothing
        return (A._ismutating ? A._fT(similar(x,promote_type(eltype(A),eltype(x)),A.N),x) : A._fT(x))
    elseif A._fC!=nothing
        return (A._ismutating ? conj!(A._fC(similar(x,promote_type(eltype(A),eltype(x)),A.N),conj(x))) : conj(A._fC(conj(x))))
    else
        error("transpose not implemented for $A")
    end
end
Base.Ac_mul_B!(y::AbstractVector,A::FunctionMap,x::AbstractVector)=begin
    A._ishermitian && return Base.A_mul_B!(y,A,x)
    (length(x)==A.M && length(y)==A.N) || throw(DimensionMismatch())
    if A._fC!=nothing
        return (A._ismutating ? A._fC(y,x) : copy!(y,A._fC(x)))
    elseif A._fT!=nothing
        return (A._ismutating ? conj!(A._fT(y,conj(x))) : copy!(y,conj(A._fT(conj(x)))))
    else
        error("ctranspose not implemented for $A")
    end
end
Base.Ac_mul_B(A::FunctionMap,x::AbstractVector)=begin
    A._ishermitian && return A*x
    length(x)==A.M || throw(DimensionMismatch())
    if A._fC!=nothing
        return (A._ismutating ? A._fC(similar(x,promote_type(eltype(A),eltype(x)),A.N),x) : A._fC(x))
    elseif A._fT!=nothing
        return (A._ismutating ? conj!(A._fT(similar(x,promote_type(eltype(A),eltype(x)),A.N),conj(x))) : conj(A._fT(conj(x))))
    else
        error("ctranspose not implemented for $A")
    end
end


text:
# Autogenerated wrapper script for MUMPS_seq_jll for aarch64-linux-gnu-libgfortran5
export libcmumps, libdmumps, libsmumps, libzmumps

using CompilerSupportLibraries_jll
using METIS_jll
using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("MUMPS_seq")
JLLWrappers.@declare_library_product(libcmumps, "libcmumps.so")
JLLWrappers.@declare_library_product(libdmumps, "libdmumps.so")
JLLWrappers.@declare_library_product(libsmumps, "libsmumps.so")
JLLWrappers.@declare_library_product(libzmumps, "libzmumps.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, METIS_jll, OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libcmumps,
        "lib/libcmumps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdmumps,
        "lib/libdmumps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsmumps,
        "lib/libsmumps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libzmumps,
        "lib/libzmumps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## Example for using DSANet

# Make sure all the required packages are available
cd(@__DIR__)
using Pkg; Pkg.activate(".")
Pkg.instantiate()

@info "Loading packages"
using FluxArchitectures
using Random
using Plots

# Load some sample data
@info "Loading data"
poollength = 10
horizon = 15
datalength = 4000
input, target = get_data(:exchange_rate, poollength, datalength, horizon) |> gpu

# Define the network architecture
inputsize = size(input, 1)
local_length = 3
n_kernels = 3
d_model = 4
hiddensize = 1
n_layers = 3
n_head = 2

# Define the neural net
@info "Creating model and loss"
Random.seed!(123)
model = DSANet(inputsize, poollength, local_length, n_kernels, d_model,
               hiddensize, n_layers, n_head) |> gpu

# MSE loss
function loss(x, y)
    Flux.reset!(model)
    return Flux.mse(model(x), permutedims(y))
end

# Callback for plotting the training
cb = function ()
    Flux.reset!(model)
    pred = model(input) |> permutedims |> cpu
    Flux.reset!(model)
    p1 = plot(pred, label="Predict")
    p1 = plot!(cpu(target), label="Data", title="Loss $(loss(input, target))")
    display(plot(p1))
end

# Training loop
@info "Starting training"
@info "Start loss" loss = loss(input, target)
Flux.train!(loss, Flux.params(model),Iterators.repeated((input, target), 50),
             ADAM(0.005), cb=cb)

@info "Finished"
@info "Final loss" loss = loss(input, target)


text:
function generate_toy_data(rng::AbstractRNG)
    x = range(0, 10; length=5)
    ys = rand(rng, GP(SEKernel())(x, 1e-6), 3)
    y1 = ys[:, 1]
    y2 = ys[:, 2]
    y3 = ys[:, 3]
    indices = randcycle(rng, 5)
    x_train = zeros(3)
    y_1_train = zeros(3)
    y_2_train = zeros(3)
    y_3_train = zeros(3)
    x_test = zeros(2)
    y_1_test = zeros(2)
    y_2_test = zeros(2)
    y_3_test = zeros(2)
    for (i, val) in enumerate(indices)
        if i <= 3
            x_train[i] = x[val]
            y_1_train[i] = y1[val]
            y_2_train[i] = y2[val]
            y_3_train[i] = y3[val]
        else
            x_test[i - 3] = x[val]
            y_1_test[i - 3] = y1[val]
            y_2_test[i - 3] = y2[val]
            y_3_test[i - 3] = y3[val]
        end
    end
    x_train = MOInputIsotopicByOutputs(x_train, 3)
    x_test = MOInputIsotopicByOutputs(x_test, 3)
    y_train = vcat(y_1_train, y_2_train, y_3_train)
    y_test = vcat(y_1_test, y_2_test, y_3_test)

    return x_train, x_test, y_train, y_test
end

# Posterior should approximately agree with observations samples from the prior, if
# observed under a sufficiently small amount of noise.
# Testing for really eggregious bugs, as opposed to numerical issues, so tolerances are
# quite loose.
function test_sampling_consistency(rng, f, x; rtol=1e-2, atol=1e2, σ²=1e-6)
    fx = f(x, σ²)
    y = rand(fx)
    f_post = posterior(fx, y)
    @test rand(rng, f_post(x, σ²)) ≈ y rtol = rtol
    @test mean(f_post(x)) ≈ y rtol = rtol
    @test var(f_post(x)) ≈ zeros(length(y)) rtol = rtol atol = atol
end

function approx_equivalent(rng::AbstractRNG, fx1::FiniteGP, fx2::FiniteGP)
    @test length(fx1) == length(fx2)

    @test mean(fx1) ≈ mean(fx2)
    @test var(fx1) ≈ var(fx2)
    @test cov(fx1) ≈ cov(fx2)

    y = rand(rng, fx1)
    @test logpdf(fx1, y) ≈ logpdf(fx2, y)
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	ucommit(EP::Model, inputs::Dict, UCommit::Int)

This function creates decision variables and cost expressions associated with thermal plant unit commitment or start-up and shut-down decisions (cycling on/off)

**Unit commitment decision variables:**

This function defines the following decision variables:

$\nu_{y,t,z}$ designates the commitment state of generator cluster $y$ in zone $z$ at time $t$;
$\chi_{y,t,z}$ represents number of startup decisions in cluster $y$ in zone $z$ at time $t$;
$\zeta_{y,t,z}$ represents number of shutdown decisions in cluster $y$ in zone $z$ at time $t$.

**Cost expressions:**

The total cost of start-ups across all generators subject to unit commitment ($y \in UC$) and all time periods, t is expressed as:
```math
\begin{aligned}
	C^{start} = \sum_{y \in UC, t \in T} \omega_t \times start\_cost_{y} \times \chi_{y,t}
\end{aligned}
```

The sum of start-up costs is added to the objective function.
"""
function ucommit(EP::Model, inputs::Dict, UCommit::Int)

	println("Unit Commitment Module")

	dfGen = inputs["dfGen"]

	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
	COMMIT = inputs["COMMIT"] # For not, thermal resources are the only ones eligible for Unit Committment

	### Variables ###

	## Decision variables for unit commitment
	# commitment state variable
	@variable(EP, vCOMMIT[y in COMMIT, t=1:T] >= 0)
	# startup event variable
	@variable(EP, vSTART[y in COMMIT, t=1:T] >= 0)
	# shutdown event variable
	@variable(EP, vSHUT[y in COMMIT, t=1:T] >= 0)

	### Expressions ###

	## Objective Function Expressions ##

	# Startup costs of "generation" for resource "y" during hour "t"
	@expression(EP, eCStart[y in COMMIT, t=1:T],(inputs["omega"][t]*inputs["C_Start"][y]*vSTART[y,t]))

	# Julia is fastest when summing over one row one column at a time
	@expression(EP, eTotalCStartT[t=1:T], sum(eCStart[y,t] for y in COMMIT))
	@expression(EP, eTotalCStart, sum(eTotalCStartT[t] for t=1:T))

	EP[:eObj] += eTotalCStart

	### Constratints ###
	## Declaration of integer/binary variables
	if UCommit == 1 # Integer UC constraints
		for y in COMMIT
			set_integer.(vCOMMIT[y,:])
			set_integer.(vSTART[y,:])
			set_integer.(vSHUT[y,:])
			if y in inputs["RET_CAP"]
				set_integer(EP[:vRETCAP][y])
			end
			if y in inputs["NEW_CAP"]
				set_integer(EP[:vCAP][y])
			end
		end
	end

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
#=

# In-depth tutorial

This objective of this in-depth tutorial is to introduce the reader to the functionality of this library and reach a level of familiarity sufficient for building more complex applications.

This tutorial is *not* intended as a replacement for a Vulkan tutorial. In particular, it is assumed that the reader has a basic understanding of Vulkan. If you are new to Vulkan, feel free to follow the [official Vulkan tutorial](https://vulkan-tutorial.com/) along with this one. The Vulkan tutorial will teach you the concepts behind Vulkan, and this tutorial how to use the API it *from Julia*.
A lot of resources are available online for learning about Vulkan, such as the [Vulkan Guide](https://github.com/KhronosGroup/Vulkan-Guide) by the Khronos Group. You can find a more detailed list [here](https://www.vulkan.org/learn).

## Initialization

The entry point of any Vulkan application is an `Instance`, so we will create one. We will use validation layers and the extension `VK_EXT_debug_utils` that will allow logging from Vulkan. We will also provide an `ApplicationInfo` parameter to request the 1.2 version of the API.

=#

using SwiftShader_jll # hide
using Vulkan
@set_driver :SwiftShader # hide

const application_info = ApplicationInfo(
    v"0.0.1", # application version
    v"0.0.1", # engine version
    v"1.2"; # requested API version
    application_name="Demo",
    engine_name="DemoEngine",
)

#=

The application and engine versions don't matter here, but must be provided regardless. We request the version 1.2 of the Vulkan API (ensure you have a driver compatible with the 1.2 version first), which is the latest minor version. Application and engine names won't matter either, but we provide them for demonstration purposes.

=#

const instance = Instance(
    ["VK_LAYER_KHRONOS_validation"],
    ["VK_EXT_debug_utils"];
    application_info,
)

#=

This simple call does a few things under the hood:
- it creates an `InstanceCreateInfo` with the provided arguments
- it calls `create_instance` with the create info, which in turn:
  - calls the API constructor `vkCreateInstance`
  - checks if an error occured; if so, return a `ResultTypes.Result` type wrapping an error
  - registers a finalizer to a newly created `Instance` that will call `destroy_instance` (forwarding to `vkDestroyInstance`) when necessary
- unwraps the result of `create_instance`, which is assumed to be a success code (otherwise an exception is thrown).

Note that this little abstraction does not induce any loss of functionality. Indeed, the `Instance` constructor has a few keyword arguments not mentioned above for a more advanced use, which simply provides default values.

Note that we pass in arrays, version numbers and strings; but the C API does not know anything about Julia types. Fortunately, these conversions are taken care of by the wrapper, so that we don't need to provide pointers for arrays and strings nor integers that act as version numbers.

We now setup a debug messenger that we'll use for logging. Its function is to process messages sent by the Vulkan API. We could use the default debug callback provided by Vulkan.jl, namely [`default_debug_callback`](@ref); but instead we will implement our own callback for educational purposes. We'll just define a function that prints whatever message is received from Vulkan.

We won't just `println`, because it does context-switching which is not allowed in finalizers (and the callback may be called in a finalizer, notably when functions like `vkDestroy...` are called). We can use `jl_safe_printf` which does not go through the Julia task system to safely print messages. The data that will arrive from Vulkan will be a `Ptr{core.VkDebugUtilsMessengerCallbackDataEXT}`

=#

function debug_callback(severity, type, p_data, p_user_data)
    p_data ≠ C_NULL || return UInt32(0) # don't print if there's no message
    data = unsafe_load(p_data)
    msg = unsafe_string(data.pMessage)
    ccall(:jl_safe_printf, Cstring, (Cstring,), string(msg, '\n'))
    return UInt32(0)
end

# Because we are passing a callback to Vulkan as a function pointer, we need to convert it to a function pointer using `@cfunction`:

const debug_callback_c = @cfunction(
    default_debug_callback,
    UInt32,
    (
        DebugUtilsMessageSeverityFlagEXT,
        DebugUtilsMessageTypeFlagEXT,
        Ptr{core.VkDebugUtilsMessengerCallbackDataEXT},
        Ptr{Cvoid},
    )
)

#=

!!! warning
    If you intend to do this inside a module that will be precompiled, you must create the function pointer in the `__init__()` stage:
    ```julia
    const debug_callback_c = Ref{Ptr{Cvoid}}()
    function __init__()
        debug_callback_c[] = @cfunction(...) # the expression above
    end
    ```
    This is because function pointers are only valid in a given runtime environment, so you need to get it at runtime (and *not* compile time).

Note that the signature uses pointers and structures from VulkanCore.jl (accessible through the exported `core` and `vk` aliases). This is because we currently don't generate wrappers for defining function pointers. The flag types are still wrapper types, because the wrapped versions share the same binary representation as the core types. Let's create the debug messenger for all message types and severities:

=#

const debug_messenger = DebugUtilsMessengerEXT(
    instance,
    |(
        DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
        DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT,
        DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT,
        DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT,
    ),
    |(
        DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT,
        DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT,
        DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT,
    ),
    debug_callback_c,
)

#=

!!! note
    `DebugUtilsMessengerEXT` is an extension-defined handle. Any extension function such as `vkCreateDebugUtilsMessengerEXT` and `vkDestroyDebugUtilsMessengerEXT` (called in the constructor and finalizer respectively) must be called using a function pointer. This detail is abstracted away with the wrapper, as API function pointers are automatically retrieved as needed and stored in a thread-safe global dispatch table. See more in the [Dispatch](@ref) section.

We can now enumerate and pick a physical device that we will use for this tutorial.

=#


text:
Base.broadcastable(win::AbstractWindow) = Ref(win)
Base.broadcastable(wm::AbstractWindowManager) = Ref(wm)

"""
Window dimensions (width, height).
"""
extent(win::AbstractWindow) = not_implemented_for(win)

"""
Get window DPI.
"""
dpi(win::AbstractWindow) = not_implemented_for(win)

"""
Resize the window with `extent` reflecting the new (width, height).
"""
set_extent(win::AbstractWindow, extent) = not_implemented_for(win)

"""
Clean up window state upon destruction. The window cannot be assumed to be mapped or to be valid.
"""
terminate_window!(wm::AbstractWindowManager, win::AbstractWindow) = not_implemented_for(wm)

"""
Map a window to the screen.
"""
map_window(win::AbstractWindow) = not_implemented_for(win)

"""
Unmap a window from the screen.
"""
unmap_window(win::AbstractWindow) = not_implemented_for(win)

"""
Set window title.
"""
set_title(win::AbstractWindow, title) = not_implemented_for(win)

"""
Set window icon.
"""
set_icon(win::AbstractWindow, icon) = not_implemented_for(win)

"""
Set window icon title.
"""
set_icon_title(win::AbstractWindow, icon_title) = not_implemented_for(win)

attach_graphics_context!(win::AbstractWindow, ctx) = not_implemented_for(win)

get_window_symbol(wm::AbstractWindowManager, win) = not_implemented_for(wm)

get_window(wm::AbstractWindowManager, id) = not_implemented_for(wm)

callbacks(wm::AbstractWindowManager, id) = not_implemented_for(wm)


text:
function is_raw_ssh_private_key(content::AbstractString)
    x1 = occursin("-", content)
    x2 = occursin(" ", content)
    x3 = occursin("BEGIN ", content)
    x4 = occursin("END ", content)
    x5 = occursin(" PRIVATE KEY", content)

    return x1 && x2 && x3 && x4 && x5
end

function decode_ssh_private_key(content::AbstractString)
    if is_raw_ssh_private_key(content)
        @info("This is a raw SSH private key.")
        return content
    end

    @info(
        "This doesn't look like a raw SSH private key. I will assume that it is a Base64-encoded SSH private key. I will now try to Base64-decode it."
    )
    decoded_content = String(Base64.base64decode(content))

    if !is_raw_ssh_private_key(decoded_content)
        msg = string(
            "Private key is not in the correct format. ",
            "It must be a private key in PEM format. ",
            "You should use the OpenSSH `ssh-keygen` command to generate the private key. ",
            "You should pass the `-m PEM` option to ensure that the private key is in PEM format.",
        )
        throw(UnableToParseSSHKey(msg))
    end

    @info(
        "This was a Base64-encoded SSH private key. I Base64-decoded it, and now I have the raw SSH private key."
    )
    return decoded_content
end


text:
module _Simple

using Brainstorm._Algorithm: binarysearch

include("binarysplitting.jl")

# Compute int(e * 10^digits)
# This is done using Taylor series with binary splitting
euler(digits::Integer) = begin
    # How many terms to compute
    #
    # http://www.johndcook.com/blog/2011/06/10/stirling-approximation/
    #
    #             k! > 10^digits
    # =>      log k! > digits
    # => k log k - k > digits
    #
    # with log n! = sum_k=1^n log k
    #            >= int_x=1^n log x dx
    #             = n log n - n + 1
    #             > n log n - n
    #
    n = binarysearch(termfunc, digits + 1, one(digits), digits + 100)[2]

    # calculate P(0, N) and Q(0, N)
    p, q = binarysplitting(zero(n), n)

    o = big(10)^digits # one
    o + (o * p) ÷ q
end

termfunc(k) = k * (ndigits(k) - 2)

compute(a, b) = big(1), big(b)

combine(am, mb) = begin
    p₁, q₁ = am
    p₂, q₂ = mb

    p = p₁ * q₂ + p₂
    q = q₁ * q₂

    p, q
end

end


text:
using RecipesBase
using StatsBase: autocor

struct ConvErrorBars
    convergence::ConvergenceHistory
    exponents::AbstractVector
    i::Int
end

@recipe function f(data::ConvErrorBars;
                   threshold_color = :green,
                   unstable_color = :darkblue,
                   stable_color = :darkred)
    history = data.convergence
    exponents = data.exponents
    i = data.i

    if isempty(history.kinds)
        label := ""
        return []  # to suppress the warning about makershape
    end

    kind_to_color = Dict(
        UnstableConvError => unstable_color,
        StableConvError => stable_color,
    )
    errorcolor = [kind_to_color[k] for k in history.kinds]
    thresholdcolor = [threshold_color for _ in errorcolor]
    invisibles = [nothing for _ in errorcolor]

    history.details
    errormarker = [
        if (k == UnstableConvError) && ! history.details[i][j].tail_ok
            :xcross
        else
            :hline
        end
        for (j, k) in enumerate(history.kinds)
    ]
    thresholdmarker = [:hline for _ in 1:length(errormarker)]

    seriestype := :scatter
    markersize --> [5 2]
    label := ""
    markershape --> [errormarker thresholdmarker]
    # markerstrokestyle := [:solid :dash]
    markercolor := [errorcolor invisibles]
    markerstrokecolor := [errorcolor thresholdcolor]
    yerror := [history.errors[i] history.thresholds[i]]
    x = history.orth
    l = [exponents[n] for n in history.orth]
    y = [l l]

    (x, y)
    # nothing
end

""" Plot `LERecordingSolver` via `RecipesBase`."""
@recipe function f(solver::Union{LESolverRecFTLE,
                                 AbstractRenormalizer{<: LESolRecFTLE},
                                 LESolRecFTLE};
                   correlation = true,
                   known_exponents=nothing)
    le_hist, ci_hist = exponents_stat_history(solver)
    dim_lyap = size(le_hist)[1]
    known_exponents = known_exponents == nothing ? [] : known_exponents

    if correlation
        layout := (dim_lyap + 1, 1)
    else
        layout := (dim_lyap, 1)
    end

    ylims = @views [[minimum(le_hist[i, :]),
                     maximum(le_hist[i, :])]
                    for i = 1:dim_lyap]
    for i in 1:min(dim_lyap, length(known_exponents))
        ylims[i][1] = min(ylims[i][1], known_exponents[i])
        ylims[i][2] = max(ylims[i][2], known_exponents[i])
    end
    for i in 1:dim_lyap
        ymin, ymax = ylims[i]
        dy = ymax - ymin
        ylims[i] = [ymin - dy * 0.05,
                    ymax + dy * 0.05]
    end

    for i in 1:dim_lyap
        @series begin
            subplot := i
            label --> ""
            color --> 1
            fillcolor --> 2
            fillalpha --> 0.5
            fillrange := hcat(le_hist[i, 2:end] - ci_hist[i, 2:end],
                              le_hist[i, 2:end] + ci_hist[i, 2:end])
            (2:size(le_hist, 2),
             [le_hist[i, 2:end] le_hist[i, 2:end]])
            # Note: I think it should work w/o "2:end" in principle
            # but it looks like that Plots.jl has some trouble
            # rendering fillrange with NaN.
        end
        if has_convergence_history(solver, i)
            @series begin
                subplot := i
                data = ConvErrorBars(
                    convergence_history(solver),
                    view(le_hist, i, :),
                    i)
                (data,)
            end
        end
        if i <= length(known_exponents)
            @series begin
                subplot := i
                linetype := :hline
                label --> ""
                [known_exponents[i]]
            end
        end
        @series begin
            subplot := i
            label := ""
            ylabel := "LE$i"
            if i == dim_lyap
                xlabel := "Number of orthonormalizations"
            end
            ylim --> ylims[i]
            xlim --> (1, Inf)
            xscale --> :log10
            []
        end
    end

    if correlation && length(ftle_history(solver)) > 1
        @series begin
            lags = 0:ceil(Int, sqrt(length(ftle_history(solver))))
            corrs = [autocor(ftle_history(solver, i), lags)
                     for i in 1:dim_lyap]

            subplot := dim_lyap + 1
            ylabel := "Correlation"
            xlabel := "Lag (orthonormalizations)"
            label := ["FTLE$i" for i in 1:dim_lyap]
            legend := :top

            (lags, corrs)
        end
    end
end

@recipe function f(demo::LEDemo)
    known_exponents --> demo.example.known_exponents
    demo.solver
end


text:
## --- Define function for picking random points

    function random_lat_lon(n)
        randlon = rand(n)*360-180
        randlat = 90 - acos.(rand(n)*2-1)*180/pi

        return (randlat, randlon)
    end

## --- Define functions for querying macrostrat

    function query_macrostrat(lat, lon, zoom)
        resp = HTTP.get("https://macrostrat.org/api/mobile/map_query?lat=$lat&lng=$lon&z=$zoom")
        str = String(resp.body)
        parsed = JSON.Parser.parse(str)
        try
            parsed["success"]["data"]["burwell"][1]["lith"]
        catch error
            resp = HTTP.get("https://macrostrat.org/api/mobile/map_query?lat=$lat&lng=$lon&z=1")
            str = String(resp.body)
            parsed = JSON.Parser.parse(str)
        end
        return parsed
    end

    function get_macrostrat_min_age(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["t_int_age"]::Number
        catch error
            return NaN
        end
    end

    function get_macrostrat_max_age(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["b_int_age"]::Number
        catch error
            return NaN
        end
    end

    function get_macrostrat_map_id(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["map_id"]::Number
        catch error
            return NaN
        end
    end

    function get_macrostrat_lith(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["lith"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_descrip(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["descrip"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_name(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["name"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_strat_name(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["strat_name"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_comments(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["comments"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_ref_title(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["ref"]["ref_title"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_ref_authors(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["ref"]["authors"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_ref_year(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["ref"]["ref_year"]
        catch error
            return "NA"
        end
    end

    function get_macrostrat_ref_doi(jobj)
        try
            return jobj["success"]["data"]["burwell"][1]["ref"]["isbn_doi"]
        catch error
            return "NA"
        end
    end

## --- Functions for dealing with SRTM15
    resourcepath = "data"

    # Read srtm15plus file from HDF5 storage, downloading from cloud if necessary
    function get_srtm15plus_aveslope(varname="")
        # Available variable names: "slope", "y_lat_cntr", "x_lon_cntr",
        # "nanval", "cellsize", "scalefactor", and "reference". Units are
        # meters of elevation and decimal degrees of latitude and longitude

        # Construct file path
        filepath = joinpath(resourcepath,"srtm15plus_aveslope.h5")

        # Download HDF5 file from Google Cloud if necessary
        if ~isfile(filepath)
            print("Downloading srtm15plus.h5 from google cloud storage\n")
            download("https://storage.googleapis.com/statgeochem/srtm15plus_aveslope.h5", filepath)
        end

        # Read and return the file
        return h5read(filepath, "vars/"*varname)
    end

    # Find the elevation of points at position (lat,lon) on the surface of the
    # Earth, using the SRTM15plus 15-arc-second elevation model.
    function find_srtm15plus_aveslope(srtm15plus,lat,lon)

        # Interpret user input
        if length(lat) != length(lon)
            error("lat and lon must be equal length\n")
        elseif isa(srtm15plus,Dict)
            data = srtm15plus["slope"]
        elseif isa(srtm15plus, Array)
            data = srtm15plus
        else
            error("wrong srtm15plus variable")
        end

        # Scale factor (cells per degree) = 60 * 4 = 240
        # (15 arc seconds goes into 1 arc degree 240 times)
        sf = 240

        # Create and fill output vector
        out=Array{Float64}(size(lat));
        for i=1:length(lat)
            if isnan(lat[i]) || isnan(lon[i]) || lat[i]>90 || lat[i]<-90 || lon[i]>180 || lon[i]<-180
                # Result is NaN if either input is NaN or out of bounds
                out[i] = NaN
            else
                # Convert latitude and longitude into indicies of the elevation map array
                # Note that STRTM15 plus has N+1 columns where N = 360*sf
                row = 1 + round(Int,(90+lat[i])*sf)
                col = 1 + round(Int,(180+lon[i])*sf)
                # Find result by indexing
                res = data[row,col]
                if res > 1000
                    out[i] = NaN
                else
                    out[i] = res
                end
            end
        end

        return out
    end

    function find_srtm15plus_aveslope_around(srtm15plus,lat,lon; halfwidth=1::Integer, max_allowed_slope=1000::Number)

        # Interpret user input
        if length(lat) != length(lon)
            error("lat and lon must be equal length\n")
        elseif isa(srtm15plus,Dict)
            data = srtm15plus["slope"]
        elseif isa(srtm15plus, Array)
            data = srtm15plus
        else
            error("wrong srtm15plus variable")
        end

        # Scale factor (cells per degree) = 60 * 4 = 240
        # (15 arc seconds goes into 1 arc degree 240 times)
        sf = 240

        # Create and fill output vector
        out=Array{Float64}(size(lat));
        for i=1:length(lat)
            if isnan(lat[i]) || isnan(lon[i]) || lat[i]>90 || lat[i]<-90 || lon[i]>180 || lon[i]<-180
                # Result is NaN if either input is NaN or out of bounds
                out[i] = NaN
            else
                # Convert latitude and longitude into indicies of the elevation map array
                # Note that STRTM15 plus has N+1 columns where N = 360*sf
                row = 1 + round(Int,(90+lat[i])*sf)
                col = 1 + round(Int,(180+lon[i])*sf)

                # Find result by indexing
                k = 0;
                out[i] = 0;
                for r = (row-halfwidth):(row+halfwidth)
                    for c = (col-halfwidth):(col+halfwidth)
                        if r>0 && r<43202
                            res = data[r,mod(c-1,86400)+1]
                            if res < max_allowed_slope
                                k +=1
                                out[i] += res
                            end
                        end
                    end
                end
                out[i] /= k
            end
        end

        return out
    end


text:
using Documenter, ParticleScattering

makedocs(
    sitename = "ParticleScattering.jl",
    authors = "Boaz Blankrot",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = Any[
        "Home" => "index.md",
        "Tutorials" => Any[
            "tutorial1.md",
            "tutorial2.md",
            "tutorial_optim_angle.md",
            "tutorial_optim_radius.md"
        ],
        "Choosing Minimal N and P" => "minimalNP.md",
        "Incident Field Types" => "incident_fields.md",
        "Adding New Shapes" => "new_shapes.md",
        "API" => "api.md"
    ]
)

deploydocs(
    repo   = "github.com/bblankrot/ParticleScattering.jl.git",
)


text:
@testset "isolateroots" begin
    prec = 128
    RR = RealField(prec)

    problems = [(cos, 0, 10, [RR(π)/2, 3RR(π)/2, 5RR(π)/2]),
                (x -> (x - 1)*(x - 2)*(x - 3), 0, 5, [1, 2, 3]),
                (x -> exp(x) - 1, -eps(), 5, [0])]

    for (f, a, b, roots) in problems
        a = RR(a)
        b = RR(b)
        roots = RR.(roots)

        for evaltype in (:ball, :taylor)
            found, flags = isolateroots(f, a, b, evaltype = evaltype)

            @test length(found) == length(roots) && length(flags) == length(roots)
            if evaltype == :ball
                @test !any(flags)
            else
                @test all(flags)
            end
            @test all(contains(ArbTools.setinterval(found[i]...), roots[i]) for i in 1:length(roots))
        end
    end
end


text:
function newton(f::Function,x::Array{T,1},tol::T,maxiters::S) where {T<: AbstractFloat, S<: Int}

  n = length(x)
  x_new = similar(x)

  iter = 0
  len = Inf
  while len > tol

    g = derivative(f,x)
  	update = g\f(x)
  	for i = 1:n
  	  x_new[i] = x[i] - update[i]
  	end

    len = maximum(abs,x_new-x)
    x = copy(x_new)

    iter += 1
    if iter >= maxiters
      break
    end

  end

  return x_new, f(x_new), iter

end


text:
abstract type AVR <: GeneratorComponent end

@def avr_ports begin
      state_input = Vector{Symbol}()
      inner_input = [V_pss_var, VI_gen_var, VR_gen_var]
end

"""
Parameters of a AVR that returns a fixed voltage to the rotor winding:
# Conmutable structor
```julia
AVRFixed(Emf)
```

# Arguments
* `Emf`   ::Float64 : Fixed voltage to the rotor winding

"""
mutable struct AVRFixed <: AVR
      Emf::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function AVRFixed(Emf::Float64)
            n_states = 0
            states = Vector{Symbol}()
            @avr_ports
            new(Emf, n_states, states, Ports(state_input, inner_input))
      end
end

"""
Parameters of a simple proportional AVR in the derivative of EMF
i.e. an integrator controller on EMF

# Conmutable structor
```julia
AVRSimple(Kv)
```

#Arguments
* `Kv`   ::Float64 : Integral Gain

"""
mutable struct AVRSimple <: AVR
     Kv::Float64
     n_states::Int64
     states::Vector{Symbol}
     ports::Ports

     function AVRSimple(Kv::Float64)
        n_states = 1
        states = [:Vf]
        @avr_ports
        new(Kv, n_states, states, Ports(state_input, inner_input))
     end

end

"""
Parameters of an Automatic Voltage Regulator Type I - Resembles IEEE Type DC1

# Conmutable structor
```julia
AVRTypeI(Ka, Ke, Kf, Ta, Tf, Te, Tr, Vr_max, Vr_min, Ae, Be)
```

# Arguments
* `Ka`::Float64 : Amplifier Gain
* `Ke`::Float64 : Field circuit integral deviation
* `Kf`::Float64 : Stabilizer Gain in s * pu/pu
* `Ta`::Float64 : Amplifier Time Constant in s
* `Te`::Float64 : Field Circuit Time Constant in s
* `Tf`::Float64 : Stabilizer Time Constant in s
* `Tr`::Float64 : Voltage Measurement Time Constant in s
* `Vr_max`::Float64 : Maximum regulator voltage in pu
* `Vr_min`::Float64 : Minimum regulator voltage in pu
* `Ae`::Float64 : 1st ceiling coefficient
* `Be`::Float64 : 2nd ceiling coefficient
"""
mutable struct AVRTypeI <: AVR
      Ka::Float64
      Ke::Float64
      Kf::Float64
      Ta::Float64
      Te::Float64
      Tf::Float64
      Tr::Float64
      Vr_max::Float64
      Vr_min::Float64
      Ae::Float64
      Be::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function AVRTypeI(Ka::Float64,
                        Ke::Float64,
                        Kf::Float64,
                        Ta::Float64,
                        Te::Float64,
                        Tf::Float64,
                        Tr::Float64,
                        Vr_max::Float64,
                        Vr_min::Float64,
                        Ae::Float64,
                        Be::Float64) where {P<:PSS}

         n_states = 4
         states = [:Vf, :Vr1, :Vr2, :Vm]
         @avr_ports
         new(Ka, Ke, Kf, Ta, Te, Tf, Tr, Vr_max,
             Vr_min, Ae, Be, n_states, states, Ports(state_input, inner_input))
      end
end

"""
Parameters of an Automatic Voltage Regulator Type II -  Typical static exciter model

# Conmutable structor
```julia
AVRTypeII(K0, T1, T2, T3, T4, Te, Tr, Vr_max, Vr_min, Ae, Be)
```

# Arguments
* `K0`::Float64 : Regulator Gain
* `T1`::Float64 : First Pole in s
* `T2`::Float64 : First zero in s
* `T3`::Float64 : First Pole in s
* `T4`::Float64 : First zero in s
* `Te`::Float64 : Field Circuit Time Constant in s
* `Tr`::Float64 : Voltage Measurement Time Constant in s
* `Vr_max`::Float64 : Maximum regulator voltage in pu
* `Vr_min`::Float64 : Minimum regulator voltage in pu
* `Ae`::Float64 : 1st ceiling coefficient
* `Be`::Float64 : 2nd ceiling coefficient
"""
mutable struct AVRTypeII <: AVR
      K0::Float64
      T1::Float64
      T2::Float64
      T3::Float64
      T4::Float64
      Te::Float64
      Tr::Float64
      Vr_max::Float64
      Vr_min::Float64
      Ae::Float64
      Be::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function AVRTypeII(K0::Float64,
                        T1::Float64,
                        T2::Float64,
                        T3::Float64,
                        T4::Float64,
                        Te::Float64,
                        Tr::Float64,
                        Vr_max::Float64,
                        Vr_min::Float64,
                        Ae::Float64,
                        Be::Float64) where {P<:PSS}

         n_states = 4
         states = [:Vf, :Vr1, :Vr2, :Vm]
         @avr_ports
         new(K0, T1, T2, T3, T4, Te, Tr, Vr_max,
             Vr_min, Ae, Be, n_states, states, Ports(state_input, inner_input))
      end
end

"""
Parameters of an Automatic Voltage Regulator Type II from PSAT Manual -  Typical static exciter model

# Conmutable structor
```julia
AVRTypeII(K0, T1, T2, T3, T4, Te, Tr, Vr_max, Vr_min, Ae, Be)
```

# Arguments
* `Kf`::Float64 : Regulator Gain
* `Ta`::Float64 : Amplifier time constant
* `Kf`::Float64 : Stabilizer gain
* `Tf`::Float64 : Stabilizer time constant
* `Te`::Float64 : Field Circuit Time Constant in s
* `Tr`::Float64 : Voltage Measurement Time Constant in s
* `Vr_max`::Float64 : Maximum regulator voltage in pu
* `Vr_min`::Float64 : Minimum regulator voltage in pu
* `Ae`::Float64 : 1st ceiling coefficient
* `Be`::Float64 : 2nd ceiling coefficient
"""
mutable struct AVRTypeIIManual <: AVR
      Ka::Float64
      Ta::Float64
      Kf::Float64
      Tf::Float64
      Te::Float64
      Tr::Float64
      Vr_max::Float64
      Vr_min::Float64
      Ae::Float64
      Be::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function AVRTypeIIManual(Ka::Float64,
                        Ta::Float64,
                        Kf::Float64,
                        Tf::Float64,
                        Te::Float64,
                        Tr::Float64,
                        Vr_max::Float64,
                        Vr_min::Float64,
                        Ae::Float64,
                        Be::Float64) where {P<:PSS}

         n_states = 4
         states = [:Vf, :Vr1, :Vr2, :Vm]
         @avr_ports
         new(Ka, Ta, Kf, Tf, Te, Tr, Vr_max,
             Vr_min, Ae, Be, n_states, states, Ports(state_input, inner_input))
      end
end


text:
using OrdinaryDiffEq, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: van

prob1 = ODEProblem(van,[0,2.],(0.0,6),inv(0.003))
# out-of-place test
function _van(u, p, t)
  μ = p[1]
  [μ*((1-u[2]^2)*u[1] - u[2]),
   1*u[1]]
end
prob2 = ODEProblem(_van,[0,2.],(0.0,6),inv(0.003))
probArr = [prob1, prob2]
# Test if switching back and forth
is_switching_fb(sol) = maximum(diff(findall(x->x==2, sol.alg_choice))) > 5
for prob in probArr
  alg = AutoTsit5(Rodas5(); maxstiffstep=5, maxnonstiffstep=5, stiffalgfirst=true)
  sol = solve(prob, alg)
  @test length(sol.t) < 280
  @test typeof(alg.algs[sol.alg_choice[1]]) <: Rodas5
  @test is_switching_fb(sol)
  sol = solve(prob, AutoDP5(Rodas5(); maxstiffstep=2, maxnonstiffstep=2,
                            stifftol=11//10, nonstifftol=9//10),
                            reltol=1e-5, abstol=1e-5)
  @test length(sol.t) < 625
  @test is_switching_fb(sol)

  sol = solve(prob,AutoVern6(Kvaerno3(); maxstiffstep=4, maxnonstiffstep=4))
  @test length(sol.t) < 690
  @test is_switching_fb(sol)
  sol = solve(prob,AutoVern7(Hairer42(); maxstiffstep=4, maxnonstiffstep=4))
  @test length(sol.t) < 540
  @test is_switching_fb(sol)
  sol = solve(prob,AutoVern8(Rosenbrock23(); maxstiffstep=4, maxnonstiffstep=4))
  @test length(sol.t) < 910
  @test is_switching_fb(sol)
  sol = solve(prob,AutoVern9(KenCarp3(); maxstiffstep=4, maxnonstiffstep=4))
  @test length(sol.t) < 470
  @test is_switching_fb(sol)
end