text:
"""
Initializes the `TabuSearchOptimizer` and the `SearchState`.
"""
function initialize_search(
    rr2::RR2,
    cons::Vector{Constraint},
    pen::Vector{Int},
    start::Matrix{Int},
    maxiter::Int,
    tabulength::Int,
    infeaslength::Int,
)
    entries = conentries(rr2, cons)
    tso = TabuSearchOptimizer(
        maxiter,
        tabulength,
        infeaslength,
        rr2,
        cons,
        entries,
    )
    hpen, spen = 0, 0
    offset = [eval_constraint(start, c) for c in cons]
    for (i, c) in enumerate(cons)
        if c.pen == 0
            hpen += offset[i] * pen[i]
        else
            spen += offset[i] * pen[i]
        end
    end
    moves = getmoves(rr2)

    tabu = Dict((m, 0) for m in moves)
    blocked = Dict((m, 0) for m in getmoves(rr2))
    optval = hpen == 0 ? spen : typemax(Int)
    state = SearchState(
        0, #it
        hpen,
        spen,
        copy(start), #sched
        optval,
        copy(start), #optsched
        0, #lastimprov
        offset,
        moves,
        tabu,
        blocked,
        pen,
    )
    return tso, state
end

"""
Initialize buffers used in the Tabu Search. `offsetbuf` is used for storing
constraint offsets, `schedbuf` is used as a schedule buffer and `difbuf` is used
for move evaluation.
"""
function init_buffers(tso::TabuSearchOptimizer, state::SearchState)
    offsetbuf = copy(state.offset)
    schedbuf = similar(state.sched)
    difbuf = BitMatrix(undef, tso.rr2.nteams, tso.rr2.nslots)
    return offsetbuf, schedbuf, difbuf
end

"""
    search_step!(
        tso::TabuSearchOptimizer,
        state::SearchState,
        offsetbuf::Vector{Int},
        schedbuf::Matrix{Int},
        difbuf::BitMatrix,
    )

Selects a suitable neighbor solution and moves to it.
"""
function search_step!(
    tso::TabuSearchOptimizer,
    state::SearchState,
    offsetbuf::Vector{Int},
    schedbuf::Matrix{Int},
    difbuf::BitMatrix;
    accept_infeas = false,
)
    state.it += 1
    opt_feas_mv, opt_infeas_mv = nothing, nothing
    opt_infeas_hpen, opt_infeas_spen = typemax(Int), typemax(Int)
    opt_feas_spen = typemax(Int)

    # compute constrained entries of violated constraints
    infeas = falses(size(schedbuf))
    for idx in 1:length(tso.constraints)
        if state.offset[idx] > 0
            infeas = infeas .| @view tso.conentries[:,:,idx]
        end
    end
    if !any(infeas)
        infeas .= true
    end

    # randomize move order
    shuffle!(state.moves)

    for mv in state.moves
        # identical move already evaluated this iteration
        state.blocked[mv] == state.it && continue
        # reset buffers
        offsetbuf .= state.offset
        schedbuf .= state.sched

        move!(schedbuf, mv)
        # check if move changes constrained entry
        difbuf .= (schedbuf .!= state.sched)
        !any(difbuf .& infeas) && continue

        # evaluate move and block equivalent moves
        eval_move!(
            offsetbuf,
            schedbuf,
            difbuf,
            state.sched,
            tso.constraints,
            tso.conentries,
        )
        block_equiv!(state.blocked, mv, difbuf, state.it)

        hpen, spen = 0, 0
        hpenold, spenold = 0, 0
        for i in 1:length(tso.constraints)
            if tso.constraints[i].pen == 0
                hpen += state.pen[i] * offsetbuf[i]
                hpenold += state.pen[i] * state.offset[i]
            else
                spen += state.pen[i] * offsetbuf[i]
                spenold += state.pen[i] * state.offset[i]
            end
        end

        if state.tabu[mv] >= state.it # move is marked tabu
            # aspiration
            if hpen == 0 && spen < state.optval
                opt_feas_mv, opt_feas_spen = mv, spen
                break
            end
        else
            if hpen == 0
                if spen < spenold
                    opt_feas_mv, opt_feas_spen = mv, spen
                    break
                elseif spen < opt_feas_spen
                    opt_feas_mv, opt_feas_spen = mv, spen
                end
            else
                if (spen < spenold && hpen <= hpenold) || (spen <= spenold && hpen < hpenold)
                    opt_infeas_mv = mv
                    opt_infeas_hpen, opt_infeas_spen = hpen, spen
                    break
                elseif hpenold == 0 && (spen, hpen) < (opt_infeas_spen, opt_infeas_hpen)
                    opt_infeas_mv = mv
                    opt_infeas_hpen, opt_infeas_spen = hpen, spen
                elseif hpenold > 0 && (hpen, spen) < (opt_infeas_hpen, opt_infeas_spen)
                    opt_infeas_mv = mv
                    opt_infeas_hpen, opt_infeas_spen = hpen, spen
                end
            end
        end
    end

    if opt_feas_mv === nothing # no feasible neighbor solution found
        @assert opt_infeas_mv !== nothing
        bestmv = opt_infeas_mv
        state.hpen, state.spen = opt_infeas_hpen, opt_infeas_spen
    else
        if accept_infeas && opt_infeas_spen < opt_feas_spen
            @assert opt_infeas_mv !== nothing
            bestmv = opt_infeas_mv
            state.hpen, state.spen = opt_infeas_hpen, opt_infeas_spen
        else
            bestmv = opt_feas_mv
            state.hpen, state.spen = 0, opt_feas_spen
        end
    end
    schedbuf .= state.sched
    move!(schedbuf, bestmv)
    difbuf .= (schedbuf .!= state.sched)
    eval_move!(
        state.offset,
        schedbuf,
        difbuf,
        state.sched,
        tso.constraints,
        tso.conentries,
    )
    move!(state.sched, bestmv)

    if state.hpen == 0 && state.spen < state.optval
        state.optsched .= state.sched
        state.optval = state.spen
        state.lastimprov = state.it
    end

    # mark bestmv and equivalent moves as tabu
    block_equiv!(state.tabu, bestmv, difbuf, state.it + tso.tabulength)
    return bestmv
end

"""
    eval_move!(
        offset::Vector{Int},
        schedbuf::Matrix{Int},
        difbuf::BitMatrix,
        sched::Matrix{Int},
        constraints::Vector{Constraint},
        entries::BitArray{3},
    )

evaluates a move and updates constraint offsets. Only iterates through
constraints that are constraining (team, slot) entries changed by move.

# Arguments

`schedbuf::Matrix{Int}` - schedule buffer
`difbuf::BitMatrix` - schedule difference buffer
`sched::Matrix{Int}` - current solution
`constraints::Vector{Constraint}` - constraints of the model
`entries::BitArray{3}` - used for determining relevant constraints
"""
function eval_move!(
    offsetbuf::Vector{Int},
    schedbuf::Matrix{Int},
    difbuf::BitMatrix,
    sched::Matrix{Int},
    constraints::Vector{Constraint},
    entries::BitArray{3},
)
    for (i, c) in enumerate(constraints)
        if any(difbuf .& entries[:, :, i])
            offsetbuf[i] = eval_constraint(schedbuf, c)
        end
    end
end

"""
Marks equivalent moves to mv as blocked.
"""
function block_equiv!(
    blocked::Dict{Move,Int},
    mv::Move,
    dif::BitMatrix,
    it::Int,
)
    blocked[mv] = it
    if mv[1] == :partial_swap_teams!
        t1, t2, s = mv[2:end]
        for slot = 1:size(dif, 2)
            if dif[t1, slot]
                blocked[(:partial_swap_teams!, t1, t2, slot)] = it
            end
        end
    elseif mv[1] == :partial_swap_rounds!
        t, s1, s2 = mv[2:end]
        for team = 1:size(dif, 1)
            if dif[team, s1]
                blocked[(:partial_swap_rounds!, team, s1, s2)] = it
            end
        end
    end
end


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

@recipe function f(pset::PointSet{Dim,T}) where {Dim,T}
  pts = [pset[i] for i in 1:nelements(pset)]

  seriestype --> :scatter
  seriescolor --> :black
  legend --> false

  if Dim == 1
    @series begin
      first.(coordinates.(pts)), fill(0, nelements(pset))
    end
  else
    aspect_ratio --> :equal
    @series begin
      @. Tuple(coordinates(pts))
    end
  end
end

@recipe function f(pset::PointSet{Dim,T}, data::AbstractVector) where {Dim,T}
  pts = [pset[i] for i in 1:nelements(pset)]

  seriestype --> :scatter

  if Dim == 1
    first.(coordinates.(pts)), data
  else
    aspect_ratio --> :equal
    marker_z --> data
    colorbar --> true
    @. Tuple(coordinates(pts))
  end
end

ext: jl


text:
using DelayEmbeddings
using StaticArrays
using Test
import Downloads

# Download some test timeseries
tsfolder = joinpath(@__DIR__, "timeseries")
todownload1 = ["$n.csv" for n in 1:4]
todownload = ["test_time_series_lorenz_standard_N_10000_multivariate.csv", "test_time_series_roessler_N_10000_multivariate.csv"]
append!(todownload, todownload1)
repo = "https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/timeseries"
mkpath(tsfolder)
for a in todownload
    Downloads.download(repo*"/"*a, joinpath(tsfolder, a))
end

ti = time()

diffeq = (atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))

@testset "DelayEmbeddings tests" begin
    include("dataset_tests.jl")
    include("embedding_tests.jl")
    include("traditional/delaytime_test.jl")
    include("traditional/embedding_dimension_test.jl")
    include("unified/test_pecora.jl")
    include("unified/uzal_cost_test.jl")
    include("unified/test_pecuzal_embedding.jl")
    include("unified/test_garcia.jl")
    include("unified/mdop_tests.jl")

end

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, sigdigits=3), " seconds or ", round(ti/60, sigdigits=3), " minutes")


text:
using Compat.Test

@testset "Test ApplicationBuilder (by compiling examples/*.jl)" begin
    include("ApplicationBuilder.jl")
end
@testset "Command-line interface (compiling examples/*.jl)" begin
    include("build_app-cli.jl")
end


text:
function TroposphericGradientConverter(arg0::SpacecraftState, arg1::jint, arg2::AttitudeProvider)
    return TroposphericGradientConverter((SpacecraftState, jint, AttitudeProvider), arg0, arg1, arg2)
end

function get_free_state_parameters(obj::TroposphericGradientConverter)
    return jcall(obj, "getFreeStateParameters", jint, ())
end

function get_parameters(obj::TroposphericGradientConverter, arg0::FieldSpacecraftState, arg1::DiscreteTroposphericModel)
    return jcall(obj, "getParameters", Vector{Gradient}, (FieldSpacecraftState, DiscreteTroposphericModel), arg0, arg1)
end

function get_state(obj::TroposphericGradientConverter, arg0::DiscreteTroposphericModel)
    return jcall(obj, "getState", FieldSpacecraftState, (DiscreteTroposphericModel,), arg0)
end


text:
abstract type Generator <: StaticInjection end
const Generators = Array{<:Generator, 1}

abstract type HydroGen <: Generator end
abstract type RenewableGen <: Generator end
abstract type ThermalGen <: Generator end

function IS.get_limits(
    valid_range::Union{NamedTuple{(:min, :max)}, NamedTuple{(:max, :min)}},
    unused::T,
) where {T <: Generator}
    # Gets min and max value defined for a field,
    # e.g. "valid_range": {"min":-1.571, "max":1.571}.
    return (min = valid_range.min, max = valid_range.max, zero = 0.0)
end

get_rating(gen::Generator) = get_rating(get_tech(gen))
get_active_power_limits_min(gen::Generator) = get_activepowerlimits(get_tech(gen)).min
get_active_power_limits_max(gen::Generator) = get_activepowerlimits(get_tech(gen)).max


text:
@coeff function K_OH(reactor::Reactor)
  sqrt(
    K_omega * K_I * D^2 /
    ( K_T^1.5 * sqrt(2 * mu_0 * sigma_max_hat) )
  )
end


text:
# Some time-steppers for AIBECS simulations

# TODO:
# - increment other methods
#     - IMEX
#     - other bashford something methods
# - Try to integrate into DiffEq implementation instead

"""
    euler_forward_step(x, p, δt, F)

Returns the Euler-forward increment `x + δt * F(x,p)` (out of place).
"""
euler_forward_step(x, p, δt, F) = x + δt * F(x, p)
"""
    euler_forward_step!(x, p, δt, F)

In place version of [`euler_forward_step`](@ref).
"""
euler_forward_step!(x, p, δt, F) = x .= euler_forward_step(x, p, δt, F)

"""
    euler_backward_step(x, p, δt, F, ∇ₓF)

Returns the Euler-backward increment (out of place).
Solves `H(y) = 0` where `H(y) = F(y,p) + (x - y) / δt`.
"""
function euler_backward_step(x, p, δt, F, ∇ₓF)
    H(y,p) = F(y,p) - y / δt + x / δt
    ∇ₓH(y,p) = ∇ₓF(y,p) - I / δt
    prob = SteadyStateProblem(H, ∇ₓH, x, p)
    return solve(prob, CTKAlg()).u
end
"""
    euler_backward_step!(x, p, δt, F, ∇ₓF)

In place version of [`euler_backward_step`](@ref).
"""
function euler_backward_step!(x, p, δt, F, ∇ₓF)
    H(y,p) = F(y,p) - y / δt + x / δt
    ∇ₓH(y,p) = ∇ₓF(y,p) - I / δt
    prob = SteadyStateProblem(H, ∇ₓH, x, p)
    x .= solve(prob, CTKAlg()).u
    return x
end

"""
    crank_nicolson_step(x, p, δt, F, ∇ₓF)

Returns the crank nicolson increment (out of place).
Solves `H(y) = 0` where `H(y) = (F(y,p) + F(x,p)) / 2 + (x - y) / δt`.
"""
function crank_nicolson_step(x, p, δt, F, ∇ₓF)
    H(y,p) = F(y,p) / 2 - y / δt + F(x,p) / 2 + x / δt
    ∇ₓH(y,p) = ∇ₓF(y,p) / 2 - I / δt
    prob = SteadyStateProblem(H, ∇ₓH, x, p)
    return solve(prob, CTKAlg()).u
end
"""
    crank_nicolson_step!(x, p, δt, F, ∇ₓF)

In place version of [`crank_nicolson_step`](@ref).

For a tracer equation of the form

`∂x/∂t = F(x,p)`

The Crank-Nicolson scheme is given by

`(xᵢ₊₁ - xᵢ) / δt = (F(xᵢ₊₁,p) + F(xᵢ,p)) / 2`

which requires solving `H(y) = 0` where `H` is defined by

`H(y) = (F(y,p) + F(xᵢ,p)) / 2 + (xᵢ - y) / δt`

where `xᵢ` is prescribed.
`H(y) = 0` is thus defined as a steady-state problem and solved via Newton's method.
"""
function crank_nicolson_step!(x, p, δt, F, ∇ₓF)
    H(y,p) = F(y,p) / 2 - y / δt + F(x,p) / 2 + x / δt
    ∇ₓH(y,p) = ∇ₓF(y,p) / 2 - I / δt
    prob = SteadyStateProblem(H, ∇ₓH, x, p)
    x .= solve(prob, CTKAlg()).u
    return x
end

"""
    crank_nicolson_leapfrog_step(xᵢ, xᵢ₋₁, p, δt, T, ∇ₓL, NL)

Returns the Crank-Nicolson-leapfrog-step increment (out of place).
See [`crank_nicolson_leapfrog_step_A⁺_and_A⁻`](@ref) for more details and a more efficient
computation.
"""
function crank_nicolson_leapfrog_step(xᵢ, xᵢ₋₁, p, δt, T, ∇ₓL, NL)
    return (I / 2δt + (T(p) - ∇ₓL(p)) / 2) \ (NL(xᵢ,p) + (I / 2δt - (T(p) - ∇ₓL(p)) / 2) * xᵢ₋₁)
end
"""
    crank_nicolson_leapfrog_step(xᵢ, xᵢ₋₁, p, A⁺::Factorization, A⁻, G)

Returns the Crank-Nicolson-leapfrog-step increment (out of place) using stored factors.
See [`crank_nicolson_leapfrog_step_A⁺_and_A⁻`](@ref) for more information.
"""
function crank_nicolson_leapfrog_step(xᵢ, xᵢ₋₁, p, A⁺::Factorization, A⁻, NL)
    return A⁺ \ (NL(xᵢ,p) + A⁻ * xᵢ₋₁)
end
"""
    crank_nicolson_leapfrog_step_A⁺_and_A⁻(p, δt, T, ∇ₓL)

Returns `A⁺` and `A⁻` for the Crank-Nicolson-leapfrog-step increment (out of place).
`A⁺` is a `Factorization` object and `A⁻` is a matrix.

For a tracer equation of the form

`∂x/∂t = F(x,p) = NL(x,p) + [∇ₓL(p) - T(p)] x`

where `F(x,p)` has been split into 3 terms:
- the linear transport, `T(p) x`,
- a non linear source/sink term, `NL(x,p)`,
- and a linear source/sink term, `L(x,p) = ∇ₓL(p) x` (note the signs of `∇ₓL` and `T`).
The Crank-Nicolson-leapfrog scheme uses this nonlinear/linear split and is given by

`(xᵢ₊₁ - xᵢ₋₁) / 2δt = NL(xᵢ,p) + [∇ₓL(p) - T(p)] (xᵢ₊₁ + xᵢ₋₁) / 2`

or, with the `xᵢ₊₁` terms on the left,

`A⁺ xᵢ₊₁ = NL(xᵢ,p) + A⁻ xᵢ₋₁`

with `2A⁺ = I / δt + T(p) - ∇ₓL(p)` and `2A⁻ = I / δt - T(p) + ∇ₓL(p)`,
which only depend on `p` (not `x`).
Thus, the Crank-Nicolson-leapfrog steps only needs to factorize `A⁺` once, making it much faster than the standard Crank-Nicolson scheme, justifying this function.
"""
function crank_nicolson_leapfrog_step_A⁺_and_A⁻(p, δt, T, ∇ₓL)
    return factorize(I / 2δt + (T(p) - ∇ₓL(p)) / 2), I / 2δt - (T(p) - ∇ₓL(p)) / 2
end


text:
module NDMCIntegrate

export SMC

    using Statistics

    function SMC(; f::Function, R₀::Real = 0.0, R₁::Real, θ₀::Real = 0.0, θ₁::Real = π, Ŋ::Integer = 10^5)
        ΔR = R₁ - R₀
        Δθ = θ₁ - θ₀
        Đist = [
            muladd.(rand(Ŋ), ΔR, R₀),
            muladd.(rand(Ŋ), Δθ, θ₀)]
        fₓ = [f(getindex.(Đist, n)...) for n ∈ 1:Ŋ]
        Avg = mean(fₓ)
        ∫fₓdx = Avg * ΔR * Δθ
        Δ = std(fₓ) / √Ŋ
        return ∫fₓdx, Δ
    end

end


text:
"""
    Solver Data
"""
struct SolverData{T}
    objective::Vector{T}                # objective value
    gradient::Vector{T}                 # Lagrangian gradient
    max_violation::Vector{T}            # maximum constraint violation

    indices_state::Vector{Vector{Int}}  # indices for state trajectory
    indices_action::Vector{Vector{Int}} # indices for control trajectory

    step_size::Vector{T}                # step length
    status::Vector{Bool}                # solver status

    iterations::Vector{Int}

    cache::Dict{Symbol,Vector{T}}       # solver stats
end

function solver_data(dynamics::Vector{Dynamics{T}};
    max_cache=1000) where T

    # indices x and u
    indices_state = Vector{Int}[]
    indices_action = Vector{Int}[]
    n_sum = 0
    m_sum = 0
    n_total = sum([d.num_state for d in dynamics]) + dynamics[end].num_next_state
    for d in dynamics
        push!(indices_state, collect(n_sum .+ (1:d.num_state)))
        push!(indices_action, collect(n_total + m_sum .+ (1:d.num_action)))
        n_sum += d.num_state
        m_sum += d.num_action
    end
    push!(indices_state, collect(n_sum .+ (1:dynamics[end].num_next_state)))

    objective = [Inf]
    max_violation = [0.0]
    step_size = [1.0]
    gradient = zeros(num_trajectory(dynamics))
    cache = Dict(:objective     => zeros(max_cache),
                 :gradient      => zeros(max_cache),
                 :max_violation => zeros(max_cache),
                 :step_size     => zeros(max_cache))

    SolverData(objective, gradient, max_violation, indices_state, indices_action, step_size, [false], [0], cache)
end

function reset!(data::SolverData)
    fill!(data.objective, 0.0)
    fill!(data.gradient, 0.0)
    fill!(data.max_violation, 0.0)
    fill!(data.cache[:objective], 0.0)
    fill!(data.cache[:gradient], 0.0)
    fill!(data.cache[:max_violation], 0.0)
    fill!(data.cache[:step_size], 0.0)
    data.status[1] = false
    data.iterations[1] = 0
end

# TODO: fix iter
function cache!(data::SolverData)
    iter = 1 #data.cache[:iter]
    # (iter > length(data[:objective])) && (@warn "solver data cache exceeded")
    data.cache[:objective][iter] = data.objective[1]
    data.cache[:gradient][iter] = data.gradient
    data.cache[:max_violation][iter] = data.max_violation
    data.cache[:step_size][iter] = data.step_size
    return nothing
end


text:
using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions,
      RecursiveArrayTools, Distributions, Test

println("One parameter case")
f1 = @ode_def begin
  dx = a*x - x*y
  dy = -3y + x*y
end a
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=50))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(1.5,0.1),1.0,1.8)]

bayesian_result = stan_inference(prob1,t,data,priors;num_samples=300,
                                 num_warmups=500,likelihood=Normal)

@test mean(get(bayesian_result.chains,:theta_1)[1]) ≈ 1.5 atol=3e-1

# Test norecompile
bayesian_result2 = stan_inference(prob1,t,data,priors,bayesian_result.model;
                                  num_samples=300,num_warmups=500,likelihood=Normal)

@test mean(get(bayesian_result2.chains,:theta_1)[1]) ≈ 1.5 atol=3e-1

priors = [truncated(Normal(1.,0.01),0.5,2.0),truncated(Normal(1.,0.01),0.5,2.0),truncated(Normal(1.5,0.01),1.0,2.0)]
bayesian_result = stan_inference(prob1,t,data,priors;num_samples=300,
                                 num_warmups=500,likelihood=Normal,sample_u0=true)

@test mean(get(bayesian_result.chains,:theta_1)[1]) ≈ 1. atol=3e-1
@test mean(get(bayesian_result.chains,:theta_2)[1]) ≈ 1. atol=3e-1
@test mean(get(bayesian_result.chains,:theta_3)[1]) ≈ 1.5 atol=3e-1

sol = solve(prob1,Tsit5(),save_idxs=[1])
randomized = VectorOfArray([(sol(t[i]) + .01 * randn(1)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(1.5,0.1),0.5,2)]
bayesian_result = stan_inference(prob1,t,data,priors;num_samples=300,
                                 num_warmups=500,likelihood=Normal,save_idxs=[1])

@test mean(get(bayesian_result.chains,:theta_1)[1]) ≈ 1.5 atol=3e-1

priors = [truncated(Normal(1.,0.01),0.5,2),truncated(Normal(1.5,0.01),0.5,2)]
bayesian_result = stan_inference(prob1,t,data,priors;num_samples=300,
                                 num_warmups=500,likelihood=Normal,save_idxs=[1],sample_u0=true)

@test mean(get(bayesian_result.chains,:theta_1)[1]) ≈ 1. atol=3e-1
@test mean(get(bayesian_result.chains,:theta_2)[1]) ≈ 1.5 atol=3e-1

println("Four parameter case")
f1 = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=50))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(1.5,0.01),0.5,2),truncated(Normal(1.0,0.01),0.5,1.5),
          truncated(Normal(3.0,0.01),0.5,4),truncated(Normal(1.0,0.01),0.5,2)]

bayesian_result = stan_inference(prob1,t,data,priors;num_samples=100,num_warmups=500,vars =(DiffEqBayes.StanODEData(),InverseGamma(4,1)))
@test mean(get(bayesian_result.chains,:theta_1)[1]) ≈ 1.5 atol=1e-1
@test mean(get(bayesian_result.chains,:theta_2)[1]) ≈ 1.0 atol=1e-1
@test mean(get(bayesian_result.chains,:theta_3)[1]) ≈ 3.0 atol=1e-1
@test mean(get(bayesian_result.chains,:theta_4)[1]) ≈ 1.0 atol=1e-1


text:
"""
    Queue

A simple queue interface based on a doubly linked list.
"""
struct Queue{T <: DLinkedList}
    data::T
end

"""
    enqueue!(q::Queue{DLinkedList{T}}, x::T)

Enqueue element (last in) the queue. Time complexity O(1).
"""
function enqueue!(q::Queue{DLinkedList{T}}, x::T) where T
    addlast!(q.data, x)
    return nothing
end

"""
    dequeue!(q::Queue)

Remove first element from the queue. Time complexity O(1).
"""
dequeue!(q::Queue) = removefirst!(q.data)

"""
    peek(q::Queue)

Lookup first queue element, don't remove.
"""
peek(q::Queue) = peekfirst(q.data)

"""
    length(q::Queue)

Get size of queue.
"""
length(q::Queue) = length(q.data)

"""
    size(q::Queue)

Get size of queue.
"""
size(q::Queue) = length(q)

"""
    isempty(q::Queue)

Check if queue is empty.
"""
isempty(q::Queue) = length(q) == 0


text:
"""
    BigArray
currently, assume that the array dimension (x,y,z,...) is >= 3
all the manipulation effects in the x,y,z dimension
"""
struct BigArray{D<:AbstractBigArrayBackend, T<:Real} <: AbstractBigArray
    kvStore     :: D
    info        :: Info{T}
    mip         :: Integer
    fillMissing :: Bool
    mode        :: Symbol
end

"""
    BigArray( d::AbstractBigArrayBackend, info::Info;
                  fillMissing::Bool=fillMissing,
                  mode::Symbol=DEFAULT_MODE)

Parameters:
    d: the bigarray storage backend
    info: the info containing the metadata
    mip: mip level. Note that mip 1 is the highest resolution.
            the mip level is like a image pyramid with difference downsampled levels.
    fillMissing: whether fill the missing blocks in the storage backend with zeros or not.
    mode: the io mode with options in {multithreading, sequential, multiprocesses, sharedarray}
"""
function BigArray( d::AbstractBigArrayBackend, info::Info;
                  mip::Integer=1, fillMissing::Bool=DEFAULT_FILL_MISSING,
                  mode::Symbol=DEFAULT_MODE)
    BigArray(d, info, mip, fillMissing, mode)
end

"""
    BigArray(layerPath::AbstractString; mip=1, fillMissing=true, mod::Symbol=DEFAULT_MODE)
"""
function BigArray(layerPath::AbstractString; mip::Integer=1,                            fillMissing::Bool=DEFAULT_FILL_MISSING,
                    mode::Symbol=DEFAULT_MODE)
    if isdir(layerPath) || startswith(layerPath, "file://")
        layerPath = replace(layerPath, "file://"=>"/", count=1)
        d = BinDict(layerPath)
    elseif startswith(layerPath, "gs://")
        d = GSDict(layerPath)
    elseif startswith(layerPath, "s3://")
        d = S3Dict(layerPath)
    else
        @error "only support protocols of {file, gs, s3}, but got: " layerPath
    end
    BigArray(d; mip=mip, fillMissing=fillMissing, mode=mode)
end

@inline function BigArray(d::AbstractBigArrayBackend; mip::Integer=1,
                          fillMissing::Bool=DEFAULT_FILL_MISSING,
                          mode::Symbol=DEFAULT_MODE)
    info = d["info"] |> Info
    return BigArray(d, info, mip, fillMissing, mode)
end

@inline function BigArray( d::AbstractBigArrayBackend, info::Vector{UInt8};
                  mip::Integer=1, fillMissing::Bool=DEFAULT_FILL_MISSING,
                  mode::Symbol=DEFAULT_MODE)
    info = Codings.decode(info, GzipCoding) |> Info
    BigArray(d, string(info), mip, fillMissing, mode)
end

@inline function BigArray( d::AbstractBigArrayBackend, info::AbstractString;
                            mip::Integer=1, fillMissing::Bool=DEFAULT_FILL_MISSING,
                            mode::Symbol=DEFAULT_MODE)
    BigArray(d, Info(info), mip, fillMissing, mode)
end

@inline function BigArray( d::AbstractBigArrayBackend, infoConfig::Dict{Symbol, Any};
                            mip::Integer=1, fillMissing::Bool=DEFAULT_FILL_MISSING,
                            mode::Symbol=DEFAULT_MODE)
    BigArray(d, Info(infoConfig)info, mip, fillMissing, mode)
end

"""
    BigArray(info::Info; mip::Integer=1, fillMissing::Bool=DEFAULT_FILL_MISSING, mode=DEFAULT_MODE)

create a new directory with random name.
this function was designed for test and benchmark.
we need another function the clear the whole array
"""
function BigArray(info::Info; mip::Integer=1,
                  fillMissing::Bool=DEFAULT_FILL_MISSING,
                  mode=DEFAULT_MODE)
    # prepare directory
    layerDir = tempname()
    datasetDir = joinpath(layerDir, string(Infos.get_key(info, 1)))
    mkdir(layerDir)
    mkdir(datasetDir)
    d = BinDict(layerDir)

    # write the info as file
    write(joinpath(layerDir, "info"), JSON.json(Dict(info)))

    return BigArray(d, info, mip, fillMissing, mode)
end

######################### base functions #######################

@inline function Base.ndims(ba::BigArray)
    info = get_info(ba)
    ndims(info)
end

@inline function Base.eltype( ba::BigArray{D,T} ) where {D, T}
    T
end

function Base.size( ba::BigArray )
    # get size according to the size in info file
    return get_volume_size(ba)
end

function Base.size(ba::BigArray, i::Int)
    size(ba)[i]
end

@inline function Base.show(io::IO, ba::BigArray)
    #show(io, ba.kvStore)
    show(io, get_chunk_size(ba))
end

function Base.display(ba::BigArray)
    for field in fieldnames(ba)
        println("$field: $(getfield(ba,field))")
    end
end

"""
    Base.setindex!( ba::BigArray{D,T}, buf::Array{T,N},

setindex with different mode: taskthreads, multithreads, sequential
"""
@inline function Base.setindex!( ba::BigArray{D,T}, buf::Array{T,N},
            idxes::Union{UnitRange, Int, Colon} ... ) where {D,T,N}
    @assert N == ndims(ba)
    if ba.mode == :taskthreads
        setindex_fun! = setindex_taskthreads!
    elseif ba.mode == :multithreads
        setindex_fun! = setindex_multithreads!
    elseif ba.mode == :sequential
        setindex_fun! = setindex_sequential!
    else
        error("only support modes of multithreads, multiprocesses, sharedarray, sequential")
    end
    setindex_fun!(ba, buf, idxes...)
end

@inline function Base.CartesianIndices(ba::BigArray)
    N = ndims(ba)
    offset = get_offset(ba)
    start = offset + one(CartesianIndex{N})
    stop = offset + CartesianIndex(get_volume_size(ba))
    ranges = map((x,y)->x:y, start.I, stop.I)
    return CartesianIndices( ranges )
end

"""
adjust the global and buffer range according to total volume size.
shrink the range stop if the ranges passes the volume boundary.
"""
function adjust_volume_boundary(ba::BigArray, chunkGlobalRange::CartesianIndices,
                                globalRange::CartesianIndices,
                                rangeInChunk::CartesianIndices,
                                rangeInBuffer::CartesianIndices)
    offset = get_offset(ba)
    volumeStop = map(+, offset.I, get_volume_size(ba))
    chunkGlobalRangeStop = [last(chunkGlobalRange).I ...,]
    globalRangeStop = [last(globalRange).I ...,]
    rangeInBufferStop = [last(rangeInBuffer).I ...,]
    rangeInChunkStop = [last(rangeInChunk).I...,]

    for (i,s) in enumerate(volumeStop)
        if chunkGlobalRangeStop[i] > s
            chunkGlobalRangeStop[i] = s
        end
        distanceOverBorder = globalRangeStop[i] - s
        if distanceOverBorder > 0
            globalRangeStop[i] -= distanceOverBorder
            @assert globalRangeStop[i] == s
            @assert globalRangeStop[i] > first(globalRange).I[i]
            rangeInBufferStop[i] -= distanceOverBorder
            rangeInChunkStop[i] -= distanceOverBorder
        end
    end
    start = first(chunkGlobalRange).I
    stop =  (chunkGlobalRangeStop...,)
    chunkGlobalRange = CartesianIndices( map((x,y)->x:y, start, stop) )

    start = first(globalRange).I
    stop = (globalRangeStop...,)
    globalRange = CartesianIndices( map((x,y)->x:y, start, stop) )

    start = first(rangeInBuffer).I
    stop = (rangeInBufferStop...,)
    rangeInBuffer = CartesianIndices( map((x,y)->x:y, start, stop) )

    start = first(rangeInChunk).I
    stop = (rangeInChunkStop...,)
    rangeInChunk = CartesianIndices( map((x,y)->x:y, start, stop) )
    return chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer
end

"""
    Base.getindex( ba::BigArray, idxes::Union{UnitRange, Int}...)

get index with different modes: taskthreads, multithreads, sequential
"""
@inline function Base.getindex( ba::BigArray, idxes::Union{UnitRange, Int}...)
    if ba.mode == :taskthreads
        getindex_fun = getindex_taskthreads
    elseif ba.mode == :multithreads
        getindex_fun = getindex_multithreads
    elseif ba.mode == :sequential
        getindex_fun = getindex_sequential
    else
        error("only support mode of (sharedarray, multi_processes, multithreads, sequential)")
    end
    getindex_fun(ba, idxes...)
end

@inline function get_key_value_store(ba::BigArray) ba.kvStore end
@inline function get_info(ba::BigArray) ba.info end
@inline function get_mip(ba::BigArray) ba.mip end
@inline function get_encoding(ba::BigArray)
    info = get_info(ba)
    mip = get_mip(ba)
    Infos.get_encoding(info, mip)
end

@inline function get_num_channels(ba::BigArray)
    Infos.get_num_channels(get_info(ba))
end

@inline function get_mip_level_name(ba::BigArray)
    Infos.get_key(get_info(ba), get_mip(ba)) |> string
end

@inline function get_chunk_size(ba::BigArray)
    chunkSize = Infos.get_chunk_size(get_info(ba), get_mip(ba))
    if ndims(ba) == 3
        return chunkSize
    else
        return (chunkSize..., get_num_channels(ba))
    end
end

@inline function set_chunk_size(ba::BigArray, chunkSize::Tuple{Int})
    info = get_info(ba)
    # if there are multiple channels, ignore the channel number
    Infos.set_chunk_size(info, chunkSize[1:3])
end

@inline function get_offset(ba::BigArray)
    offset = Infos.get_offset(get_info(ba), get_mip(ba))
    if ndims(ba) == 3
        return offset
    else
        return CartesianIndex(offset.I..., 0)
    end
end

@inline function get_volume_size(ba::BigArray)
    volumeSize = Infos.get_volume_size(get_info(ba), get_mip(ba))
    if ndims(ba) == 3
        return volumeSize
    else
        return (volumeSize..., get_num_channels(ba))
    end
end

@inline function get_mode(ba::BigArray) ba.mode end
@inline function set_mode(ba::BigArray, mode::Symbol)
    ba.mode = mode
end

###################### utils ####################
"""
    get_num_chunks(ba::BigArray, idxes::Union{UnitRange,Int}...)
get number of chunks needed to do cutout from this range
"""
function get_num_chunks(ba::BigArray, idxes::Union{UnitRange, Int}...)
    chunkNum = 0
    baIter = ChunkIterator(idxes, get_chunk_size(ba); offset=get_offset(ba))
	for (blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in baIter
        chunkNum += 1
	end
    chunkNum
end

"""
    list_missing_chunks(ba::BigArray, idxes::Union{UnitRange, Int}...)
list the non-existing keys in the index range
if the returned list is empty, then all the chunks exist in the storage backend.
"""
function list_missing_chunks(ba::BigArray, idxes::Union{UnitRange, Int}...)
    t1 = time()
    sz = map(length, idxes)
    missingChunkList = Vector{CartesianIndices}()
    baIter = ChunkIterator(idxes, get_chunk_size(ba); offset=get_offset(ba))
    @sync begin
        for (blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in baIter
            @async begin
                if !haskey(ba.kvStore, cartesian_range2string(chunkGlobalRange))
                    push!(missingChunkList, chunkGlobalRange)
                end
            end
        end
    end
    missingChunkList
end

function list_missing_chunks(ba::BigArray, keySet::Set{String},
                             idxes::Union{UnitRange, Int}...)
    t1 = time()
    sz = map(length, idxes)
    missingChunkList = Vector{CartesianIndices}()
    baIter = ChunkIterator(idxes, get_chunk_size(ba); offset=get_offset(ba))
    for (blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in baIter
        if !(cartesian_range2string(chunkGlobalRange) in keySet)
            push!(missingChunkList, chunkGlobalRange)
        end
    end
    missingChunkList
end

"""
    commit_info(ba::BigArray)

write info to the storage backend
"""
function commit_info(ba::BigArray)
    info = get_info(ba)
    d = get_key_value_store(ba)
    d["info"] = string(info)
end


text:
#This example has been adopted from https://github.com/IntelLabs/ParallelAccelerator.jl/blob/master/examples/black-scholes/black-scholes.jl

using ArrayFire

function blackscholes_serial(sptprice::AbstractArray{Float32},
                           strike::AbstractArray{Float32},
                           rate::AbstractArray{Float32},
                           volatility::AbstractArray{Float32},
                           time::AbstractArray{Float32})
    logterm = log10( sptprice ./ strike)
    powterm = .5f0 .* volatility .* volatility
    den = volatility .* sqrt(time)
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den
    d2 = d1 .- den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike .* exp(- rate .* time)
    c1 = futureValue .* NofXd2
    call_ = sptprice .* NofXd1 .- c1
    put  = call_ .- futureValue .+ sptprice
end

@inline function cndf2(in::AbstractArray{Float32})
    out = 0.5f0 .+ 0.5f0 .* erf(0.707106781f0 .* in)
    return out
end

function run(iterations)
    sptprice   = Float32[ 42.0 for i = 1:iterations ]
    initStrike = Float32[ 40.0 + (i / iterations) for i = 1:iterations ]
    rate       = Float32[ 0.5 for i = 1:iterations ]
    volatility = Float32[ 0.2 for i = 1:iterations ]
    time       = Float32[ 0.5 for i = 1:iterations ]

    sptprice_gpu = AFArray(sptprice)
    initStrike_gpu = AFArray(initStrike)
    rate_gpu = AFArray(rate)
    volatility_gpu = AFArray(volatility)
    time_gpu = AFArray(time)

    tic()
    put1 = blackscholes_serial(sptprice, initStrike, rate, volatility, time)
    t1 = toq()
    println("Serial checksum: ", sum(put1))
    tic()
    put2 = blackscholes_serial(sptprice_gpu, initStrike_gpu, rate_gpu, volatility_gpu, time_gpu)
    t2 = toq()
    println("Parallel checksum: ", sum(put2))
    return t1, t2
end

function driver()
    srand(0)
    tic()
    iterations = 10^7
    blackscholes_serial(Float32[], Float32[], Float32[], Float32[], Float32[])
    blackscholes_serial(AFArray(Float32[1., 2.]), AFArray(Float32[1., 2.]),
                        AFArray(Float32[1., 2.]), AFArray(Float32[1., 2.]), AFArray(Float32[1., 2.]))
    println("SELFPRIMED ", toq())
    tserial, tparallel = run(iterations)
    println("Time taken for CPU = $tserial")
    println("Time taken for GPU = $tparallel")
    println("Speedup = $(tserial / tparallel)")
    println("CPU rate = ", iterations / tserial, " opts/sec")
    println("GPU rate = ", iterations / tparallel, " opts/sec")
end
driver()


text:
using SparseArrays
println("testing registration")

src_nodes = [1.0 1.0;
			2.0 3.0;
			3.0 3.0]
dst_nodes = [5.0 1.0;
			6.0 3.0;
			7.0 3.0]
# incidence = [1 -1 0;
# 			1 0 -1;
# 			0 1 -1]
edges = spzeros(Int64, 3, 3)
edges[1:2, 1] = [1, -1]
edges[1, 2] = 1
edges[3, 2] = -1
edges[2:3, 3] = [1, -1]

mesh = Mesh(src_nodes, dst_nodes, edges)
mesh_indices = [1, 2]
matches = Matches(mesh_indices, src_nodes[1:2,:], dst_nodes[1:2,:])
new_mesh = matches_to_mesh(matches, mesh)
true_mesh = Mesh(src_nodes[1:2,:], dst_nodes[1:2,:], edges[1:2, 1])
@test new_mesh.src_nodes == true_mesh.src_nodes
@test new_mesh.dst_nodes == true_mesh.dst_nodes
@test new_mesh.edges == true_mesh.edges


text:
using BAT
using Test

using Distributions
using ValueShapes
using IntervalSets
using LinearAlgebra: Diagonal, ones

@testset "bridge_sampling_integration" begin
    function test_integration(algorithm::IntegrationAlgorithm, title::String,
                              dist::Distribution; val_expected::Real=1.0,
                              val_rtol::Real=3.5, err_max::Real=0.2)
        @testset "$title" begin
            samplingalg = MCMCSampling(
                mcalg = MetropolisHastings(),
                trafo = NoDensityTransform(),
                nsteps = 2*10^5,
                burnin = MCMCMultiCycleBurnin(nsteps_per_cycle = 10^5, max_ncycles = 60)
            )
            sample = bat_sample(dist, samplingalg).result
            sd = SampledDensity(dist, sample)
            sample_integral = bat_integrate(sd, algorithm).result

            @test isapprox(sample_integral.val, val_expected, atol=val_rtol*sample_integral.err)
            @test sample_integral.err < err_max
        end
    end
    test_integration(BridgeSampling(trafo=NoDensityTransform()), "funnel distribution", BAT.FunnelDistribution(), val_rtol = 15)
    test_integration(BridgeSampling(trafo=NoDensityTransform()), "multimodal student-t distribution", BAT.MultimodalStudentT(), val_rtol = 50)
    test_integration(BridgeSampling(trafo=NoDensityTransform()), "Gaussian shell", BAT.GaussianShell(), val_rtol = 15)
    test_integration(BridgeSampling(trafo=NoDensityTransform()), "MvNormal", MvNormal(Diagonal(ones(5))), val_rtol = 15)

end


text:
include("algebra.jl")
include("CachedOperator.jl")
include("FiniteOperator.jl")
include("TransposeOperator.jl")
include("InterlaceOperator.jl")
include("OperatorFunction.jl")


text:
function test_convert(mdp::PedCarMDP)
    @showprogress for i=1:n_states(mdp)
        s = ind2state(mdp, i)
        svec = convert_s(Vector{Float64}, s, mdp)
        s_reconstructed = convert_s(PedCarMDPState, svec, mdp)
        if !isapprox(s_reconstructed, s)
            @debug("index $i, state $s")
            @debug("reconstructed state: $s_reconstructed")
            return false
        end
    end
    return true
end

function Base.isapprox(s1::PedCarMDPState, s2::PedCarMDPState, kwargs...)
    return isapprox(s1.crash, s2.crash, kwargs...) &&
           isapprox(s1.ego, s2.ego, kwargs...) &&
           isapprox(s1.car, s2.car, kwargs...) &&
           isapprox(s1.ped, s2.ped, kwargs...) &&
           s1.route == s2.route
end

@testset begin "convert_s"
    params = UrbanParams(nlanes_main=1,
                        crosswalk_pos =[VecSE2(6, 0., pi/2), VecSE2(-6, 0., pi/2), VecSE2(0., -5., 0.)],
                        crosswalk_length =  [14.0, 14., 14.0],
                        crosswalk_width = [4.0, 4.0, 3.1],
                        stop_line = 22.0)
    env = UrbanEnv(params=params)
    mdp = PedCarMDP(env=env, pos_res=4.0, vel_res=3., ped_birth=0.3, car_birth=0.3)

    @test test_convert(mdp)
end


text:
SX(x::T) where {T <: Number} = casadi.SX(x)
SX(x::T) where {T <: AbstractVecOrMat{SX}} = convert(SX, x)
SX(x::AbstractVecOrMat{T}) where {T <: AbstractFloat} = casadi.SX(x)
SX(x::AbstractVecOrMat{T}) where {T <: Integer} = casadi.SX(x)
SX(x::AbstractString) = casadi.SX.sym(x)
SX(x::AbstractString, i1::Integer) = casadi.SX.sym(x, i1)
SX(x::AbstractString, i1::Integer, i2::Integer) = casadi.SX.sym(x, i1, i2)
SX(i1::Integer, i2::Integer) = casadi.SX(i1, i2)

MX(x::T) where {T <: Number} = casadi.MX(x)
MX(x::T) where {T <: AbstractVecOrMat{MX}} = convert(MX, x)
MX(x::AbstractVecOrMat{T}) where {T <: AbstractFloat} = casadi.MX(x)
MX(x::AbstractVecOrMat{T}) where {T <: Integer} = casadi.MX(x)
MX(x::AbstractString) = casadi.MX.sym(x)
MX(x::AbstractString, i1::Integer) = casadi.MX.sym(x, i1)
MX(x::AbstractString, i1::Integer, i2::Integer) = casadi.MX.sym(x, i1, i2)
MX(i1::Integer, i2::Integer) = casadi.MX(i1, i2)

convert(::Type{C}, s::AbstractString) where C <: CasadiSymbolicObject = C(s)


text:
# Project Files for MAC wave project

#
# include project function files
include("macWave.jl")
include("grids.jl")
include("diffusion.jl")
include("stratification.jl")
include("dbdt.jl")
include("magnetic_field.jl")
include("restoring.jl")
include("velocityphi.jl");
include("flow.jl");

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
# Test basic variable definition
@randomvariable(m2, x, mean = 0)
@recoursevariable(m2, y)
@variable(m2, z)
true


text:
export MaxTimeoutEnv

mutable struct MaxTimeoutEnv{E<:AbstractEnv} <: AbstractEnvWrapper
    env::E
    max_t::Int
    current_t::Int
end

"""
    MaxTimeoutEnv(env::E, max_t::Int; current_t::Int = 1)

Force `is_terminated(env)` return `true` after `max_t` interactions.
"""
MaxTimeoutEnv(env::E, max_t::Int; current_t::Int = 1) where {E<:AbstractEnv} =
    MaxTimeoutEnv(env, max_t, current_t)

function (env::MaxTimeoutEnv)(args...; kwargs...)
    env.env(args...; kwargs...)
    env.current_t = env.current_t + 1
end

RLBase.is_terminated(env::MaxTimeoutEnv) =
    (env.current_t > env.max_t) || is_terminated(env.env)

function RLBase.reset!(env::MaxTimeoutEnv)
    env.current_t = 1
    RLBase.reset!(env.env)
end


text:
module Jagot

include("misc.jl")
include("plotting.jl")
include("spectra.jl")
include("utils.jl")

export plotting, physics, spectra

end # module


text:
@i function (:-=)(gaussian_log)(y!::GVar{T}, x::GVar{T}) where T
	y!.x -= gaussian_log(x.x)
	@routine @invcheckoff begin
		exp_x ← zero(x)
		jac ← zero(x)
		exp_x += exp(-x)
		exp_x += 1
		jac += 1/exp_x
	end
	x.g += y!.g * jac
	~@routine
end

@i function (:-=)(gaussian_nlog)(y!::GVar{T}, x::GVar{T}) where T
	y!.x -= gaussian_nlog(x.x)
	@routine @invcheckoff begin
		exp_x ← zero(x)
		jac ← zero(x)
		exp_x += exp(-x)
		exp_x -= 1
		jac -= 1/exp_x
	end
	x.g += y!.g * jac
	~@routine
end

@i function :(-=)(convert)(out!::GVar{Tx, Tg}, y::ULogarithmic) where {Tx, Tg}
	out! -= exp(y.log)
end


text:
function sledge_down(trees, dx, dy)
    is_tree(i, j) = trees[i, j] == '#'

    x, y = 1, 1
    cnt = 0
    nr_rows, nr_cols = size(trees)
    while x <= nr_rows
        if is_tree(x, y)
            cnt += 1
        end
        x += dx
        y = ((y + dy - 1) % nr_cols) + 1
    end
    cnt
end

function tree_product(trees)
    slopes = [(1, 1), (1, 3), (1, 5), (1, 7), (2, 1)]
    prod(map(((dx, dy),) -> sledge_down(trees, dx, dy), slopes))
end

function main()
    lines = readlines("../resources/input_03.txt")
    trees = permutedims(hcat(map(collect, lines)...))
    result = tree_product(trees)
    println(result)
end

main()


text:
"""
    trapezint(x::AbstractVector, y::AbstractVector)
    trapezint(step::Real, y::AbstractVector)
Calculate the integral of `y` versus `x` using the trapezoidal rule.
The second method assumes equally-spaced `y`.
"""
function trapezint{X<:Real, Y<:Real}(x::AbstractVector{X}, y::AbstractVector{Y})
  integral::Float64 = 0.0
  @inbounds for i in 1:(length(x)-1)
    integral += (x[i+1] - x[i])*(y[i+1] + y[i])
  end
  integral *= 0.5
end
function trapezint{Y<:Real}(step::Real, y::AbstractVector{Y})
  integral::Float64 = 0.0
  @inbounds for i in 1:(length(y)-1)
    integral += step*(y[i+1] + y[i])
  end
  integral *= 0.5
end

"""
    fdderiv(x::AbstractVector, y::AbstractVector)
    fdderiv(step::Real, y::AbstractVector)
Calculate the derivative of `y` versus `x` using the forward difference approximation.
The second method assumes equally-spaced `y`. Notice that the last element of the
derivative cannot be computed and therefore it is equalled to the second-to-last.
"""
function fdderiv(x::AbstractVector, y::AbstractVector)
  der = zeros(y)
  @inbounds @simd for i in 1:length(y)-1
    der[i] = (y[i+1]-y[i])/(x[i+1] - x[i])
  end
  der[end] = der[end-1]
  return der
end
function fdderiv(x::Real, y::AbstractVector)
  der = zeros(y)
  @inbounds @simd for i in 1:length(y)-1
    der[i] = (y[i+1]-y[i])/(x)
  end
  der[end] = der[end-1]
  return der
end

"""
    fdderiv!(der::AbstractVector, x::AbstractVector, y::AbstractVector)
    fdderiv!(der::AbstractVector, step::Real, y::AbstractVector)
Calculate the derivative of `y` versus `x` using the forward difference approximation
and write it in-place for `der`.
The second method assumes equally-spaced `y`. Notice that the last element of
the derivative
cannot be computed and therefore it is equalled to the second-to-last.
"""
function fdderiv!(der::AbstractVector, x::AbstractVector, y::AbstractVector)
  @inbounds @simd for i in 1:length(y)-1
    der[i] = (y[i+1]-y[i])/(x[i+1] - x[i])
  end
  der[end] = der[end-1]
  return der
end
function fdderiv!(der::AbstractVector, x::Real, y::AbstractVector)
  @inbounds @simd for i in 1:length(y)-1
    der[i] = (y[i+1]-y[i])/(x)
  end
  der[end] = der[end-1]
  return der
end

"""Return the order of magnitude of a number"""
mag(x::Real) = ceil(log10(x))


text:
mutable struct Vertex
    vector
    index
    body
    internal
    contact
end

Vertex(vector, index, body) = Vertex(vector, index, body, false, nothing)

mutable struct Contact
    vertex
    normal
    tangent
end

Contact(vertex) = Contact(vertex, 0., 0.)

function vertices(vectors, body)
    vertices = Array{Vertex,1}()
    for (index, vector) in enumerate(vectors)
        vertex = Vertex(vector, index, body)
        contact = Contact(vertex)
        vertex.contact = contact
        push!(vertices, vertex)
    end
    return vertices
end

"""
    center(vertices::Array{Vertex,1})

Compute the first moment ("center of mass") of the vertices.
"""
function center(vertices::Array{Vertex,1})
    A = area(vertices)
    M = Vector(0., 0.)
    for index in eachindex(vertices)
        vector = vertices[index].vector
        next_vector = vertices[index % length(vertices) + 1].vector
        cross_product = det(vector, next_vector)
        M += (vector + next_vector) * cross_product
    end
    return M / 6. / A
end

"""
    mean(vertices::Array{Vertex,1})

Compute the average position of the vertices.
"""
function mean(vertices::Array{Vertex,1})
    average = Vector(0., 0.)
    for vertex in vertices
        average += vertex.vector
    end
    return average / length(vertices)
end

"""
    area(vertices::Array{Vertex,1})

Compute the area bound by the convex hull of the vertices.
"""
function area(vertices::Array{Vertex,1})
    A = 0.
    prev_vector = vertices[end].vector
    for index in eachindex(vertices)
        vector = vertices[index].vector
        A += det(prev_vector, vector)
        prev_vector = vector
    end
    return A / 2.
end

"""
    inertia(vertices::Array{Vertex,1})

Compute the second moment ("moment of inertia") of the vertices.
"""
function inertia(vertices::Array{Vertex,1}, mass)
    numerator = 0.
    denominator = 0.
    for index in eachindex(vertices)
        vector = vertices[index].vector
        next_vector = vertices[index % length(vertices) + 1].vector
        cross_product = det(next_vector, vector)
        numerator += abs(cross_product) * (dot(next_vector, next_vector) + dot(next_vector, vector) + dot(vector, vector))
        denominator += cross_product
    end
    return mass / 6 * numerator / denominator
end

"""
    translate!(vertices::Array{Vertex,1}, translation::Vector)

Move vertices by translation.
"""
function translate!(vertices::Array{Vertex,1}, translation::Vector)
    for index in eachindex(vertices)
        vertices[index].vector += translation
    end
end

"""
    rotate!(vertices::Array{Vertex,1}, angle::Vector, point::Vector)

Rotate vertices by angle about a point.
"""
function rotate!(vertices::Array{Vertex,1}, angle, point::Vector)
    angle == 0. && return
    for index in eachindex(vertices)
        rotate!(vertices[index].vector, angle, point)
    end
end

"""
    contains(vertices::Array{Vertex,1}, point::Vector)

Check if a point is inside `vertices`.
"""
function contains(vertices::Array{Vertex,1}, point::Vector)
    for index in eachindex(vertices)
        vector = vertices[index].vector
        next_vector = vertices[index % length(vertices) + 1].vector
        if det(point - vector, next_vector - vector) > 0
            return false
        end
    end
    return true
end


text:
module GaussianSlab
function get_fn(argd::Dict)
    center = [float(x) for x in argd["center"]]
    sigma  = float(argd["sigma"])
    k = float(argd["k"])
    upper = float(argd["upper"])
    lower = float(argd["lower"])
    raw_direction = [float(x) for x in argd["direction"]]
    direction = raw_direction/norm(raw_direction)
    function fn_compute(coordinates)
        d = dot(coordinates - center, direction)
        if lower <= d <= upper
            return k * exp(-(d)^2 /(2 *sigma))
        else
            return 0.0
        end
    end

    return fn_compute
end

export get_fn
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape