text:
using Plots
using DifferentialEquations

function eulers_method(f::Function, α::Real, a::Real, b::Real, N::Int64)

    n1 = N + 1
    u = zeros(n1,2)

    h = (b - a) / N
    u[1,1] = a
    u[1,2] = α

    for i in 2:n1
        u[i,2] = u[i-1,2] + h * f(u[i-1,1], u[i-1,2])
        u[i,1] = a + (i - 1) * h
    end

    return u
end

function rk4_method(f::Function, α::Real, a::Real, b::Real, N::Int64)

    n1 = N + 1
    u = zeros(n1,2)

    h = (b - a) / N
    u[1,1] = a
    u[1,2] = α

    for i in 2:n1
        t = u[i-1,1]
        w = u[i-1,2]

        k1 = h * f(t, w)
        k2 = h * f(t + h / 2, w + k1 / 2)
        k3 = h * f(t + h / 2, w + k2 / 2)
        k4 = h * f(t + h, w + k3)

        u[i,2] = w + (k1 + 2 * k2 + 2 * k3 + k4) / 6
        u[i,1] = a + (i - 1) * h
    end

    return u
end

function main()

    # one_ode()
    system_ode()

end

function one_ode()

    df(t,y) = y - t^2 + 1
    f(t) = -0.5*exp(t) + t^2 + 2t + 1

    u0 = 0.5
    a = 0
    b = 2
    n = 10
    tlin = 0:0.2:2

    sol1 = eulers_method(df, u0, a, b, n)
    sol2 = rk4_method(df, u0, a, b, n)

    # DE library
    df(u,p,t) = u - t^2 + 1
    tspan = (0., 2.)
    prob = ODEProblem(df, u0, tspan)
    sol3 = solve(prob)

    plot(tlin, f.(tlin), label="exact", legend=:bottomright, dpi=150)

    plot!(sol3.t, sol3.u, markershape=:x, label="DE Lib")

    plot!(sol2[:,1], sol2[:,2], markershape=:+, label="rk4")
    plot!(sol1[:,1], sol1[:,2], markershape=:o, label="eulers")

end

function system_ode()

    # du - output solutions
    # u  - input points
    # p  - parameters
    # t  - time
    function df(du, u, p, t)

        du[1] = -4u[2] + cos(t)
        du[2] = u[1]

    end

    u0 = [0, 0]
    tspan = (0, 1.0)
    prob = ODEProblem(df, u0, tspan)
    sol = solve(prob)

    # plot(sol, label=["x" "y"], dpi=150)
    plot(sol, vars=(1,2), label="f(x,y,t)", dpi=150)
end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
# ---
# title: 552. Student Attendance Record II
# id: problem552
# author: Indigo
# date: 2021-08-18
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/student-attendance-record-ii/description/>
# hidden: true
# ---
#
# Given a positive integer **n** , return the number of all possible attendance
# records with length n, which will be regarded as rewardable. The answer may be
# very large, return it after mod 109 \+ 7.
#
# A student attendance record is a string that only contains the following three
# characters:
#
#   1. **'A'** : Absent.
#   2. **'L'** : Late.
#   3. **'P'** : Present.
#
# A record is regarded as rewardable if it doesn't contain **more than one 'A'
# (absent)** or **more than two continuous 'L' (late)**.
#
# **Example 1:**
#
#
#
#     Input: n = 2
#     Output: 8
#     Explanation:
#     There are 8 records with length 2 will be regarded as rewardable:
#     "PP" , "AP", "PA", "LP", "PL", "AL", "LA", "LL"
#     Only "AA" won't be regarded as rewardable owing to more than one absent times.
#
#
# **Note:** The value of **n** won't exceed 100,000.
#
#
## @lc code=start
using LeetCode

function check_record(n::Int)
    P = [1 1 0 1 0 0
        1 0 1 1 0 0
        1 0 0 1 0 0
        0 0 0 1 1 0
        0 0 0 1 0 1
        0 0 0 1 0 0]
    res = mat_fast_mul(P, n, Int(1e9 + 7))
    sum(@view res[1, :]) % Int(1e9 + 7)
end
## @lc code=end

ext: jl


text:
@recipe function f(
    result::AbstractRanks
)
    names = string.(result.names)
    ylabel = rankingcriterion(result)

    # prevent pyplot from converting underscores to subscripts
    if _plots_module().backend() == _plots_module().PyPlotBackend()
        names = replace.(names, "_" => "-")
        ylabel = replace(ylabel, "_" => "-")
    end

    @series begin
        seriestype --> :bar
        legend --> false
        yguide --> ylabel
        xrotation --> 0

        names, result.values
    end
end

function rankingcriterion(result::MeasurementRanks)
    if isa(result.criterion, SumOfSmallestIntervals)
        return "sum of rel. increases of smallest "*string(result.criterion.p*100)*"% intervals"

    elseif isa(result.criterion, SmallestInterval)
        return "rel. increase of smallest "*string(result.criterion.p*100)*"% interval of "*string(result.criterion.key)

    elseif isa(result.criterion, HighestDensityRegion)
        return "rel. increase of "*string(result.criterion.p*100)*"% HDR("*join(map(string, result.criterion.keys), ", ")*")"
    end

    return ""
end

function rankingcriterion(result::UncertaintyRanks)
    if isa(result.criterion, SumOfSmallestIntervals)
        return "sum of rel. decreases of smallest "*string(result.criterion.p*100)*"% intervals"

    elseif isa(result.criterion, SmallestInterval)
        return "rel. decrease of smallest "*string(result.criterion.p*100)*"% interval of "*string(result.criterion.key)

    elseif isa(result.criterion, HighestDensityRegion)
        return "rel. decrease of "*string(result.criterion.p*100)*"% HDR of "*string(result.criterion.keys)
    end

    return ""
end


text:
@inline function _iota!(array::Vector{Int64})
  for i in eachindex(array)
    @inbounds array[i] = i
  end
end

@inline function _intermediateOutputSerial!(smcio::SMCIO, p::Int64)
  if length(smcio.allZetas) > 0
    @inbounds vec = smcio.allZetas[p]
    _copyParticles!(vec, smcio.zetas)
    @inbounds smcio.allWs[p] .= smcio.ws
    @inbounds smcio.allEves[p] .= smcio.eves
    @inbounds p != smcio.n && (smcio.allAs[p] .= smcio.internal.as)
  end
end

function _Vhat1Serial(smcio::SMCIO)
  N::Int64 = smcio.N
  eves::Vector{Int64} = smcio.eves
  ws::Vector{Float64} = smcio.ws
  q::Int64 = 1 + smcio.internal.nresamples
  i::Int64 = 1
  v::Float64 = 0.0
  while i <= N
    @inbounds e = eves[i]
    tmp::Float64 = 0.0
    @inbounds while i <= N && eves[i] == e
      @inbounds tmp += ws[i]
      i += 1
    end
    tmp /= N
    v += tmp*tmp
  end
  return 1 - (1-v)*(N/(N-1))^q
end

@inline function _resampleSerial!(smcio::SMCIO, p::Int64, conditional::Bool)
  if p != smcio.n-1 && smcio.esses[p] > smcio.essThreshold
    smcio.resample[p] = false
    _iota!(smcio.internal.as)
    return
  end
  smcio.resample[p] = true
  smcio.internal.nresamples += 1
  if conditional
    smcio.internal.as[1] = 1
    sampleSerial!(smcio.internal.as, getRNG(), smcio.ws, smcio.N-1,
      smcio.internal.scratch1, smcio.internal.scratch2, 1)
  else
    sampleSerial!(smcio.internal.as, getRNG(), smcio.ws, smcio.N,
      smcio.internal.scratch1, smcio.internal.scratch2)
  end

end

# Serial implementation of SMC and cSMC
@inline function _smcSerial!(model::SMCModel, smcio::SMCIO{Particle},
  ref::Union{Nothing, Vector{Particle}} = nothing,
  refout::Union{Nothing, Vector{Particle}} = nothing) where Particle
  zetas = smcio.zetas
  zetaAncs = smcio.internal.zetaAncs
  lws = smcio.internal.lws
  ws = smcio.ws
  as = smcio.internal.as
  engine = getRNG()
  pScratch = smcio.internal.particleScratch
  smcio.internal.nresamples = 0

  logZhats = smcio.logZhats
  lZ = 0.0

  _iota!(smcio.eves)

  for p = 1:smcio.n
    if p > 1
      _copyParticles!(zetaAncs, zetas, as)
      smcio.internal.oldEves .= smcio.eves
      _setEves!(smcio.eves, smcio.internal.oldEves, as)
    end
    if ref == nothing
      _mutateParticles!(zetas, engine, p, model.M!, zetaAncs, pScratch)
    else
      _mutateParticles!(zetas, engine, p, model.M!, zetaAncs, pScratch, ref[p])
    end

    _logWeightParticles!(lws, p, model.lG, zetas, pScratch,
      p == 1 || smcio.resample[p-1], ws)

    smcio.internal.maxlw = maximum(lws)
    ws .= exp.(lws .- smcio.internal.maxlw)
    smcio.internal.mws = mean(ws)
    lZ += smcio.internal.maxlw + log(smcio.internal.mws)
    ws ./= smcio.internal.mws
    @inbounds logZhats[p] = lZ
    @inbounds smcio.esses[p] = smcio.N / mapreduce(x -> x*x, +, ws)

    smcio.Vhat1s[p] = _Vhat1Serial(smcio)

    p < smcio.n && _resampleSerial!(smcio, p, ref != nothing)

    smcio.fullOutput && _intermediateOutputSerial!(smcio, p)
  end

  refout != nothing && SequentialMonteCarlo.pickParticle!(refout, smcio)
  return
end


text:
"
Defines a SEIR epidemic model and dataset.

##Fields

- `yobs` Vector of observed recovery times for each population member.
- `popsize` Number of population members.
- `get_pDist` A function `f(parameters::Array{Float64, 1})` returning a `Distribution` for pressure thresholds
- `get_eDist` A function `f(parameters::Array{Float64, 1})` returning a `Distribution` on time spent in the exposed stage
- `get_iDist` A function `f(parameters::Array{Float64, 1})` returning a `Distribution` on time spent in the infectious stage
- `round_times` A function `f(x::Float64)` mapping a removal time to rounded `Float64` time which will be used in the distance function. This is the identity mapping by default.
"
type SEIRDetails <: ApplicationDetails
    yobs::Array{Float64, 1}
    popsize::Int
    get_pDist::Function
    get_eDist::Function
    get_iDist::Function
    round_times::Function
end

function SEIRDetails(yobs::Array{Float64, 1}, popsize::Int, get_pDist::Function, get_eDist::Function, get_iDist::Function)
    SEIRDetails(yobs, popsize, get_pDist, get_eDist, get_iDist, x->x)
end

"
Summary of a single epidemic, including latent variables.

##Fields

- `popsize` Size of the population
- `pDist` Distribution of pressure thresholds
- `eDist` Distribution of time spent in the exposed stage
- `iDist` Distribution of time spent in the infectious stage
- `compartmentTimes` Element `[i,j]` is when individual i leaves compartment j
- `u_pressures` Uniform draws associated with `pressure` (see below). n.b. This is 1 shorter than pressures since the first pressure is fixed at 0.
- `u_eDuration` Uniform draws associated with `eDuration` (see below)
- `u_iDuration` Uniform draws associated with `iDuration` (see below)
- `pressures` Infectious pressures associated with each individual
- `eDuration` Duration in the E state (or hypothetical duration if that state not reached)
- `iDuration` Duration in the I state (or hypothetical duration if that state not reached)
- `nInfectiousAtExposure` How many individuals are infectious at this individual's exposure time (or zero if this individual never exposed)
- `totalPressure' Infectious pressure at end of epidemic
- `loglikelihood` The log likelihood of this epidemic
"
type EpidemicSummary
    popsize::Int
    pDist::Distributions.ContinuousUnivariateDistribution
    eDist::Distributions.ContinuousUnivariateDistribution
    iDist::Distributions.ContinuousUnivariateDistribution
    compartmentTimes::Array{Float64, 2}
    u_pressures::Array{Float64, 1}
    u_eDuration::Array{Float64, 1}
    u_iDuration::Array{Float64, 1}
    pressures::Array{Float64, 1}
    eDuration::Array{Float64, 1}
    iDuration::Array{Float64, 1}
    nInfectiousAtExposure::Array{Int, 1}
    totalPressure::Float64
    loglikelihood::Float64
end

"
Plot details of an epidemic `e`.
For SIR models, `hideExposed=true` omits details of the exposed stage.
"
function plotEpidemic(e::EpidemicSummary; hideExposed=false)
    timeMatrix = e.compartmentTimes
    whichInfected = !isinf(timeMatrix[:,3])
    eTimes = timeMatrix[whichInfected, 1]
    iTimes = timeMatrix[whichInfected, 2]
    rTimes = timeMatrix[whichInfected, 3]
    times = sort([eTimes; iTimes; rTimes])
    nS = map(t -> (e.popsize - sum(eTimes.<=t)), times)
    nE = map(t -> (sum(eTimes.<=t) - sum(iTimes.<=t)), times)
    nI = map(t -> (sum(iTimes.<=t) - sum(rTimes.<=t)), times)
    nR = map(t -> sum(rTimes.<=t), times)
    d = vcat(DataFrame(Time = times, Count=nS, Stage="Susceptible"),
             DataFrame(Time = times, Count=nI, Stage="Infectious"),
             DataFrame(Time = times, Count=nR, Stage="Recovered"))
    if !hideExposed
        d = vcat(d, DataFrame(Time = times, Count=nI, Stage="Exposed"))
    end
    plot(d, x=:Time, y=:Count, colour=:Stage, Geom.step)
end

"
A `SEIRParticles` variable stores the particles used in an SMC algorithm for a SEIR application.

##Fields

- `nparticles` How many particles are used
- `epidemics` Array of `EpidemicSummary` objects for each particle
- `llest` Log of likelihood estimate
- `ε` The current acceptance threshold
- `iteration` How many SMC iterations have been performed
- `reachedBound` Whether the lower bound on log likelihood has been reached
- `schedule` Schedule of epsilon values used (or, for an adaptive algorithm, schedule so far)
- `maxz` Maximum slice sampling step size in previous iteration
"
type SEIRParticles <: Particles
    nparticles::Int
    epidemics::Array{EpidemicSummary, 1}
    llest::Float64
    ε::Float64
    iteration::Int
    reachedBound::Bool
    schedule::Array{Float64, 1}
    maxz::Float64
end

function Particles(a::SEIRDetails, nparticles::Int, parameters::Array{Float64, 1})
    epidemics = EpidemicSummary[simulateEpidemic(a, parameters) for i in 1:nparticles]
    SEIRParticles(nparticles, epidemics, 0.0, Inf, 0, false, Float64[], 0.5)
end

"
Return an `EpidemicSummary` variable given model details `a` and parameters `θ`.
"
function simulateEpidemic(a::SEIRDetails, θ::Array{Float64, 1})
    λscaled = θ[1]/a.popsize
    u_pressures = rand(a.popsize-1)
    u_eDuration = rand(a.popsize)
    u_iDuration = rand(a.popsize)
    pDist::Distributions.ContinuousUnivariateDistribution = a.get_pDist(θ)
    eDist::Distributions.ContinuousUnivariateDistribution = a.get_eDist(θ)
    iDist::Distributions.ContinuousUnivariateDistribution = a.get_iDist(θ)
    makeEpidemic(λscaled, u_pressures, u_eDuration, u_iDuration, pDist, eDist, iDist)
end

"
Return an `EpidemicSummary` variable corresponding to exposure parameter `λscaled`, and specified `pressureThresholds`, `eDuration` and `iDuration` arrays.
The `pDist`, `eDist` and `iDist` distributions must also be supplied.
Note it is assumed that `pressureThresholds[1]==0.0` i.e. individual 1 is exposed at time zero.
"
function makeEpidemic(λscaled::Float64, u_pressures::Array{Float64, 1}, u_eDuration::Array{Float64, 1}, u_iDuration::Array{Float64, 1}, pDist::Distributions.ContinuousUnivariateDistribution, eDist::Distributions.ContinuousUnivariateDistribution, iDist::Distributions.ContinuousUnivariateDistribution)
    ##Initialisation
    pressureThresholds = zeros(1+length(u_pressures))
    pressureThresholds[2:end] = quantile(pDist, u_pressures) ##nb pressureThresholds[1] is fixed at 0.0 as it's the initial infective
    eDuration = quantile(eDist, u_eDuration)
    iDuration = quantile(iDist, u_iDuration)
    popsize = length(pressureThresholds)
    compartmentTimes = fill(Inf, (popsize, 3))
    nInfectiousAtExposure = [1; zeros(Int, popsize-1)]
    compartmentTimes[1, :] = [0.0, eDuration[1], eDuration[1]+iDuration[1]]
    orderedFutureEIs = Float64[compartmentTimes[1,2]]
    orderedFutureIRs = Float64[compartmentTimes[1,3]]
    orderedThresholds = sort(pressureThresholds[2:end], rev=true)
    futureSEIndices = sortperm(pressureThresholds, rev=true)
    pop!(futureSEIndices) ##Get rid of initial infective
    t = 0.0
    pressure = 0.0
    nInfectious = 0
    if length(orderedThresholds) > 0
        nextThreshold = pop!(orderedThresholds)
        nextExposed = pop!(futureSEIndices)
    else
        nextThreshold = Inf ##Don't update nextExposed as value irrelevant
    end

    ##Main loop
    while length(orderedFutureEIs)+length(orderedFutureIRs)>0
        nextEITime = length(orderedFutureEIs) > 0 ? orderedFutureEIs[end] : Inf
        nextIRTime = length(orderedFutureIRs) > 0 ? orderedFutureIRs[end] : Inf
        nextEventTime = min(nextEITime, nextIRTime)
        nextEventType = nextEITime < nextIRTime ? 1 : -1 ##Change in number infectious
        nextPressure = pressure + (nextEventTime - t)*nInfectious*λscaled
        if nextThreshold < nextPressure
            timeExposed = t + (nextThreshold - pressure)/(nInfectious*λscaled)
            timeInfectious = timeExposed + eDuration[nextExposed]
            timeRecovered = timeInfectious + iDuration[nextExposed]
            compartmentTimes[nextExposed, :] = [timeExposed, timeInfectious, timeRecovered]
            insertsorted!(orderedFutureEIs, timeInfectious, rev=true)
            insertsorted!(orderedFutureIRs, timeRecovered, rev=true)
            nInfectiousAtExposure[nextExposed] = nInfectious
            pressure = nextThreshold
            if length(orderedThresholds) > 0
                nextThreshold = pop!(orderedThresholds)
                nextExposed = pop!(futureSEIndices)
            else
              nextThreshold = Inf ##Don't update nextExposed as value irrelevant
            end
            t = timeExposed
        else
            pressure = nextPressure
            t = nextEventTime
            nInfectious += nextEventType
            if nextEventType == 1
                pop!(orderedFutureEIs)
            else
                pop!(orderedFutureIRs)
            end
        end
    end
    e = EpidemicSummary(popsize, pDist, eDist, iDist, compartmentTimes, u_pressures, u_eDuration, u_iDuration, pressureThresholds, eDuration, iDuration, nInfectiousAtExposure, pressure, 0.0)
    e.loglikelihood = logLikelihood(e)
    return e
end

"
Returns `SEIRParticles` variable which is a copy of `p` except that its ith particle is a copy of the `indices[i]`th old particle.
"
function shuffle!(p::SEIRParticles, indices::Array{Int, 1})
    epidemics_new = Array(EpidemicSummary, p.nparticles)
    for i in 1:p.nparticles
        epidemics_new[i] = deepcopy(p.epidemics[indices[i]])
    end
    p.epidemics[:] = epidemics_new[:]
    return nothing
end

"
Returns distance between observation time vector `y` and observations specified in `a`.

The distance is the sum of two components.
The first is Euclidean distance between sorted values, excluding any which are infinite in either vector after sorting.
The second is the difference in the number of finite values multiplied by `penalty`.
`penalty` should be chosen so that the second component is much larger.
"
function getDistance(a::SEIRDetails, e::EpidemicSummary; penalty=1000.0)
    n = length(a.yobs)
    robs = sort(a.yobs) ##robs = observed removal times
    rsim = sort(e.compartmentTimes[:,3]) ##rsim = simulated removal times
    psim = sort(e.pressures) ##psim = observed pressure thresholds
    ##Convert to inter-removal times
    if robs[1] < Inf
        robs .-= robs[1]
    end
    if rsim[1] < Inf
        rsim .-= rsim[1]
    end
    ##Find sum square difference of rounded removal times, where both are finite
    nFinite_obs = sum(isfinite(robs))
    nFinite_sim = sum(isfinite(rsim))
    nFinite_both = min(nFinite_obs, nFinite_sim)
    robs_round = map(a.round_times, robs[1:nFinite_both])
    rsim_round = map(a.round_times, rsim[1:nFinite_both])
    dist = sqrt(sum((robs_round .- rsim_round) .^ 2))
    ##Penalise for different numbers of removals
    dist += penalty*abs(nFinite_obs - nFinite_sim)
    ##Penalise for pressures
    if (nFinite_sim > nFinite_obs)
        dist += e.totalPressure*(nFinite_sim-nFinite_obs) - sum(psim[(nFinite_obs+1):nFinite_sim])
    else
        dist += sum(psim[(nFinite_sim+1):nFinite_obs])
    end
    return dist
end

"
Returns vector of Euclidean distances between each particle and the observations.
"
function getDistances(a::SEIRDetails, p::SEIRParticles)
    Float64[getDistance(a, e) for e in p.epidemics]
end

"
Performs an MCMC update on all particles, using slice sampling.

The new pressures and durations are inputed to the `makeEpidemic` function to derive the resulting event times.
"
function updateParticles!(a::SEIRDetails, p::SEIRParticles, θ::Array{Float64, 1})
    n = a.popsize
    λscaled = θ[1]/n
    pDist::Distributions.ContinuousUnivariateDistribution = p.epidemics[1].pDist
    eDist::Distributions.ContinuousUnivariateDistribution = p.epidemics[1].eDist
    iDist::Distributions.ContinuousUnivariateDistribution = p.epidemics[1].iDist
    function dosim(u::Array{Float64, 1})
        u_pressures = u[1:(n-1)]
        u_eDuration = u[n:(2*n-1)]
        u_iDuration = u[(2*n):(3*n-1)]
        e = makeEpidemic(λscaled, u_pressures, u_eDuration, u_iDuration, pDist, eDist, iDist)
        return (e, getDistance(a, e))
    end
    z = zeros(p.nparticles)
    for i in 1:p.nparticles
        ecurr = p.epidemics[i]
        newu, p.epidemics[i], z[i] = slicesample([ecurr.u_pressures; ecurr.u_eDuration; ecurr.u_iDuration], dosim, p.ε, w=2.0*p.maxz)
    end
    p.maxz = maximum(abs(z))
    return nothing
end

"
Return log likelihood of epidemic `e`
"
function logLikelihood(e::EpidemicSummary)
    all(e.pressures .>= 0.0) || return(-Inf)
    ll = 0.0
    for i in 1:e.popsize
        ll += logpdf(e.pDist, e.pressures[i])
        ll += logpdf(e.eDist, e.eDuration[i])
        ll += logpdf(e.iDist, e.iDuration[i])
    end
    return ll
end


text:
"""
    savemethodmodel(path, method, model[; force = false])

Save a trained `model` along with a `method` to `path` for later inference.
Use [`loadmethodmodel`](#) for loading both back into a session. If `path`
already exists, only write to it if `force = true`.

If `model` weights are on a GPU, they will be moved to the CPU before saving
so they can be loaded in a non-GPU environment.

[JLD2.jl](https://github.com/JuliaIO/JLD2.jl) is used for serialization.
"""
function savemethodmodel(path, method::LearningMethod, model; force = false)
    if !force && isfile(path)
        error("$path already exists. Use `force = true` to overwrite.")
    end
    jldsave(string(path); model=cpu(model), method=method)
end

"""
    loadmethodmodel(path) -> (method, model)

Load a trained `model` along with a `method` from `path` that were saved
using [`savemethodmodel`](#).

[JLD2.jl](https://github.com/JuliaIO/JLD2.jl) is used for serialization.
"""
function loadmethodmodel(path)
    isfile(path) || error("\"$path\" is not an existing file.")
    method, model = jldopen(string(path), "r") do f
        f["method"], f["model"]
    end
    return method, model
end


text:
using DCD
using Test

@testset "DCD.File" begin
    dcd = load_dcd("data/test.dcd")
    @test_throws String load_dcd("data/test.txt")
    @test dcd isa DCD.File
    @test natoms(dcd) == 3072
    @test nframes(dcd) == 1000
    @test timestep(dcd) == 200.0
    @test hascell(dcd) == true
    @test DCD.nbytes_header(dcd) == 276
    @test DCD.nbytes_frame(dcd) == 36944
    frame, state = iterate(dcd)
    @test frame isa DCD.Frame
    @test state == 2
    @test length(collect(dcd)) == nframes(dcd)
    @test dcd[begin] isa DCD.Frame
    @test dcd[end] isa DCD.Frame
end

@testset "DCD.Frame" begin
    dcd = load_dcd("data/test.dcd")
    f1 = DCD.Frame(dcd, 1)
    f2 = DCD.Frame(dcd, 2)
    @test cell(f1) == [30.73, 30.73, 30.73, 90.0, 90.0, 90.0]
    @test cell(f2) == [30.73, 30.73, 30.73, 90.0, 90.0, 90.0]
    @test size(positions(f1)) == (3, natoms(dcd))
    @test positions(f1)[1] ≈ 4.3010964
    @test timepassed(f1) == 200.0
end


text:
using NCTiles,NCDatasets

# Climatology example: This example has two parts. First, we take a NetCDF file with climatology_bounds
# and rewrite it as is (6a). Then we are write a new NetCDF file with updated dimensions (6b).
# This example is based on the processing that had to be done to put this data into
# LAS/OpenDap. In order to display LAS prefers the data have dimensions lon,lat(,dep),time and that each
# of these be one dimensional (linear), whereas the data is originally formatted to have dimensions
# i,j(,k),t with 2D latitude and longitude added as variables (curvelinear). Additionally, the attribute
# "positive" needed to be added to the depth dimension. This is climatology data and includes
# the climatology_bounds.

# Point to NetCDF file
examplesdir = joinpath("data","ex6")
indir = joinpath(examplesdir,"input")
savedir = joinpath(examplesdir,"output")
if ~ispath(savedir); mkpath(savedir); end
fldname = "FeT"
fname = joinpath(indir,fldname*".0001.nc")
README = [fldname*" -- Source: Gael Forget; version: alpha."]

# Read in current NetCDF File
ncvars,ncdims,fileatts = readncfile(fname)

# Ex 6a: Rewrite as-is
rm(joinpath(savedir,"ex6a.nc"),force=true)
write(ncvars,joinpath(savedir,"ex6a.nc"),README=README)

# Ex 6b: Rewrite NetCDF file with updated dimensions

# Reworking lat,lon,dep Dimensions
dims = []
dimdict = Dict(["i_" => "lon", "j_" => "lat", "k_" => "dep"])
for d in ncvars[fldname].dims
    global dims
    if haskey(dimdict,d.name[1:end-1])
        newdim = dimdict[d.name[1:end-1]]
        newdimname = join([newdim, d.name[end]],"_")
        attribs = ncvars[newdim].atts
        if newdim == "dep" # Add "positive" => "down" to depth
            attribs["positive"] = "down"
            attribs["standard_name"] = "depth"
        end
        # Grab one-dimensional values for lat/lon
        vals = newdim == "lat" ? ncvars[newdim].values[1,:] : ncvars[newdim].values[:,1]

        # Add updated dimension to dims array
        dims = vcat(dims,NCvar(newdimname,
                                ncvars[newdim].units,
                                ncdims[d.name].dims,
                                vals,
                                attribs,
                                NCDatasets))
        pop!(ncvars,newdim) #remove extra dimension variable
    end
end

# Reworking Time Dimension
dims = vcat(dims,NCvar("tim",
                        ncvars["tim"].units,
                        ncdims["t"].dims,
                        Array{DateTime}(ncvars["tim"].values[:]),
                        ncvars["tim"].atts,
                        NCDatasets))
pop!(ncvars,"tim")

# Replacing dimensions in the NCvars
ncvars[fldname] = NCvar(ncvars[fldname].name,
                                ncvars[fldname].units,
                                dims,
                                ncvars[fldname].values,
                                ncvars[fldname].atts,
                                ncvars[fldname].backend)
ncvars["climatology_bounds"] = NCvar(ncvars["climatology_bounds"].name,
                                        ncvars["climatology_bounds"].units,
                                        [ncdims["tcb"],dims[end]],
                                        ncvars["climatology_bounds"].values,
                                        ncvars["climatology_bounds"].atts,
                                        ncvars["climatology_bounds"].backend)

# Write result to a new NetCDF file
rm(joinpath(savedir,"ex6b.nc"),force=true)
write(ncvars,joinpath(savedir,"ex6b.nc"),README=README)


text:
module Zabbix

# Imports
using HTTP
using JSON

"""
# Zabbix API

#### create the zabbix object

julia>z = ZabbixAPI("http://SERVER_URL/zabbix/api_jsonrpc.php","USERNAME","PASSWORD", true)
ZabbixAPI("http://SERVER_URL/zabbix/api_jsonrpc.php", "USERNAME", "PASSWORD", true, 0, Dict("Content-Type"=>"application/json-rpc"), "2.0")

### Note:

Here the verbosity is set to true by default. You may wish to set up by passing false instead of true.
"""
mutable struct ZabbixAPI

    # To be supplied fields
    server_url::String
    username::String
    password::String

    # By default set, but can be changed while instantiating
    verbose::Bool

    # By default set
    id::Int64
    headers::Dict
    jsonrpc::String

    # Inner constructor
    ZabbixAPI(server_url,username,password,id=1,verbose=true,headers=Dict("Content-Type"=>"application/json-rpc"),
    jsonrpc="2.0") = new(server_url,username,password,id,verbose,headers,jsonrpc)
end

"""
# Gets the apiinfo.version data from zabbix server

#### create the zabbix object

julia>z = ZabbixAPI("http://SERVER_URL/zabbix/api_jsonrpc.php","USERNAME","PASSWORD", true)

#### get the api version

julia>api_version(z)

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

v"3.2.11"

# the type returned is actually a VersionNumber
julia>typeof(Zabbix.api_version(zabbix))

VersionNumber

"""
function api_version(z::ZabbixAPI)
    verbose = z.verbose

    # contruct data
    dict_data = Dict("jsonrpc"=>z.jsonrpc,"id"=>z.id,"method"=>"apiinfo.version","params"=>Dict())
    json_data = JSON.json(dict_data)

    # requests data from zabbix
    if verbose @info("Hitting $(z.server_url) ...") end
    output = HTTP.post(z.server_url,body=json_data,headers=z.headers)

    # return response
    VersionNumber(JSON.parse(String(output.body))["result"])
end

"""
# Gets the authentication token from zabbix server

#### create the zabbix object

julia>z = ZabbixAPI("http://SERVER_URL/zabbix/api_jsonrpc.php","USERNAME","PASSWORD", true)

#### get the auth token

julia> auth_token(z)

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

"bc86891a1c6c9ef5d41d640eb258a81a"

"""
function auth_token(z::ZabbixAPI)
    verbose = z.verbose

    # construct data
    dict_data=Dict("jsonrpc"=>z.jsonrpc,"id"=>z.id,"method"=>"user.login",
            "params"=>Dict("user"=>z.username,"password"=>z.password))
    json_data = JSON.json(dict_data)

    # requests data from zabbix
    if verbose @info("Hitting $(z.server_url) ...") end
    output = HTTP.post(z.server_url,body=json_data,headers=z.headers)

    # return token
    JSON.parse(String(output.body))["result"]

end

"""
# Gets all host for a user

#### create the zabbix object

julia>z = ZabbixAPI("http://SERVER_URL/zabbix/api_jsonrpc.php","USERNAME","PASSWORD", false)

#### get all hosts for a user

julia>get_all_hosts(z)

INFO: Getting authentication token ...

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

Dict{String,Any} with 3 entries:
  "id"      => 1
  "jsonrpc" => "2.0"
  "result"  => Any[Dict{String,Any}(Pair{String,Any}("host", "localhost"),Pair{…

"""
function get_all_hosts(z::ZabbixAPI)
    verbose = z.verbose

    # get token
    if verbose @info("Getting authentication token ... ") end
    token = auth_token(z)

    # construct data
    dict_data = Dict("jsonrpc"=>z.jsonrpc,"id"=>z.id,"method"=>"host.get","auth"=>token,
                        "params"=>Dict("output"=>["hostid","host"],"selectInterfaces"=>["interfaceid","ip"]))
    json_data = JSON.json(dict_data)

    # requests data from zabbix
    if verbose @info("Hitting $(z.server_url) ...") end
    output = HTTP.post(z.server_url,body=json_data,headers=z.headers)

    # return token
    JSON.parse(String(output.body))

end

"""
# Make request to the zabbix server

The make_request function requires you to pass methods(aka Zabbix methods like hosts.get etc) and params ie. parameters in a form of a Dict() object. A easy sample is given on Zabbix's official website
For references on various methods available in Zabbix,
head over to https://www.zabbix.com/documentation/2.2/manual/api/reference

#### create the zabbix object

julia>z = ZabbixAPI("http://SERVER_URL/zabbix/api_jsonrpc.php","USERNAME","PASSWORD", true)

#### create the zabbix method

julia> method = "host.get"

"host.get"

#### construct the params dict object

julia> params = Dict("output"=>"extend", "filter"=>Dict("host"=>["localhost"]))

Dict{String,Any} with 2 entries:
  "output" => "extend"
  "filter" => Dict("host"=>String["localhost"])

#### create the zabbix object

julia>z = ZabbixAPI("http://SERVER_URL/zabbix/api_jsonrpc.php","USERNAME","PASSWORD", true)

#### finally make request

julia> Zabbix.make_request(z, method, params)

INFO: Getting authentication token ...

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

INFO: Updating request id for next API call ...

Dict{String,Any} with 3 entries:
  "id"      => 1
  "jsonrpc" => "2.0"
  "result"  => Any[Dict{String,Any}(Pair{String,Any}("lastaccess", "0"),Pair{String,Any}("ipmi_privilege", "2"),Pair{String,Any}("ipmi_error…

#### display all entries for a host

julia> Zabbix.make_request(z, method, params)["result"][1]

INFO: Getting authentication token ...

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

INFO: Updating request id for next API call ...

Dict{String,Any} with 39 entries:
  "lastaccess"         => "0"
  "ipmi_privilege"     => "2"
  "ipmi_errors_from"   => "0"
  "snmp_available"     => "0"
  "templateid"         => "0"
  "disable_until"      => "0"
  "jmx_available"      => "0"
  "maintenance_from"   => "0"
  "tls_psk_identity"   => ""
  "available"          => "1"
  "ipmi_password"      => ""
  "tls_accept"         => "1"
  "name"               => "localhost"
  "tls_issuer"         => ""
  "status"             => "0"
  "maintenance_status" => "0"
  "hostid"             => "10084"
  "tls_connect"        => "1"
  "ipmi_available"     => "0"
  "description"        => ""
  "errors_from"        => "0"
  "maintenance_type"   => "0"
  "error"              => ""
  "ipmi_username"      => ""
  "snmp_disable_until" => "0"
  "snmp_error"         => ""
  "tls_subject"        => ""
  "maintenanceid"      => "0"
  "host"               => "localhost"
  "jmx_error"          => ""
  "ipmi_disable_until" => "0"
  "snmp_errors_from"   => "0"
  ⋮                    => ⋮

#### get the host id of the zabbix host

julia> Zabbix.make_request(z, method, params)["result"][1]["hostid"]

INFO: Getting authentication token ...

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

INFO: Hitting http://SERVER_URL/zabbix/api_jsonrpc.php ...

INFO: Updating request id for next API call ...

"10084"

"""
function make_request(z::ZabbixAPI, method::String, params=Dict())
    verbose = z.verbose

    if method == "apiinfo.version" && params == Dict()
        return api_version(z)
    else
        # get token
        if verbose @info("Getting authentication token ... ") end
        token = auth_token(z)

        # construct data
        dict_data = Dict("jsonrpc"=>z.jsonrpc,"id"=>z.id,"method"=>method,"auth"=>token,"params"=>params)
        json_data = JSON.json(dict_data)

        # requests data from zabbix
        if verbose @info("Hitting $(z.server_url) ...") end
        output = HTTP.post(z.server_url,body=json_data,headers=z.headers)

        # increment the id for further calls
        if verbose @info("Updating request id for next API call ...") end
        setfield!(z,:id, z.id+1)

        # return token
        JSON.parse(String(output.body))
    end
end

# Exports
export ZabbixAPI, api_version, get_all_hosts, make_request

end


text:
using DifferentialDynamicProgramming
using Test, Statistics, LinearAlgebra

@info("Compile time is high for this package, this is expected and is not an error.")
# write your own tests here

include("test_readme.jl")
demo_linear()
demo_linear_kl(kl_step=100)
demo_pendcart()

demoQP()

# include(Pkg.dir("GuidedPolicySearch","examples","bb.jl"))


text:
using SimulationLogger
using Zygote
using Test
using BangBang

function main()
    println(">"^10 * " log true " * "<"^10)
    _main(true)
    println(">"^10 * " log false " * "<"^10)
    _main(false)
end

"""
log = true or log = false
"""
function _main(log::Bool)
    test_log(log)
    test_nested_log(log)
    test_zygote()
    nothing
end

function test_zygote()
    println("#"^10 * " test_zygote " * "#"^10)
    function my_func(y)
        __LOGGER__ = NamedTuple()
        @show @log y
        @show @log x = y * 2
        result = __LOGGER__
    end
    @show my_func(1.0)  # function works?
    @test gradient(y -> my_func(y) |> sum, 1.0) == (3.0,)
end

function test_log(log)
    println("#"^10 * " test_log " * "#"^10)
    y = 2
    if log
        __LOGGER__ = NamedTuple()
        @show @log y
        @show @log x = y * 2
        @show __LOGGER__
        @test __LOGGER__ == (; y=2, x=4)
    else
        @show @log y
        @show @log x = y * 2
        @test !(@isdefined(__LOGGER__))
    end
end

function test_nested_log(log)
    println("#"^10 * " test_nested_log " * "#"^10)
    y1 = 2
    y2 = 4
    yy = 6
    if log
        __LOGGER__ = NamedTuple()
        @show @nested_log :values x = 1
        @test __LOGGER__ == (; values = (; x = 1))
        # @show @nested_log :values x = 2  # new assignment with the same key will yield an error
        @show @nested_log :values y1, y2
        @test __LOGGER__ == (; values = (; x = 1, y1=2, y2=4))
        @show @nested_log :values yy
        @test __LOGGER__ == (; values = (; x = 1, y1=2, y2=4, yy=6))
        @show @nested_log my_logging()
        @test __LOGGER__ == (; values = (; x = 1, y1=2, y2=4, yy=6, z=1))
        @show _w = @nested_log :good my_logging2()
        @test __LOGGER__ == (; values = (; x = 1, y1=2, y2=4, yy=6, z=1), good= (; w=3))
        @test _w == 3
        @show @nested_onlylog k = 1
        @test __LOGGER__ == (; values = (; x = 1, y1=2, y2=4, yy=6, z=1), good= (; w=3), k=1)
        @show @nested_onlylog :values j = 1
        @test __LOGGER__ == (; values = (; x = 1, y1=2, y2=4, yy=6, z=1, j=1), good= (; w=3), k=1)
        @show k
        @show __LOGGER__
    else
        @show @nested_log :values x = 1
        @show @nested_log :values y1, y2
        @show @nested_log :values yy
        @show @nested_log my_logging()
        @show _w = @nested_log :good my_logging2()
        @show @nested_onlylog :values j = 1  # @nested_onlylog produces `nothing`
        # @show k  # ERROR: UndefVarError: k not defined
        @test !(@isdefined(__LOGGER__))
    end
end

@Loggable function my_logging()
    @nested_log :values z = 1
end

@Loggable function my_logging2()
    @log w = 3
    w
end

@testset "macros" begin
    main()
end


text:
# Profiles transit_poly.jl

r = 0.1  # radius ratio is 0.1
u = [0.2,0.3]; n = length(u)
N_c = n+2

const c_n = zeros(typeof(r),n+3)
const dfdrbc = zeros(typeof(r),n+3)
const a_n = zeros(typeof(r),n+1)
const dadu = zeros(typeof(r),n+1,n)
const dcdu = zeros(typeof(r),n+3,n)
const sn = zeros(typeof(r),N_c+1)
const dsndr = zeros(typeof(r),N_c+1)
const dsndb = zeros(typeof(r),N_c+1)
if iseven(N_c)
  v_max = round(Int64,N_c/2)+2
else
  v_max = round(Int64,(N_c-1)/2)+2
end
const Iv = zeros(typeof(r),v_max+1)
const s2_grad = zeros(typeof(r),2)
const Jv = zeros(typeof(r),v_max+1)
const dIvdk = zeros(typeof(r),v_max+1)
const dJvdk = zeros(typeof(r),v_max+1)

#include("../../../src/transit_poly_prealloc.jl")
include("../../../src/transit_poly_struct.jl")
#include("../../../src/dJv_seriesdk.jl")
using PyPlot

# Now, run the same with transit_poly:

function profile_transit_poly()
nb = 1000001; b=zeros(nb)
flux = zeros(nb); fgrad = zeros(4)
flux_big = zeros(nb); fgrad_big = zeros(BigFloat,4)
flux_grad = zeros(nb,4)
flux_grad_big = zeros(nb,4)
tic()
for i=0:1000000
  b[i+1] = sqrt(((i-500000.0)/500000.0*(1.0+2.0*trans.r))^2)
  flux[i+1]= transit_poly!(r,b[i+1],u,fgrad)
  flux_grad[i+1,:] = fgrad
end
toc()
plot(b,flux)
println("b: ",minimum(b)," ",maximum(b))
println("f: ",minimum(flux)," ",maximum(flux))
return
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
# Conquer with words

module ChainForth

export interpreter, interpret, ForthEngine

const VERSION = "0 dev"

mutable struct ExecutionToken # todo immutable
    code::Union{Function, Vector, Number}
    immediate::Bool
    name::String
    ExecutionToken(name, code = [], immediate = false) = new(code, immediate, name)
end

Base.show(io::IO, w::ExecutionToken) = print(io, w.name)

@enum EngineMode::Int8 MODE_INTERPRET=0 MODE_COMPILE=-1

mutable struct ForthEngine
    input::IO
    out::IO
    stack::Vector{Any}
    memory::Vector{Any} # data-space, linear memory. Words are currently not defined here
    here::Int # "data-space pointer"
    dictionary::Dict{String, ExecutionToken}
    mode::EngineMode
    latest::Union{ExecutionToken, Nothing}
    ForthEngine(input, output) = new(input, output, [], [], 0, Dict(), MODE_INTERPRET, nothing)
end

getword(engine, wordstr) = get(engine.dictionary, wordstr, nothing)

function codeof(engine, wordstr)
    word = getword(engine, wordstr)
    if isnothing(word)
        return tryparse(Int64, wordstr)
    end
    return word
end

function define(engine, wordstr, definition, immediate = false)
    engine.dictionary[wordstr] = ExecutionToken(wordstr, definition, immediate)
    return engine
end

function interpreter(input = stdin, output = stdout)
    return define_stdlib(ForthEngine(input, output))
end

include("stackops.jl")
include("arithmetic.jl")
include("compile.jl")
include("control.jl")
include("memory.jl")
include("time.jl")

function define_env(machine)
    define_stackops(machine)
    define_arithmetic(machine)
    define_controlstructures(machine)
    define_compiler(machine)
    define_memory(machine)
    define_time(machine)
end

function define_stdlib(machine)
    define_env(machine)
    interpret(machine,
        op_ifthenelse * # control.jl TODO: separate the standard lib
        op_memory *
        op_time
    )
    return machine
end

include("interpret.jl")

end # module


text:
"""
Abstract type for defining the strategy for updating the `shift` with
[`update_shift()`](@ref). Implemented strategies:

* [`DontUpdate`](@ref)
* [`DoubleLogUpdate`](@ref) - default in [`lomc!()`](@ref)
* [`LogUpdate`](@ref)
* [`LogUpdateAfterTargetWalkers`](@ref) - FCIQMC standard
* [`DoubleLogUpdateAfterTargetWalkers`](@ref)
"""
abstract type ShiftStrategy end

"""
    update_shift(s <: ShiftStrategy, shift, shiftMode, tnorm, pnorm, dτ, step, df, v_new, v_old)
Update the shift according to strategy `s`. See [`ShiftStrategy`](@ref).
"""
update_shift

"""
    DontUpdate(; targetwalkers = 1_000_000) <: ShiftStrategy
Don't update the `shift`.  Return when `targetwalkers` is reached.
"""
@with_kw struct DontUpdate <: ShiftStrategy
    targetwalkers::Int = 1_000_000
end

@inline function update_shift(s::DontUpdate, shift, _, tnorm, args...)
    return shift, false, tnorm, tnorm < s.targetwalkers
end

"""
    LogUpdateAfterTargetWalkers(targetwalkers, ζ = 0.08) <: ShiftStrategy
Strategy for updating the shift: After `targetwalkers` is reached, update the
shift according to the log formula with damping parameter `ζ`.
See [`LogUpdate`](@ref).
"""
@with_kw struct LogUpdateAfterTargetWalkers <: ShiftStrategy
    targetwalkers::Int
    ζ::Float64 = 0.08 # damping parameter, best left at value of 0.3
end
@inline function update_shift(s::LogUpdateAfterTargetWalkers,
                        shift, shiftMode, tnorm, args...)
    if shiftMode || real(tnorm) > s.targetwalkers
        return update_shift(LogUpdate(s.ζ), shift, true, tnorm, args...)
    end
    return shift, false, tnorm, true
end

"""
    LogUpdate(ζ = 0.08) <: ShiftStrategy
Strategy for updating the shift according to the log formula with damping
parameter `ζ`.

```math
S^{n+1} = S^n -\\frac{ζ}{dτ}\\ln\\left(\\frac{\\|Ψ\\|_1^{n+1}}{\\|Ψ\\|_1^n}\\right)
```
"""
@with_kw struct LogUpdate <: ShiftStrategy
    ζ::Float64 = 0.08 # damping parameter, best left at value of 0.3
end

@inline function update_shift(s::LogUpdate,
                        shift, shiftMode,
                        tnorm, pnorm, dτ, args...)
    # return new shift and new shiftMode
    return shift - s.ζ/dτ * log(tnorm/pnorm), true, tnorm, true
end

"""
    DoubleLogUpdate(; targetwalkers = 1000, ζ = 0.08, ξ = ζ^2/4) <: ShiftStrategy
Strategy for updating the shift according to the log formula with damping
parameter `ζ` and `ξ`.

```math
S^{n+1} = S^n -\\frac{ζ}{dτ}\\ln\\left(\\frac{\\|Ψ\\|_1^{n+1}}{\\|Ψ\\|_1^n}\\right)-\\frac{ξ}{dτ}\\ln\\left(\\frac{\\|Ψ\\|_1^{n+1}}{\\|Ψ\\|_1^\\text{target}}\\right)
```
When ξ = ζ^2/4 this corresponds to critical damping with a damping time scale
T = 2/ζ.
"""
struct DoubleLogUpdate{T} <: ShiftStrategy
    targetwalkers::T
    ζ::Float64 # damping parameter, best left at value of 0.08
    ξ::Float64  # restoring force to bring walker number to the target
end
function DoubleLogUpdate(;targetwalkers = 1000,  ζ = 0.08, ξ = ζ^2/4)
    return DoubleLogUpdate(targetwalkers, ζ, ξ)
end

@inline function update_shift(s::DoubleLogUpdate,
                        shift, shiftMode,
                        tnorm, pnorm, dτ, args...)
    new_shift = shift - s.ξ/dτ * log(tnorm/s.targetwalkers) - s.ζ/dτ * log(tnorm/pnorm)
    return new_shift, true, tnorm, true
end

"""
    DoubleLogUpdateAfterTargetWalkers(targetwalkers, ζ = 0.08, ξ = 0.0016) <: ShiftStrategy
Strategy for updating the shift: After `targetwalkers` is reached, update the
shift according to the log formula with damping parameter `ζ` and `ξ`.
See [`DoubleLogUpdate`](@ref).
"""
@with_kw mutable struct DoubleLogUpdateAfterTargetWalkers <: ShiftStrategy
    targetwalkers::Int
    ζ::Float64 = 0.08 # damping parameter, best left at value of 0.3
    ξ::Float64 = 0.0016 # restoring force to bring walker number to the target
end

@inline function update_shift(s::DoubleLogUpdateAfterTargetWalkers,
                        shift, shiftMode,
                        tnorm, pnorm, dτ, step, df, args...)
    if shiftMode || real(tnorm) > s.targetwalkers
            return update_shift(
                DoubleLogUpdate(s.targetwalkers,s.ζ,s.ξ),
                shift,
                true,
                tnorm,
                pnorm,
                dτ,
                step,
                df
            )
    end
    return shift, false, tnorm, true
end

# more experimental strategies from here on:

"""
    DoubleLogSumUpdate(; targetwalkers = 1000, ζ = 0.08, ξ = ζ^2/4, α = 1/2) <: ShiftStrategy
Strategy for updating the shift according to the log formula with damping
parameters `ζ` and `ξ`.

```math
S^{n+1} = S^n -\\frac{ζ}{dτ}\\ln\\left(\\frac{N_\\mathrm{w}^{n+1}}{N_\\mathrm{w}^n}\\right)
- \\frac{ξ}{dτ}\\ln\\left(\\frac{N_\\mathrm{w}^{n+1}}{N_\\mathrm{w}^\\text{target}}\\right),
```
where ``N_\\mathrm{w} =`` `(1-α)*walkernumber() + α*UniformProjector()⋅ψ` computed with
[`walkernumber()`](@ref) and [`UniformProjector()`](@ref).
When ξ = ζ^2/4 this corresponds to critical damping with a damping time scale
T = 2/ζ.
"""
struct DoubleLogSumUpdate{T} <: ShiftStrategy
    targetwalkers::T
    ζ::Float64 # damping parameter, best left at value of 0.08
    ξ::Float64  # restoring force to bring walker number to the target
    α::Float64  # mixing angle for (1-α)*walkernumber + α*UniformProjector()⋅ψ
end
function DoubleLogSumUpdate(;targetwalkers = 1000,  ζ = 0.08, ξ = ζ^2/4, α = 1/2)
    DoubleLogSumUpdate(targetwalkers,  ζ, ξ, α)
end

@inline function update_shift(s::DoubleLogSumUpdate,
                        shift, shiftMode,
                        tnorm, pnorm, dτ, step, df, v_new, v_old
)
    tp = DictVectors.UniformProjector() ⋅ v_new
    pp = DictVectors.UniformProjector() ⋅ v_old
    twn = (1 - s.α) * tnorm + s.α * tp
    pwn = (1 - s.α) * pnorm + s.α * pp
    # return new shift and new shiftMode
    new_shift = shift - s.ξ/dτ * log(twn/s.targetwalkers) - s.ζ/dτ * log(twn/pwn)
    return new_shift, true, tnorm, true
end

"""
    TripleLogUpdate(; targetwalkers = 1000, ζ = 0.08, ξ = ζ^2/4, η = 0.01) <: ShiftStrategy
Strategy for updating the shift according to the extended log formula with damping
parameters `ζ`, `ξ`, and `η`.

```math
S^{n+1} = S^n -\\frac{ζ}{dτ}\\ln\\left(\\frac{N_\\mathrm{w}^{n+1}}{N_\\mathrm{w}^n}\\right)
- \\frac{ξ}{dτ}\\ln\\left(\\frac{N_\\mathrm{w}^{n+1}}{N_\\mathrm{w}^\\text{target}}\\right)
- \\frac{η}{dτ}\\ln\\left(\\frac{\\|ℜ(Ψ^{n+1})\\|_1^2 + \\|ℑ(Ψ^{n+1})\\|_1^2}
{\\|ℜ(Ψ^{n})\\|_1^2 + \\|ℑ(Ψ^{n})\\|_1^2}\\right),
```
where ``N_\\mathrm{w}`` is the [`walkernumber()`](@ref).
When ξ = ζ^2/4 this corresponds to critical damping with a damping time scale
T = 2/ζ.
"""
struct TripleLogUpdate{T} <: ShiftStrategy
    targetwalkers::T
    ζ::Float64 # damping parameter
    ξ::Float64  # restoring force to bring walker number to the target
    η::Float64
end
function TripleLogUpdate(;targetwalkers = 1000,  ζ = 0.08, ξ = ζ^2/4, η = 0.01)
    return TripleLogUpdate(targetwalkers, ζ, ξ, η)
end

@inline function update_shift(s::TripleLogUpdate,
                        shift, shiftMode,
                        tnorm, pnorm, dτ, step, df, v_new, v_old
)
    tp = abs2(DictVectors.Norm1ProjectorPPop() ⋅ v_new)
    pp = abs2(DictVectors.Norm1ProjectorPPop() ⋅ v_old)
    # return new shift and new shiftMode
    new_shift = shift - s.ξ/dτ * log(tnorm/s.targetwalkers) - s.ζ/dτ * log(tnorm/pnorm)
    # new_shift -= s.η/dτ * log(tp/pp)
    new_shift -= s.η/dτ * log(tp/s.targetwalkers)
    return new_shift, true, tnorm, true
end

"""
    DoubleLogProjected(; target, projector, ζ = 0.08, ξ = ζ^2/4) <: ShiftStrategy
Strategy for updating the shift according to the log formula with damping
parameter `ζ` and `ξ` after projecting onto `projector`.

```math
S^{n+1} = S^n -\\frac{ζ}{dτ}\\ln\\left(\\frac{P⋅Ψ^{(n+1)}}{P⋅Ψ^{(n)}}\\right)-\\frac{ξ}{dτ}\\ln\\left(\\frac{P⋅Ψ^{(n+1)}}{\\text{target}}\\right)
```
"""
struct DoubleLogProjected{P} <: ShiftStrategy
    target::Float64
    projector::P
    ζ::Float64 # damping parameter, best left at value of 0.08
    ξ::Float64 # restoring force to bring walker number to the target
end
DoubleLogProjected(; target, projector, ζ = 0.08, ξ = ζ^2/4) = DoubleLogProjected(target, projector, ζ, ξ)

@inline function update_shift(s::DoubleLogProjected,
                        shift, shiftMode,
                        tnorm, pnorm, dτ, step, df, v_new, v_old)
    # return new shift and new shiftMode
    tp = s.projector⋅v_new
    pp = s.projector⋅v_old
    new_shift = shift - s.ζ/dτ * log(tp/pp) - s.ξ/dτ * log(tp/s.target)
    return new_shift, true, tnorm, true
end


text:
using Dates: DatePeriod, TimePeriod, firstdayofweek, lastdayofweek,
    firstdayofmonth, lastdayofmonth, firstdayofyear, lastdayofyear,
    firstdayofquarter, lastdayofquarter

# Truncation
# TODO: Just utilize floor code for truncation?
function Base.trunc(zdt::ZonedDateTime, ::Type{P}) where P <: DatePeriod
    ZonedDateTime(trunc(DateTime(zdt), P), timezone(zdt))
end
function Base.trunc(zdt::ZonedDateTime, ::Type{P}) where P <: TimePeriod
    local_dt = trunc(DateTime(zdt), P)
    utc_dt = local_dt - zdt.zone.offset
    ZonedDateTime(utc_dt, timezone(zdt); from_utc=true)
end
Base.trunc(zdt::ZonedDateTime, ::Type{Millisecond}) = zdt

# Adjusters
for prefix in ("firstdayof", "lastdayof"), suffix in ("week", "month", "year", "quarter")
    func = Symbol(prefix * suffix)
    @eval function Dates.$func(zdt::ZonedDateTime)
        ZonedDateTime($func(DateTime(zdt)), zdt.timezone)
    end
end


text:
"""
`_spectral(G)` Gives the graph held an
embedding based on the eigenvectors of the Laplacian matrix of the
graph. Specifically, the `x`-coordinates come from the eigenvector
associated with the second smallest eigenvalue, and the
`y`-coordinates come from the eigenveector associated with the third
smallest.

This may also be invoked as `_spectral(G,xcol,ycol)` to choose other
eigenvectors to use for the x and y coordinates of the embedding.
"""
function _spectral(G::SimpleGraph, xcol::Int = 2, ycol::Int = 3)
    L = laplace(G)
    EV = eigvecs(L)
    x = EV[:, xcol]
    y = EV[:, ycol]

    VV = vlist(G)
    n = length(VV)

    for k = 1:n
        v = VV[k]
        G.cache[:xy][v] = [x[k], y[k]]
    end
    scale(G)
end


text:
using LearnBase
using SparseArrays
using Test

include("aggmode.jl")
include("obsdim.jl")
include("costs.jl")
include("other.jl")


text:
module NoveltyFiltering

using Scruff
using Scruff.Models
using Scruff.SFuncs
using Scruff.Algorithms
import Scruff: make_initial, make_transition

struct NoveltySetup
    known_velocities::Vector{Float64}
    known_probs::Vector{Float64}
    novelty_prob::Float64
    novelty_prior_mean::Float64
    novelty_prior_sd::Float64
    transition_sd::Float64
    observation_sd::Float64
end

struct PositionModel <: VariableTimeModel{Tuple{}, Tuple{Float64, Float64}, Float64}
    setup::NoveltySetup
end
function make_initial(::PositionModel, ::Float64)::Dist{Float64}
    return Constant(0.0)
end
function make_transition(posmod::PositionModel, parenttimes::Tuple{Float64, Float64}, time::Float64)::SFunc{Tuple{Float64, Float64}, Float64}
    function f(pair)
        (prevval, velocity) = pair
        Normal(prevval + t * velocity, t * posmod.setup.transition_sd)
    end
    t = time - parenttimes[1]
    return Chain(Tuple{Float64, Float64}, Float64, f)
end

function novelty_network(setup::NoveltySetup, numobs::Int)::DynamicNetwork
    known_velocity = StaticModel(Cat(setup.known_velocities, setup.known_probs))(:known_velocity)
    is_novel = StaticModel(Flip(setup.novelty_prob))(:is_novel)
    novel_velocity = StaticModel(Normal(setup.novelty_prior_mean, setup.novelty_prior_sd))(:novel_velocity)
    velocity = StaticModel(If{Float64}())(:velocity)

    position = PositionModel(setup)(:position)
    observation = SimpleModel(LinearGaussian((1.0,), 0.0, setup.observation_sd))(:observation)

    variables = [known_velocity, is_novel, novel_velocity, velocity, position, observation]
    initial_graph = VariableGraph(velocity => [is_novel, novel_velocity, known_velocity], observation => [position])
    transition_graph = VariableGraph(known_velocity => [known_velocity], is_novel => [is_novel], novel_velocity => [novel_velocity],
                                     velocity => [velocity], position => [position, velocity], observation => [position])

    return DynamicNetwork(variables, initial_graph, transition_graph)
end

obsname(i) = Symbol("obs", i)

function do_experiment(setup::NoveltySetup, obs::Vector{Tuple{Float64, Float64}}, alg::Filter)
    net = novelty_network(setup, length(obs))
    runtime = Runtime(net, 0.0) # Set the time type to Float64 and initial time to 0
    init_filter(alg, runtime)

    is_novel = get_node(net, :is_novel)
    velocity = get_node(net, :velocity)
    observation = get_node(net, :observation)

    for (time, x) in obs
        evidence = Dict{Symbol, Score}(:observation => HardScore(x))
        println("Observing ", x, " at time ", time)
        # At a minimum, we need to include query and evidence variables in the filter step
        filter_step(alg, runtime, Variable[is_novel, velocity, observation], time, evidence)

        println("Probability of novel = ", probability(alg, runtime, is_novel, true))
        println("Posterior mean of velocity = ", mean(alg, runtime, velocity))
    end
end

# Known velocities are 0 and 1, novelty has mean 0 and standard deviation 10
setup = NoveltySetup([0.0, 1.0], [0.7, 0.3], 0.1, 0.0, 10.0, 1.0, 1.0)
obs1 = [(1.0, 2.1), (3.0, 5.8), (3.5, 7.5)] # consistent with velocity 2
obs2 = [(1.0, 4.9), (3.0, 17.8), (3.5, 20.5)] # consistent with velocity 6

println("Particle filter")
println("Smaller velocity")
# CoherentPF is a kind of asynchronous filter that makes sure all relevant variables are included in the filter step.
# In this example, it will ensure that position is included.
do_experiment(setup, obs1, CoherentPF(1000))
println("\nLarger velocity")
do_experiment(setup, obs2, CoherentPF(1000))

# The BP filter is included to show it's use, but it doesn't work well here.
# It only passes marginal distributions from one time step to the next,
# but with static variables, joint distributions are needed.
println("\nBP filter")
println("Smaller velocity")
do_experiment(setup, obs1, CoherentBP())
println("\nLarger velocity")
do_experiment(setup, obs2, CoherentBP())

end


text:
# import Unreal: unreal

"""
    Eisenstein{T <: Real} <: Number

An immutable pair of real numbers describing an Eisenstein number.
"""
struct Eisenstein{T <: Real} <: Number
    l::T
    r::T

    Eisenstein{T}(l::T, r::T) where {T <: Real} = new(l, r)
end

const Eisenstein16 = Eisenstein{Int8}
const Eisenstein32 = Eisenstein{Int16}
const Eisenstein64 = Eisenstein{Int32}
const Eisenstein128 = Eisenstein{Int64}
const Eisenstein256 = Eisenstein{Int128}
const EisensteinBigInt = Eisenstein{BigInt}

# Constructors
function Eisenstein(a::T, b::T) where T <: Real
    Eisenstein{T}(a, b)
end

function Eisenstein(a::T) where T <: Real
    Eisenstein{T}(a, zero(T))
end

function Eisenstein(a::Real, b::Real)
    Eisenstein(promote(a, b)...)
end

# Eisenstein unit
const ω = Eisenstein(false, true)

# Real part
function real(z::Eisenstein)
    z.l - (z.r / 2)
end

# Imaginary part
function imag(z::Eisenstein)
    z.r * halfsqrt3
end

function Complex(z::Eisenstein)
    Complex(real(z), imag(z))
end

function complex(z::Eisenstein)
    Complex(z)
end

# Unreal part
function unreal(z::Eisenstein)
    z.r
end

function isreal(z::Eisenstein)
    iszero(z.r)
end

function asarray(z::Eisenstein)
    vcat(z.l, z.r)
end

function iszero(z::Eisenstein)
    iszero(asarray(z))
end

function zero(::Type{Eisenstein{T}}) where T <: Real
    Eisenstein{T}(zero(T), zero(T))
end

function zero(::Eisenstein{T}) where T <: Real
    zero(Eisenstein{T})
end

function one(::Type{Eisenstein{T}}) where T <: Real
    Eisenstein{T}(one(T), zero(T))
end

function one(::Eisenstein{T}) where T <: Real
    one(Eisenstein{T})
end

function omega(::Type{Eisenstein{T}}) where T <: Real
    Eisenstein{T}(zero(T), one(T))
end

function omega(::Eisenstein{T}) where T <: Real
    omega(Eisenstein{T})
end

function show(io::IO, z::Eisenstein)
    print(io, "[1: ")
    print(io, z.l)
    print(io, ", ω: ")
    print(io, z.r)
    print(io, "]")
end

function conj(z::Eisenstein)
    Eisenstein(z.l - z.r, -(z.r))
end

function (+)(z::Eisenstein)
    z
end

function (+)(x::Eisenstein, y::Eisenstein)
    Eisenstein(x.l + y.l, x.r + y.r)
end

function (+)(a::Real, z::Eisenstein)
    Eisenstein(a + z.l, z.r)
end

function (+)(z::Eisenstein, a::Real)
    Eisenstein(z.l + a, z.r)
end

function (-)(z::Eisenstein)
    Eisenstein(-(z.l), -(z.r))
end

function (-)(x::Eisenstein, y::Eisenstein)
    Eisenstein(x.l - y.l, x.r - y.r)
end

function (-)(a::Real, z::Eisenstein)
    Eisenstein(a - z.l, -(z.r))
end

function (-)(z::Eisenstein, a::Real)
    Eisenstein(z.l - a, z.r)
end

function (*)(x::Eisenstein, y::Eisenstein)
    Eisenstein(
        (x.l * y.l) - (x.r * y.r),
        (x.l * y.r) + (x.r * y.l) - (x.r * y.r),
    )
end

function (*)(a::Real, z::Eisenstein)
    Eisenstein(a * z.l, a * z.r)
end

function (*)(z::Eisenstein, a::Real)
    Eisenstein(z.l * a, z.r * a)
end

function abs2(z::Eisenstein)
    (z.l)^(2) + (z.r)^(2) - (z.l * z.r)
end

function inv(z::Eisenstein)
    if iszero(z)
        error(ZeroInverse)
    end

    conj(z) / abs2(z)
end

function (/)(x::Eisenstein, y::Eisenstein)
    x * inv(y)
end

function (\)(y::Eisenstein, x::Eisenstein)
    inv(y) * x
end

function (/)(a::Real, z::Eisenstein)
    a * inv(z)
end

function (\)(z::Eisenstein, a::Real)
    inv(z) * a
end

function (/)(z::Eisenstein, a::Real)
    if iszero(a)
        error(ZeroDenominator)
    end

    Eisenstein(z.l / a, z.r / a)
end

function (\)(a::Real, z::Eisenstein)
    if iszero(a)
        error(ZeroDenominator)
    end

    Eisenstein(a \ z.l, a \ z.r)
end

function associates(z::Eisenstein{T}) where T <: Real
    a = z
    b = -(z)
    u = omega(Eisenstein{T})
    c = z * u
    u = -(u)
    d = z * u
    u = u * u
    f = z * u
    u = -(u)
    g = z * u
    a, b, c, d, f, g
end

function random(::Type{Eisenstein{T}}) where T <: Real
    Eisenstein{T}(rand(T), rand(T))
end


text:
function DecimateSym!(var::PlmVar, alg::PlmAlg, dec::DecVar{1})

    Jvec,pseudolike = MinimizePLSym(alg, var)
    Nc2 = binomial(var.N,2)
    LL = Nc2 * var.q2
#    results = (Int,Float64,Float64,Float64)[]
    iteration = 1
    while true
        SetAllCouplingToZero!(Jvec,  var, dec)
        nzeros = sum(!dec.dmask)
        fraczeros = nzeros / LL
        println("# set to zero ", nzeros ," over ", LL ," equal to ", 100.0*fraczeros," % of the couplings. Converging ...")
        pseudolike, Jvec = MinimizePLSym(alg, var, dec, Jvec, x0=Jvec)
        println("# Converged!")
#        push!(results,(iteration,1.0-fraczeros, fmin,fminloc))
        fraczeros >= dec.fracmax && break
        iteration += 1
    end
    return Jvec, pseudolike
end

function SetAllCouplingToZero!(vecJ::Array{Float64,1}, var::PlmVar, dec::DecVar{1})

    N = var.N
    q  = var.q
    q2 = var.q2

    LL = length(vecJ)
    @assert LL == (binomial(N,2)*q2 + N*q)
    println("LL = $LL dmask = ", length(dec.dmask))

@compat    vecval = Tuple{Int,Float64}[]
@compat    sizehint!(vecval,LL)
    for i=1:LL-N*q
        dec.dmask[i] && push!(vecval,(i, vecJ[i]*vecJ[i]))
    end

    lvec = length(vecval)
    idxperm = sortperm(Float64[ vecval[i][2] for i=1:lvec])

@compat    ndec = round(Int, dec.fracdec * (LL-N*q))
    for i=1:ndec
        site = vecval[idxperm[i]][1]
        dec.dmask[site] = false
    end
end

function SetCouplingToZero!(FN::Array{Float64,2},var::PlmVar, dec::DecVar{1})

    N = size(FN,1)

    q  = var.q
    q2 = var.q2

@compat    vecval = Tuple{Int,Float64}[]
    counter = 0
    for i=1:(N-1), j=(i+1):N
        counter += 1
        dec.dmask[counter] && push!(vecval,(counter, FN[i,j]))
    end

    lvec = length(vecval)
    idxperm = sortperm(Float64[ vecval[i][2] for i=1:lvec])

@compat    ndec = round(Int,dec.fracdec * counter)
    for i=1:ndec
        site = vecval[idxperm[i]][1]
        dec.dmask[site] = false
    end
end

function MinimizePLSym(alg::PlmAlg, var::PlmVar, dec::DecVar{1}, Jvec::Array{Float64,1}; x0 = zeros(Float64,length(Jvec)))

    N  = var.N
    q  = var.q
    q2 = var.q2
    Nc2 = binomial(N,2)

    function f(x::Vector, g::Vector)
        eltime = @elapsed begin
        g === nothing && (g = zeros(Float64, length(x)))
        a = PLsiteAndGradSym!(x, g, var)
        end
        println("time spent in iteration = $eltime . PL = $a")
        return a
    end

    opt = Opt(alg.method, length(x0))

    ftol_abs!(opt, alg.epsconv)

    lb,ub = ComputeConstraintsSym!(dec, var.q, var.N, x0)

    lower_bounds!(opt, lb)

    upper_bounds!(opt, ub)

    maxeval!(opt, alg.maxit)

    min_objective!(opt, f)
    elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)

    alg.verbose && @printf("pl = %.4f\t time = %.4f\n",  minf, elapstime)
    return minf,minx
end

function ComputeConstraintsSym!(decvar::DecVar{1}, q::Int, N::Int, x0::Array{Float64,1} )

    q2 = q*q
    dmask = decvar.dmask
    Nc2 = binomial(N,2)
    lb = -Inf * ones(Float64, Nc2*q2 + N*q)
    ub =  Inf * ones(Float64, Nc2*q2 + N*q)
    LL = length(x0)

    for i=1:LL-N*q
        if dmask[i] == false
            lb[ i ] = -1.0e-6
            ub[ i ] =  1.0e-6
            x0[ i ] =  0.0
        else
            x0[ i ] += (1.0 - 2.0*rand())*1e-3
        end
    end

    return lb,ub
end


text:
"""
    mergesorted(x::Vector{T}, y::T)

"""
function mergesorted(x::Vector{T}, y::T) where T
    z, indexes = mergesorted(x, [y])
    return z, indexes[1]
end

"""
    mergesorted(x::Vector{T}, y::Vector{S})

"""
function mergesorted(x::Vector{T}, y::Vector{S}) where {T,S}
    nx = length(x)
    ny = length(y)
    nz = nx + ny
    z  = Array{promote_type(T,S)}(undef, nz)

    ind_x   = 1
    ind_y   = 1
    indexes = ones(Int64, ny)

    @inbounds for ind_z in 1:nz
        val_x = x[ind_x]
        val_y = y[ind_y]
        if val_x <= val_y
            z[ind_z]  = val_x
            ind_x    += 1
            if ind_x > nx
                indexes[ind_y:end] .= (ind_z + 1):nz
                z[(ind_z + 1):end] .= y[ind_y:end]
                break
            end
        else
            z[ind_z]        = val_y
            indexes[ind_y]  = ind_z
            ind_y          += 1
            if ind_y > ny
                z[(ind_z+1):end] .= x[ind_x:end]
                break
            end
        end
    end
    return z, indexes
end

auc_trapezoidal((x,y)::Tuple{RealVector, RealVector}) = auc_trapezoidal(x, y)

"""
    auc(x::RealVector, y::RealVector)

Computes the area under curve `(x,y)` using trapezoidal rule.
"""
function auc_trapezoidal(x::RealVector, y::RealVector)
    n = length(x)
    n == length(y) || throw(DimensionMismatch("Inconsistent lengths of `x` and `y`."))

    if issorted(x)
        prm = 1:n
    elseif issorted(x, rev = true)
        prm = n:-1:1
    else
        throw(ArgumentError("x must be sorted."))
    end

    val = zero(promote_type(eltype(x), eltype(y)))
    for i in 2:n
        Δx   = x[prm[i]]  - x[prm[i-1]]
        Δy   = (y[prm[i]] + y[prm[i-1]])/2
        if !(isnan(Δx) || isnan(Δy) || Δx == 0)
            val += Δx*Δy
        end
    end
    return val
end


text:
function part1(input)
    lines = parse_input(input)
    filter!(line -> line[1] == line[3] || line[2] == line[4], lines)
    map = place_lines(lines)
    return count(>(1), values(map))
end

function part2(input)
    lines = parse_input(input)
    map = place_lines(lines)
    return count(>(1), values(map))
end

parse_input(input) = parse_line.(readlines(input))
parse_line(line) =  parse.(Int, split(replace(line, " -> " => ","), ","))

function place_lines(lines)
    map = Dict{Tuple{Int, Int}, Int}()
    for (x1, y1, x2, y2) in lines
        Δx = sign(x2 - x1)
        Δy = sign(y2 - y1)
        for i = 0:max(abs(x2 - x1), abs(y2 - y1))
            coord = (x1 + i * Δx, y1 + i * Δy)
            map[coord] = get(map, coord, 0) + 1
        end
    end
    return map
end


text:
function _map_to_interval_if_needed(model, x::VI, y::VI, l::Float64, u::Float64)
    if l == 0 && u == 1
        return x, y
    else
        Δ = u - l
        @assert Δ > 0
        x̃_v, x̃_c = MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0))
        ỹ_v, ỹ_c = MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0))
        MOIU.normalize_and_add_constraint(
            model,
            SV(x) - (l + Δ * SV(x̃_v)),
            MOI.EqualTo(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            SV(y) - (l^2 + 2l * Δ * SV(x̃_v) + Δ^2 * SV(ỹ_v)),
            MOI.EqualTo(0.0),
        )
        return x̃_v, ỹ_v
    end
end

# Collects all (scalar) quadratic constraints, normalizes them to form q >= 0,
# vectorizes them, deletes the original quadratic constraints, and then returns
# the vector of q's.
function _collect_and_canonicalize_quadratic_constraints!(
    model::MOIU.Model{Float64},
)
    scalar_quad_cons = model.moi_scalarquadraticfunction
    @assert isempty(scalar_quad_cons.moi_equalto)
    @assert isempty(scalar_quad_cons.moi_interval)
    @assert isempty(scalar_quad_cons.moi_semicontinuous)
    @assert isempty(scalar_quad_cons.moi_semiinteger)
    @assert isempty(scalar_quad_cons.moi_zeroone)

    num_gt = length(scalar_quad_cons.moi_greaterthan)
    num_lt = length(scalar_quad_cons.moi_lessthan)

    constants = Array{Float64}(undef, 0)
    affine_terms = Array{VAT}(undef, 0)
    quad_terms = Array{VQT}(undef, 0)

    func_count = 0
    for (q_ci, q_f, q_s) in scalar_quad_cons.moi_greaterthan
        func_count += 1
        push!(constants, -(q_f.constant - q_s.lower))
        for aff in q_f.affine_terms
            push!(
                affine_terms,
                VAT(func_count, SAT(-aff.coefficient, aff.variable_index)),
            )
        end
        for quad in q_f.quadratic_terms
            push!(
                quad_terms,
                VQT(
                    func_count,
                    SQT(
                        -quad.coefficient,
                        quad.variable_index_1,
                        quad.variable_index_2,
                    ),
                ),
            )
        end
    end

    @assert func_count == num_gt

    for (q_ci, q_f, q_s) in scalar_quad_cons.moi_lessthan
        func_count += 1
        push!(constants, q_f.constant - q_s.upper)
        for aff in q_f.affine_terms
            push!(
                affine_terms,
                VAT(func_count, SAT(aff.coefficient, aff.variable_index)),
            )
        end
        for quad in q_f.quadratic_terms
            push!(
                quad_terms,
                VQT(
                    func_count,
                    SQT(
                        quad.coefficient,
                        quad.variable_index_1,
                        quad.variable_index_2,
                    ),
                ),
            )
        end
    end

    @assert func_count == num_gt + num_lt

    empty!(scalar_quad_cons.moi_greaterthan)
    empty!(scalar_quad_cons.moi_lessthan)
    return VQF(affine_terms, quad_terms, constants)
end

function _indices_and_coefficients(
    I::AbstractVector{Int},
    J::AbstractVector{Int},
    V::AbstractVector{Float64},
    indices::AbstractVector{Int},
    coefficients::AbstractVector{Float64},
    f::SQF,
    canonical_index::Dict{VI,Int},
)
    variables_seen = 0
    for (i, term) in enumerate(f.quadratic_terms)
        vi_1 = term.variable_index_1
        vi_2 = term.variable_index_2
        if !haskey(canonical_index, vi_1)
            variables_seen += 1
            canonical_index[vi_1] = variables_seen
        end
        if !haskey(canonical_index, vi_2)
            variables_seen += 1
            canonical_index[vi_2] = variables_seen
        end

        I[i] = canonical_index[term.variable_index_1]
        J[i] = canonical_index[term.variable_index_2]
        V[i] = term.coefficient
        # Gurobi returns a list of terms. MOI requires 0.5 x' Q x. So, to get
        # from
        #   Gurobi -> MOI => multiply diagonals by 2.0
        #   MOI -> Gurobi => multiply diagonals by 0.5
        # Example: 2x^2 + x*y + y^2
        #   |x y| * |a b| * |x| = |ax+by bx+cy| * |x| = 0.5ax^2 + bxy + 0.5cy^2
        #           |b c|   |y|                   |y|
        #   Gurobi needs: (I, J, V) = ([0, 0, 1], [0, 1, 1], [2, 1, 1])
        #   MOI needs:
        #     [SQT(4.0, x, x), SQT(1.0, x, y), SQT(2.0, y, y)]
        if I[i] == J[i]
            V[i] *= 0.5
        end
    end
    for (i, term) in enumerate(f.affine_terms)
        indices[i] = term.variable_index.value
        coefficients[i] = term.coefficient
    end
    return
end

function _indices_and_coefficients(f::SQF, canonical_index::Dict{VI,Int})
    f_canon = MOI.Utilities.canonical(f)
    nnz_quadratic = length(f_canon.quadratic_terms)
    nnz_affine = length(f_canon.affine_terms)
    I = Vector{Int}(undef, nnz_quadratic)
    J = Vector{Int}(undef, nnz_quadratic)
    V = Vector{Float64}(undef, nnz_quadratic)
    indices = Vector{Int}(undef, nnz_affine)
    coefficients = Vector{Float64}(undef, nnz_affine)
    _indices_and_coefficients(
        I,
        J,
        V,
        indices,
        coefficients,
        f_canon,
        canonical_index,
    )
    return indices, coefficients, I, J, V
end

function _get_Q_matrix(q_func::SQF, canonical_index::Dict{VI,Int}, n::Int)
    indices, coefficients, I, J, V =
        _indices_and_coefficients(q_func, canonical_index)
    Q = Matrix{Float64}(SparseArrays.sparse(I, J, V, n, n))
    Q = 1 / 2 * (Q + Q')

    @assert n == size(Q, 1) == size(Q, 2)
    @assert LinearAlgebra.issymmetric(Q)
    return Q
end


text:
struct EnumSASPermission
    READ::String
    ADD::String
    CREATE::String
    WRITE::String
    DELETE::String
    LIST::String
    PROCESS::String
end

const SASPermission = EnumSASPermission("r", "a", "c", "w", "d", "l", "p")

sas_permissions(perms...) = join(perms)

const SAS_VERSION = "2017-11-09"

sas_date_format(d::DateTime) = Dates.format(d, DateFormat("yyyy-mm-dd\\THH:MM:SS")) * "Z" #*"UTC"

sas_canonicalize(raw_url::String) = sas_canonicalize(URIs.URI(raw_url))
function sas_canonicalize(url::URIs.URI)
    urlparts = split(url.host, ".")
    accountname = String(urlparts[1])
    restype = String(urlparts[2])
    respath = url.path[2:end]

    canonicalized_url = "/" * joinpath(restype, accountname, respath)
    canonicalized_url, accountname, restype
end

sas_filename(raw_url::String) = basename(URIs.URI(raw_url).path)
sas_content_disposition(raw_url::String) = "attachment; filename=" * sas_filename(raw_url)

function appendSAS(ctx, subscription_id::String, resource_group_name::String, raw_url::String; kwargs...)
    accountkey = extract_account_and_key(ctx, subscription_id, resource_group_name, raw_url)
    appendSAS(raw_url, accountkey.key; kwargs...)
end

function appendSAS(raw_url::String, key::String; permissions::String=SASPermission.READ,
        start::Union{Nothing,DateTime}=nothing, expiry_hours::Int=24, expiry::Union{Nothing,DateTime}=nothing,
        identifier::String="", ip::String="",
        cache_control::String="no-cache, no-store, max-age=0",
        content_disposition::String=sas_content_disposition(raw_url),
        content_encoding::String="",
        content_language::String="",
        content_type::String="application/octet-stream")
    url = URIs.URI(raw_url)
    canonicalized_url, accountname, restype = sas_canonicalize(url)

    if restype == "blob"
        sr = "b"
    else
        error("Not implemented yet for resources of type $restype")
    end

    if start === nothing
        st = ""
        start = Dates.now(Dates.UTC)
    else
        st = sas_date_format(start)
    end

    if expiry === nothing
        se = sas_date_format(start + Dates.Hour(expiry_hours))
    else
        se = sas_date_format(expiry)
    end

    str_to_sign = join([
                        permissions, st, se,
                        canonicalized_url,
                        identifier,
                        ip,
                        url.scheme,
                        SAS_VERSION,
                        cache_control,
                        content_disposition,
                        content_encoding,
                        content_language,
                        content_type
                    ], "\n")
    sig = base64encode(digest(MD_SHA256, str_to_sign, base64decode(key)))

    query_params =  [
                        "sv"=>SAS_VERSION,
                        "st"=>st,
                        "se"=>se,
                        "sr"=>sr,
                        "si"=>identifier,
                        "sp"=>permissions,
                        "sip"=>ip,
                        "spr"=>url.scheme,
                        "sig"=>sig,
                        "rscc"=>cache_control,
                        "rscd"=>content_disposition,
                        "rsce"=>content_encoding,
                        "rscl"=>content_language,
                        "rsct"=>content_type
                    ]
    filter!(kv->!isempty(kv[2]), query_params)
    query_string = join(map((kv)->kv[1]*"="*URIs.escapeuri(kv[2]), query_params), "&")
    raw_url * "?" * query_string
end

export SASPermission, appendSAS


text:
module BenchTeeRFFilter

using BenchmarkTools
using Transducers

const SUITE = BenchmarkGroup()

let xs = 1:1000
    rf = TeeRF(Filter(isodd)'(+), Filter(iseven)'(+))

    SUITE["noinit"] = @benchmarkable foldxl($rf, $xs)
    SUITE["withinit"] = @benchmarkable foldxl($rf, $xs; init = (0, 0))
end

end  # module
BenchTeeRFFilter.SUITE


text:
# ------------------------------------------------------------------------------
# linear elasticity example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh3D(1.0, 1.0, 1.0)
finep = 2
coarsep = 1
numbercomponents = 3
dimension = 3
finebasis = TensorH1LagrangeBasis(finep + 1, finep + 1, numbercomponents, dimension)
coarsebasis = TensorH1LagrangeBasis(coarsep + 1, finep + 1, numbercomponents, dimension)
ctofbasis = TensorH1LagrangeBasis(
    coarsep + 1,
    finep + 1,
    numbercomponents,
    dimension,
    collocatedquadrature = true,
)

# constants
E = 1E6                             # Young's modulus
ν = 0.3                             # Poisson's ratio
λ = E * ν / ((1 + ν) * (1 - 2 * ν)) # Lamé parameters
μ = E / (2 * (1 + ν))

function linearelasticityweakform(deltadu::Array{Float64}, w::Array{Float64})
    # strain
    dϵ = (deltadu + deltadu') / 2
    # strain energy
    dσ_11 = (λ + 2μ) * dϵ[1, 1] + λ * dϵ[2, 2] + λ * dϵ[3, 3]
    dσ_22 = λ * dϵ[1, 1] + (λ + 2μ) * dϵ[2, 2] + λ * dϵ[3, 3]
    dσ_33 = λ * dϵ[1, 1] + λ * dϵ[2, 2] + (λ + 2μ) * dϵ[3, 3]
    dσ_12 = μ * dϵ[1, 2]
    dσ_13 = μ * dϵ[1, 3]
    dσ_23 = μ * dϵ[2, 3]
    dσ = [dσ_11 dσ_12 dσ_13; dσ_12 dσ_22 dσ_23; dσ_13 dσ_23 dσ_33]

    # delta dv
    deltadv = dσ * w[1]

    return [deltadv']
end

# linear elasticity operators
function makeoperator(basis::TensorBasis)
    inputs = [
        OperatorField(basis, [EvaluationMode.gradient], "gradent of deformation"),
        OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights"),
    ]
    outputs = [
        OperatorField(
            basis,
            [EvaluationMode.gradient],
            "test function gradient of deformation",
        ),
    ]
    return Operator(linearelasticityweakform, mesh, inputs, outputs)
end
fineoperator = makeoperator(finebasis)
coarseoperator = makeoperator(coarsebasis)

# Chebyshev smoother
chebyshev = Chebyshev(fineoperator)

# p-multigrid preconditioner
multigrid = PMultigrid(fineoperator, coarseoperator, chebyshev, [ctofbasis])

# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------


text:
import BioSimulator: build_output, update_samplepath!, save_state, get_regular_path

function make_sample_path(f, x_initial, tdata, states, save_points)
    xw = build_output(f, nothing, x_initial, nothing)

    for (t, state) in zip(tdata, states)
        update_samplepath!(f, xw, nothing, t, state, nothing, save_points)
    end

    return xw
end

function check_sample_path(xw, ts, states)
    # always save initial state?
    @test length(xw) == length(ts) + 1

    # data should have been copied in the same order
    @assert length(states) == length(ts)

    for (i, state) in enumerate(states)
        @test xw.u[i+1] == state
        @test xw.t[i+1] == ts[i]
    end
end

function save_ones(simulator, state, model)
    return ones(length(state))
end

@testset "SamplePath" begin
    n = 50
    d = 3
    tdata = sort!(10 * rand(n))
    states = [rand(0:n, d) for _ in tdata]

    # initialization
    x_initial = zeros(Int, d)

    @testset "save_points = nothing" begin
        xw = make_sample_path(save_state, x_initial, tdata, states, nothing)
        check_sample_path(xw, tdata, states)
    end

    @testset "save_points = iterable" begin
        stops_array = [0.5, 6.0]
        stops_range = 0.0:0.5:5
        # stops_tuple = (0.5, 6.0)

        # cases = (stops_array, stops_range, stops_tuple)
        cases = (stops_array, stops_range)

        @testset "$(Base.typename(typeof(stops)))" for stops in cases
            ts = filter(!isequal(0), stops)
            idxs = [searchsortedlast(tdata, t)+1 for t in ts]
            xw = make_sample_path(save_state, x_initial, tdata, states, stops)
            check_sample_path(xw, ts, states[idxs])
        end
    end

    @testset "save_function" begin
        xw = make_sample_path(save_ones, x_initial, tdata, states, nothing)
        # not happy with this test
        expected = save_ones.(nothing, states, nothing)
        check_sample_path(xw, tdata, expected)
    end

    @testset "regular paths" begin
        stops = 0.0:0.5:5
        ts = filter(!isequal(0), stops)
        idxs = [searchsortedlast(tdata, t)+1 for t in ts]

        # make SamplePath from full observations
        xw = make_sample_path(save_state, x_initial, tdata, states, nothing)

        # make SamplePath from partial observations
        yw = make_sample_path(save_state, x_initial, tdata, states, stops)

        # transform full SamplePath to 'regular' path
        zw = get_regular_path(xw, stops)

        @test yw == zw
    end

    @testset "simulate" begin
        network = Network("extinction test")
        network <= Species("X", 1)
        network <= Reaction("birth", 1e-6, "X --> X + X")
        network <= Reaction("death", 1e2, "X --> 0")

        # case: save_points = nothing
        # expected behavior: tfinal is included
        trajectory = simulate(network, Direct(), tfinal = 10.0)

        @test trajectory.t[end] == 10.0
    end
end


text:
"""
    reorder_policy(env::SupplyChainEnv, param1::Dict, param2::Dict,
                        policy_type::Union{Dict, Symbol} = :rQ,
                        review_period::Union{Int, StepRange, Vector, Dict} = 1,
                        min_order_qty::Union{Real, Dict} = 0)

Apply an inventory policy to specify the replinishment orders for each material
    throughout the `SupplyChainEnv`.

# Arguments
- `env::SupplyChainEnv`: inventory management environment
- `param1::Dict`: the `s` or `r` parameter in each node for each material in the system. The `keys` are of the form `(node, material)`.
- `param2::Dict`: the `S` or `Q` parameter in each node for each material in the system. The `keys` are of the form `(node, material)`.
- `policy_type::Union{Symbol, Dict}`: `:rQ` for an `(r,Q)` policy, or `:sS` for an `(s,S)` policy. If passing a `Dict`, the policy type should be specified for each node (keys).
- `review_period::Union{Int, StepRange, Vector, Dict}`: number of periods between each inventory review (Default = `1` for continuous review.). If a `StepRange` or `Vector` is used, the `review_period` indicates which periods the review is performed on. If a `Dict` is used, the review period should be specified for each `(node, material)` `Tuple` (keys). The values of this `Dict` can be either `Int`, `StepRange`, or `Vector`. Any missing `(node, material)` key will be assigned a default value of 1.
- `min_order_qty::Union{Real, Dict}`: minimum order quantity (MOQ) at each supply node. If a `Dict` is passed, the MOQ should be specified for each `(node, material)` `Tuple` (keys). The values should be `Real`. Any missing key will be assigned a default value of 0.
"""
function reorder_policy(env::SupplyChainEnv, param1::Dict, param2::Dict;
                        policy_variable::Union{Dict,Symbol} = :inv_position,
                        policy_type::Union{Dict, Symbol} = :rQ,
                        review_period::Union{Int, StepRange, Vector, Dict} = 1,
                        min_order_qty::Union{Real, Dict} = 0)

    #check review period
    null_action = zeros(length(env.materials)*ne(env.network))
    if review_period isa Int && !iszero(mod(env.period,review_period)) #if not in review period, send null action
        return null_action
    elseif review_period isa Union{StepRange,Vector} && !in(env.period, review_period)
        return null_action
    end

    #read parameters
    nodes = [n for n in vertices(env.network) if !isempty(inneighbors(env.network, n))] #all non-source nodes can place orders
    supply_nodes = [n for n in vertices(env.network) if !isempty(outneighbors(env.network, n))] #all non-sink nodes can place orders
    arcs = [(e.src, e.dst) for e in edges(env.network)]
    mats = env.materials

    #check inputs
    for p in mats
        for n in nodes
            for param in [param1, param2]
                if !in((n,p), keys(param))
                    param[(n,p)] = -1 #if no policy given for a node/material, set the params to -1 so that a reorder is never triggered
                end
            end
            if review_period isa Dict && !in((n,p), keys(review_period))
                review_period[(n,p)] = 1 #set to default value of 1
            end
        end
        for n in supply_nodes
            if min_order_qty isa Dict && !in((n,p), keys(min_order_qty))
                min_order_qty[(n,p)] = 0 #set to default value of 0
            end
        end
    end

    #filter data for policy
    state1_df = filter(:period => j -> j == env.period, env.inv_position, view=true) #inventory position
    state1_grp = groupby(state1_df, [:node, :material]) #group by node and material
    state2_df = filter(:period => j -> j == env.period, env.ech_position, view=true) #inventory position
    state2_grp = groupby(state2_df, [:node, :material]) #group by node and material

    #create action matrix
    action = zeros(length(mats), length(arcs)) #initialize action matrix
    for n in nodes, (k, p) in enumerate(mats)
        #if trigger level is negative, skip it
        param1[n,p] < 0 && continue
        #if not a review period, skip
        if review_period isa Dict
            review_period[n,p] isa Int && !iszero(mod(env.period, review_period[n,p])) && continue #trigger if not in review period
            review_period[n,p] isa Union{StepRange,Vector} && !in(env.period, review_period[n,p]) && continue
        end
        #get MOQ
        supplier = get_prop(env.network, n, :supplier_priority)[p][1] #get first priority supplier
        MOQ = min_order_qty isa Real ?
                min_order_qty :
                min_order_qty[(supplier,p)]
        #initialize reorder quantity
        reorder = 0
        #check if reorder is triggered & set reorder policy
        pol_var = policy_variable isa Dict ? policy_variable[n] : policy_variable
        if pol_var == :inv_position
            state = state1_grp[(node = n, material = p)].level[1]
        elseif pol_var == :ech_position
            state = state2_grp[(node = n, material = p)].level[1]
        else
            @assert pol_var in [:inv_position, :ech_position] "Policy variable must be either `:inv_position` or `:ech_position`."
        end
        if state <= param1[n,p]
            pol_type = policy_type isa Dict ? policy_type[n] : policy_type
            if pol_type == :rQ #rQ policy
                reorder = max(param2[n,p], MOQ)
            elseif pol_type == :sS #sS policy
                reorder = state < param2[n,p] ? max(param2[n,p] - state, MOQ) : 0
            else
                @assert pol_type in [:rQ, :sS] "The policy type must be either `:rQ` or `:sS`."
            end
        else
            continue
        end
        #assign reorder quantity
        j = findfirst(i -> i == (supplier, n), arcs) #find index in action matrix
        action[k, j] = reorder #sate request quantity
    end

    return collect(Iterators.flatten(action))
end

"""
    simulate_policy!(env::SupplyChainEnv, args...)

Step through a simulation using a specified reorder policy. `args` are the
arguments that are passed to the `reorder_policy` function.
"""
function simulate_policy!(env::SupplyChainEnv, args...; kwargs...)
    for t in 1:env.num_periods
        action = reorder_policy(env, args...; kwargs...)
        (env)(action)
    end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to TLE parser.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Macros tle_str and tlenc_str
# ============================

@testset "Macros tle_str and tlenc_str" begin
    # Read the SCDs TLE from the file.
    tles_file = read_tle("./SCDs.tle")

    # Read the same TLEs from a string.
    tles_str = tle"""
        SCD 1
        1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
        2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066

        SCD 2
        1 25504U 98060A   18165.15074951  .00000201  00000-0  55356-5 0  9994
        2 25504  24.9961  80.1303 0017060 224.4822 286.6438 14.44043397 37312
        """

    # Read the same TLES from a string with wrong checksums.
    # This should not output any exceptions.
    tles_str_nc = tlenc"""
        SCD 1
        1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9990
        2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338060

        SCD 2
        1 25504U 98060A   18165.15074951  .00000201  00000-0  55356-5 0  9990
        2 25504  24.9961  80.1303 0017060 224.4822 286.6438 14.44043397 37310
        """

    # Compare the TLEs.
    @test length(tles_file) == length(tles_str)
    @test length(tles_file) == length(tles_str_nc)

    for i = 1:length(tles_file)
        for sym in fieldnames(TLE)
            @test getfield(tles_file[i], sym) == getfield(tles_str[i], sym)

            # Skip the comparison of checksums for `tles_str_nc`.
            ( (sym == :checksum_l1) || (sym == :checksum_l2) ) && continue
            @test getfield(tles_file[i], sym) == getfield(tles_str_nc[i], sym)
        end
    end
end

# Function: read_tle_from_string
# ==============================

@testset "Function read_tle_from_string" begin
    # Read the SCD 1 TLE from the file.
    tle_scd1_expected = read_tle("./SCDs.tle")[1]

    # Parse the same TLE using the function `read_tle_from_string`.
    tle_scd1_result = read_tle_from_string(
        "1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991",
        "2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066")[1]

    # Compare the TLEs.
    for sym in fieldnames(TLE)
        # Skip name and satellite number.
        ( (sym == :name) || (sym == :sat_num) ) && continue

        @test getfield(tle_scd1_result, sym) == getfield(tle_scd1_expected, sym)
    end
end

# Date conversion
# ===============

@testset "Epoch year conversion" begin
    tles = tlenc"""
    1 22490U 93009B   75165.62596833  .00000225  00000-0  11410-4 0  9990
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338060

    1 22490U 93009B   76165.62596833  .00000225  00000-0  11410-4 0  9990
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338060
    """

    date_1 = JDtoDate(tles[1].epoch)
    date_2 = JDtoDate(tles[2].epoch)

    @test date_1[1] == 2075
    @test date_2[1] == 1976
end

# Test errors
# ===========

@testset "TLE format errors" begin

    @test_throws ErrorException read_tle_from_string("""
    SCD 1
    1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066
    3 123
    """)

    # Errors in line 1
    # --------------------------------------------------------------------------

    @test_throws ErrorException read_tle_from_string("""
    SCD 1
    1 022490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066
    """)

    @test_throws ErrorException read_tle_from_string("""
    SCD 1
    1 022490U 93009B   181656.2596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066
    """)

    # TODO: Errors in line 2.
end

@testset "Checksum errors" begin
    @test_throws ErrorException read_tle("./SCDs-wrong_checksum.tle")
    @test_throws LoadError @eval tle"""
    SCD 1
    1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991
    2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338066

    SCD 2
    1 25504U 98060A   18165.15074951  .00000201  00000-0  55356-5 0  9994
    2 25504  24.9961  80.1303 0017060 224.4822 286.6438 14.44043397 37313
    """
    @test_throws ErrorException read_tle_from_string(
        "1 22490U 93009B   18165.62596833  .00000225  00000-0  11410-4 0  9991",
        "2 22490  24.9690 231.7852 0042844 200.7311 292.7198 14.44524498338063",
        verify_checksum = true)
end


text:
using Test
using POMDPTesting
using FiniteHorizonPOMDPs
using POMDPs
using POMDPModelTools
import POMDPModels: SimpleGridWorld, BabyPOMDP
import POMDPPolicies: FunctionPolicy
import POMDPSimulators: stepthrough

@testset "interface" begin
    @test HorizonLength(SimpleGridWorld()) == InfiniteHorizon()
    @test HorizonLength(BabyPOMDP()) == InfiniteHorizon()
end

@testset "fixhorizon" begin
    include("fixhorizon.jl")
end


text:
"""
    function bunny()

A cute little wabbit

# Example

```jldoctest
julia> cowsay("Look at that tortoise way off in the distance.", cow=Cowsay.bunny)
 ______________________________________
/ Look at that tortoise way off in the \\
\\ distance.                            /
 --------------------------------------
  \\
   \\   \\
        \\ /\\
        ( )
      .( o ).

```
"""
function bunny(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
  $thoughts
   $thoughts   \\
        \\ /\\
        ( )
      .( o ).
"""
return the_cow
end


text:
using Plots
# include(joinpath(@__DIR__, "car_obstacles.jl"))

X̄, Ū = unpack(Z̄, prob)

# Position trajectory
Px = [X̄[t][1] for t = 1:T]
Py = [X̄[t][2] for t = 1:T]

pts = Plots.partialcircle(0.0, 2.0 * π, 100, 0.1)
cx, cy = Plots.unzip(pts)
cx1 = [_cx + circles[1][1] for _cx in cx]
cy1 = [_cy + circles[1][2] for _cy in cy]
cx2 = [_cx + circles[2][1] for _cx in cx]
cy2 = [_cy + circles[2][2] for _cy in cy]
cx3 = [_cx + circles[3][1] for _cx in cx]
cy3 = [_cy + circles[3][2] for _cy in cy]
cx4 = [_cx + circles[4][1] for _cx in cx]
cy4 = [_cy + circles[4][2] for _cy in cy]

plt = plot(Shape(cx1, cy1), color = :red, label = "", linecolor = :red)
plt = plot!(Shape(cx2, cy2), color = :red, label = "", linecolor = :red)
plt = plot!(Shape(cx3, cy3), color = :red, label = "", linecolor = :red)
plt = plot!(Shape(cx4, cy4), color = :red, label = "", linecolor = :red)
plt = plot!(Px, Py, aspect_ratio = :equal, xlabel = "x", ylabel = "y",
    width = 4.0, label = "TO", color = :purple, legend = :topleft)

#DPO
X̄_dpo, Ū_dpo = unpack(Z[prob_dpo.prob.idx.nom], prob_dpo.prob.prob.nom)
X̄_dpo_mean, Ū_dpo_mean = unpack(Z[prob_dpo.prob.idx.mean], prob_dpo.prob.prob.mean)

X̄[end]
X̄_dpo[end]
X̄_dpo_mean[end]
X_sample = []
U_sample = []
for i = 1:N
	x, u = unpack(Z[prob_dpo.prob.idx.sample[i]],
		prob_dpo.prob.prob.sample[i])
	push!(X_sample, x)
	push!(U_sample, u)
end

# Position trajectory
Px_dpo = [X̄_dpo[t][1] for t = 1:T]
Py_dpo = [X̄_dpo[t][2] for t = 1:T]

Px_dpo_mean = [X̄_dpo_mean[t][1] for t = 1:T]
Py_dpo_mean = [X̄_dpo_mean[t][2] for t = 1:T]
plt = plot!(Px_dpo_mean, Py_dpo_mean, width = 2.0, label = "DPO (mean)", color = :green)

Px_sample = []
Py_sample = []
for i = 1:N
	px = [X_sample[i][t][1] for t = 1:T]
	py = [X_sample[i][t][2] for t = 1:T]
	push!(Px_sample, px)
	push!(Py_sample, py)
end

plt = plot!(Px_dpo, Py_dpo, width = 4.0, label = "DPO (nominal)", color = :orange)
for i = 1:N
	plt = plot!(Px_sample[i], Py_sample[i], width = 1.0, label = "", color = :grey)
end
display(plt)

# using PGFPlots
# const PGF = PGFPlots
#
# # TO trajectory
# p_nom = PGF.Plots.Linear(hcat(X̄...)[1,:], hcat(X̄...)[2,:],
#     mark="", style="color=cyan, line width=3pt, solid", legendentry = "TO")
#
# # DPO trajectory
# p_dpo = PGF.Plots.Linear(hcat(X̄_dpo...)[1,:], hcat(X̄_dpo...)[2,:],
#     mark="", style = "color=orange, line width=3pt, solid",legendentry = "DPO")
#
# # Obstacles
# p_circle = [PGF.Plots.Circle(circle..., style = "color=black,fill=black") for circle in circles]
#
# a = Axis([p_circle;
#     p_nom;
#     p_dpo
#     ],
#     xmin = -0.4, ymin = -0.1, xmax = 1.4, ymax = 1.1,
#     axisEqualImage = true,
#     hideAxis = false,
# 	ylabel = "y",
# 	xlabel = "x",
# 	legendStyle = "{at={(0.01,0.99)},anchor=north west}",
# 	)
#
# # Save to tikz format
# dir = joinpath(@__DIR__, "results")
# PGF.save(joinpath(dir, "car_obstacles.tikz"), a, include_preamble = false)

# # Animation
# include(joinpath(pwd(),"dynamics/visualize.jl"))
#
# vis = Visualizer()
# open(vis)
# visualize!(vis,model,[X_nom...,[X_nom[end] for t = 1:T]...],Δt=H_nom_sample[1])
# # settransform!(vis["/Cameras/default"], compose(LinearMap(RotY(-pi/2.0)*RotX(pi/2)),Translation(-1, 0,0)))
#
# for i = 1:4
#     cyl = Cylinder(Point3f0(xc[i],yc[i],0),Point3f0(xc[i],yc[i],0.1),convert(Float32,0.035))
#     setobject!(vis["cyl$i"],cyl,MeshPhongMaterial(color=RGBA(1,0,0,1.0)))
# end
#
# q_to = deepcopy(X_nom)
# for t = 1:T
# 	setobject!(vis["traj_to$t"], Sphere(Point3f0(0),
# 		convert(Float32,0.075)),
# 		MeshPhongMaterial(color=RGBA(0.0,255.0/255.0,255.0/255.0,0.75)))
# 	settransform!(vis["traj_to$t"], Translation((q_to[t][1],q_to[t][2],-0.1)))
# 	setvisible!(vis["traj_to$t"],true)
# end
#
# q_dpo = deepcopy(X_nom_sample)
# for t = 1:T
# 	setobject!(vis["traj_dpo$t"], Sphere(Point3f0(0),
# 		convert(Float32,0.075)),
# 		MeshPhongMaterial(color=RGBA(255.0/255.0,127.0/255.0,0.0,0.75)))
# 	settransform!(vis["traj_dpo$t"], Translation((q_dpo[t][1],q_dpo[t][2],-0.05)))
# 	setvisible!(vis["traj_dpo$t"],false)
# end
#
# q = X_nom_sample
# obj_path = joinpath(pwd(),"/home/taylor/Research/direct_policy_optimization/dynamics/cybertruck/cybertruck.obj")
# mtl_path = joinpath(pwd(),"/home/taylor/Research/direct_policy_optimization/dynamics/cybertruck/cybertruck.mtl")
#
# ctm = ModifiedMeshFileObject(obj_path,mtl_path,scale=0.05)
# t = 1
# setobject!(vis["ct1"],ctm)
# settransform!(vis["ct1"], compose(Translation([q[t][1];q[t][2];0.0]),LinearMap(RotZ(q[t][3]+pi)*RotX(pi/2.0))))
# setvisible!(vis["ct1"],false)
# t = 1
# setobject!(vis["ct2"],ctm)
# settransform!(vis["ct2"], compose(Translation([q[t][1];q[t][2];0.0]),LinearMap(RotZ(q[t][3]+pi)*RotX(pi/2.0))))
# setvisible!(vis["ct2"],false)
# t = 20#17#20
# setobject!(vis["ct3"],ctm)
# settransform!(vis["ct3"], compose(Translation([q[t][1];q[t][2];0.0]),LinearMap(RotZ(q[t][3]+pi)*RotX(pi/2.0))))
# setvisible!(vis["ct3"],false)
#
# t = 34#36#34
# setobject!(vis["ct4"],ctm)
# settransform!(vis["ct4"], compose(Translation([q[t][1];q[t][2];0.0]),LinearMap(RotZ(q[t][3]+pi)*RotX(pi/2.0))))
# setvisible!(vis["ct4"],false)
# t = T
# setobject!(vis["ct5"],ctm)
# settransform!(vis["ct5"], compose(Translation([q[t][1];q[t][2];0.0]),LinearMap(RotZ(q[t][3]+pi)*RotX(pi/2.0))))
# setvisible!(vis["ct5"],false)
#
# # setvisible!(vis["/Background"], true)
# # setprop!(vis["/Background"], "top_color", [135,206,250])
# # setprop!(vis["/Background"], "bottom_color", [135,206,250])
# # vis["/Background"].core

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
export LazyObservable, lazy, set!

import Base: show

# TODO: Untested and undocumented

struct LazyObservable{D} <: Subscribable{D}
    pending :: PendingSubjectInstance{ Any, Subject{Any,AsapScheduler,AsapScheduler} }
end

function LazyObservable(::Type{T}, pending) where T
    return LazyObservable{T}(pending)
end

set!(lazy::LazyObservable, observable::S) where S = on_lazy_set!(lazy, as_subscribable(S), observable)

on_lazy_set!(lazy::LazyObservable{D},  ::InvalidSubscribableTrait,       observable) where D                = throw(InvalidSubscribableTraitUsageError(observable))
on_lazy_set!(lazy::LazyObservable{D1}, ::SimpleSubscribableTrait{D2},    observable) where { D1, D2 <: D1 } = _on_lazy_set!(lazy, observable)
on_lazy_set!(lazy::LazyObservable{D1}, ::ScheduledSubscribableTrait{D2}, observable) where { D1, D2 <: D1 } = _on_lazy_set!(lazy, observable)

function _on_lazy_set!(lazy, observable)
    next!(lazy.pending, observable)
    complete!(lazy.pending)
    return nothing
end

function on_subscribe!(observable::LazyObservable{D}, actor) where D
    return subscribe!(observable.pending |> switch_map(D), actor)
end

lazy(::Type{T} = Any) where T = LazyObservable(T, PendingSubject(Any))

Base.show(io::IO, ::LazyObservable{D}) where D = print(io, "LazyObservable($D)")