text:
include("LatticeQCD.jl")
using .LatticeQCD
run_wizard()


text:
mutable struct Font{T} <: AbstractResource
    ptr::Ptr{T}
    render_ptr::Ptr{SDL.Renderer}
    filename::String
    fontsize::Int
    cache::Dict{UInt64,Texture}

    function Font{T}(ptr, render_ptr, filename, fontsize) where {T}
        self = new{T}(ptr, render_ptr, filename, fontsize, Dict{UInt64,Texture}())
        finalizer(destroy!, self)
        return self
    end
end
Font(ptr::Ptr{T}, render_ptr, filename, fontsize) where {T} = Font{T}(ptr, render_ptr, filename, fontsize)

function Font(resources::Resources, filename::String; fontsize=12)
    if filename in keys(resources)
        @debug("resource already loaded: '$filename'")
        return resources[filename]::Font
    end

    ptr = SDL.TTF_OpenFont(filename, fontsize)
    self = Font(ptr, resources.render_ptr, filename, fontsize)
    finalizer(destroy!, self)
    resources[filename] = self
    return self
end

function destroy!(font::Font{SDL.TTF_Font})
    SDL.TTF_CloseFont(font.ptr)
    font.ptr = C_NULL
    empty!(font.cache)
    return nothing
end

function text(font::Font, txt::String; color=colorant"#303030", halign=0.0, valign=0.0)
    _hash = hash((txt, color, halign, valign))
    if _hash in keys(font.cache)
        return font.cache[_hash]::Texture
    end

    texture = Texture(font, txt, color=color, halign=halign, valign=valign)
    font.cache[_hash] = texture
    return texture
end

function Texture(font::Font, text::String; color=colorant"#303030", halign=0.0, valign=0.0)
    sdl_surface = SDL.TTF_RenderText_Blended(font.ptr, text, SDL.Color(rgba_fromcolor(color)...))
    return Texture(font.render_ptr, sdl_surface, halign=halign, valign=valign)
end

ext: jl


text:
module INP

export InpStiffness

using ...TopOptProblems: Metadata, StiffnessTopOptProblem, QuadraticHexahedron
using Ferrite
using ....TopOpt.Utilities: find_black_and_white, find_varind
import ...TopOptProblems:
    nnodespercell,
    getE,
    getν,
    getgeomorder,
    getdensity,
    getpressuredict,
    getcloaddict,
    getfacesets

include(joinpath("Parser", "Parser.jl"))
using .Parser

include("inpstiffness.jl")

end


text:
using KrylovKit

struct EquippedOperator{O,W}
    operator::O
    working_memory::W
end
Base.eltype(::EquippedOperator{O}) where {O} = eltype(O)
Base.eltype(::Type{<:EquippedOperator{O}}) where {O} = eltype(O)

function equip(operator, vector; style=IsDeterministic())
    if eltype(operator) === valtype(vector)
        wm = WorkingMemory(vector; style)
    else
        wm = WorkingMemory(similar(vector, eltype(operator)); style)
    end
    return EquippedOperator(operator, wm)
end

function (eo::EquippedOperator)(src::TVec)
    op = eo.operator
    wm = eo.working_memory
    style = wm.style
    dst = similar(src)
    mul!(dst, op, src, wm, style)
end

function KrylovKit.eigsolve(
    ham::AbstractHamiltonian, dv::TVec, howmany::Int, which::Symbol=:SR;
    issymmetric=eltype(ham) <: Real && LOStructure(ham) === IsHermitian(),
    ishermitian=LOStructure(ham) === IsHermitian(),
    verbosity=0,
    style=IsDeterministic(),
    kwargs...
)
    eo = equip(ham, dv; style)
    return eigsolve(eo, dv, howmany, which; issymmetric, verbosity, kwargs...)
end


text:
function count_nucleotides(strand::AbstractString)

    a = 0
    c = 0
    g = 0
    t = 0

    for i in strand
        if i == 'A'
            a += 1
        elseif i == 'C'
            c += 1
        elseif i == 'G'
            g += 1
        elseif i == 'T'
            t += 1
        else
            return throw(DomainError())
        end
    end
    return Dict( 'A' => a, 'C' => c, 'G' => g, 'T' => t)
end

count_nucleotides("GAGAGACGTT")

max_stars_repo_name: Devlious/LearningJulia

max_issues_repo_name: Devlious/LearningJulia

max_forks_repo_name: Devlious/LearningJulia


text:
using Polylogarithms
using Documenter

DocMeta.setdocmeta!(Polylogarithms, :DocTestSetup, :(using Polylogarithms); recursive=true)

makedocs(;
    modules=[Polylogarithms],
    authors="Matthew Roughan <matthew.roughan@adelaide.edu.au> and contributors",
    repo="https://github.com/mroughan/Polylogarithms.jl/blob/{commit}{path}#L{line}",
    sitename="Polylogarithms.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://matthew.roughan@adelaide.edu.au.github.io/Polylogarithms.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mroughan/Polylogarithms.jl",
)


text:
#///////////////////////////////////////
#// File Name: bic_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/02/25
#// Description: Test code for src/bic.jl
#///////////////////////////////////////

using Distributions
using GeometricalPredicates
using LinearAlgebra
using Random
using RobotOS
using VoronoiCells

@testset "Buffered Input Cell Test" begin
# # vertex_direction test
generator = IndexablePoint2D(0.0, 0.0, 1);
vertex = Point2D(1.0, 1.0);
@test vertex_direction(generator, vertex) ≈ pi/4;
# # sort_vertices test
vertices = [Point2D(1.0, 1.0),
            Point2D(1.0, -1.0),
            Point2D(-1.0, -1.0),
            Point2D(-1.0, 1.0)];
vertices_sorted = sort_vertices(generator, vertices);
@test vertices_sorted == [Point2D(-1.0, -1.0),
                          Point2D(1.0, -1.0),
                          Point2D(1.0, 1.0),
                          Point2D(-1.0, 1.0)];
# # scale_and_shift test
x_test = 10.0;
@test scale_and_shift(x_test, scale_factor=0.1;
                      subtract=2.0, offset=1.0) ≈ 1.8;
# # find_extrema test
input_positions_test = [[-1.0, -1.0],
                        [1.0,  -1.0],
                        [1.0,   1.0],
                        [-1.0,  1.0]];
xmin, xmax, ymin, ymax = find_extrema(input_positions_test);
@test (xmin, xmax) == (-8.0, 8.0);
@test (ymin, ymax) == (-8.0, 8.0);
# # retract_cell test
vertices = input_positions_test
cell = retract_cell(generator, vertices, 0.2)
points = getpoints(cell);
@test points == [Point2D(-0.8,-0.8),
                 Point2D(0.8, -0.8),
                 Point2D(0.8,  0.8),
                 Point2D(-0.8, 0.8)];
# # get_input_positions test
e_init = RobotState([0.0, 0.0, 0.0, 0.0]);
ap_dict = Dict{String, Vector{Float64}}()
ap_dict["Pedestrian/1"] = [1.0,   1.0];
ap_dict["Pedestrian/2"] = [1.0,  -1.0];
ap_dict["Pedestrian/3"] = [-1.0,  1.0];
ap_dict["Pedestrian/4"] = [-1.0, -1.0];
w_init = WorldState(e_init, ap_dict);
input_positions = get_input_positions(w_init);
@test input_positions == [[0.0,   0.0],
                          [1.0,   1.0],
                          [1.0,  -1.0],
                          [-1.0,  1.0],
                          [-1.0, -1.0]];
# # compute_buffered_voronoi_cells test
bvc = compute_buffered_voronoi_cells(w_init, 0.2)
points = getpoints(bvc[1]);
@test all(isapprox.((getx(points[1]), gety(points[1])), (0.0, -(1 - 0.2*sqrt(2))), atol=1e-8));
@test all(isapprox.((getx(points[2]), gety(points[2])), (1 - 0.2*sqrt(2), 0.0), atol=1e-8));
@test all(isapprox.((getx(points[3]), gety(points[3])), (0.0, 1 - 0.2*sqrt(2)), atol=1e-8));
@test all(isapprox.((getx(points[4]), gety(points[4])), (-(1 - 0.2*sqrt(2)), 0.0), atol=1e-8));
# # compute_buffered_voronoi_cells test (empty cell case)
bvc2 = compute_buffered_voronoi_cells(w_init, 1.0)
points2 = getpoints(bvc2[1]);
@test length(points2) == 1;
@test all((getx(points2[1]), gety(points2[1])) .≈ (0.0, 0.0))
# # compute_buffered_input_cell test
bic = compute_buffered_input_cell(bvc[1], w_init.e_state, 0.4);
points = getpoints(bic)
@test all(isapprox.((getx(points[1]), gety(points[1])), (0.0, -2/(0.4^2)*(1 - 0.2*sqrt(2))), atol=1e-8));
@test all(isapprox.((getx(points[2]), gety(points[2])), (2/(0.4^2)*(1 - 0.2*sqrt(2)), 0.0), atol=1e-8));
@test all(isapprox.((getx(points[3]), gety(points[3])), (0.0, 2/(0.4^2)*(1 - 0.2*sqrt(2))), atol=1e-8));
@test all(isapprox.((getx(points[4]), gety(points[4])), (-2/(0.4^2)*(1 - 0.2*sqrt(2)), 0.0), atol=1e-8));

projected = project([100., 0.], bic);
@test all(isapprox.(projected, [2/(0.4^2)*(1 - 0.2*sqrt(2)), 0.], atol=1e-8));
end

@testset "LQ Tracking Test" begin
# # lq_tracking test
target_trajectory = Trajectory2D()
target_trajectory[Time(0.0)] = [0.0, 0.0];
target_trajectory[Time(10.0)] = [10.0, 0.0];

e_init = RobotState([0.0, 0.0, 1.0, 0.0]);
ap_dict = Dict{String, Vector{Float64}}()
ap_dict["Pedestrian/1"] = [1.0,   1.0];
w_init = WorldState(e_init, ap_dict);

Cep = Matrix(0.5I, 2, 2);
Cu = Matrix(0.2I, 2, 2);
β_pos = 0.6;
β_col = 0.4;
α_col = 100.0;
λ_col = 0.2;
σ_risk = 0.0;
cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);

dtc = 0.02;
dto = 0.4;
prediction_steps = 12;
num_samples = 30;
sim_param = SimulationParameter(dtc, dto, prediction_steps, num_samples, cost_param)

u_array, x_array = lq_tracking(w_init, target_trajectory, sim_param);
for ii = 1:length(u_array)
    @test all(isapprox.(u_array[ii], [0.0, 0.0], atol=1e-8))
    @test all(isapprox.(x_array[ii], [0.4*(ii - 1),  0.0, 1.0, 0.0], atol=1e-8))
end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

julia_version = v"1.5.3"

name = "jlqml"
version = v"0.1.5"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/barche/jlqml.git", "a36226c77ad712611ae142b9a5855e33eaddbb53"),
]

# Bash recipe for building across all platforms
script = raw"""
if test -f "$prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake"; then
    sed -i 's/_qt5gui_find_extra_libs.*AGL.framework.*//' $prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake
fi

# Override compiler ID to silence the horrible "No features found" cmake error
if [[ $target == *"apple-darwin"* ]]; then
  macos_extra_flags="-DCMAKE_CXX_COMPILER_ID=AppleClang -DCMAKE_CXX_COMPILER_VERSION=10.0.0 -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
fi

mkdir build
cd build
cmake \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DJulia_PREFIX=${prefix} \
    -DQt5Core_DIR=$prefix/lib/cmake/Qt5Core \
    -DQt5Quick_DIR=$prefix/lib/cmake/Qt5Quick \
    -DQt5Svg_DIR=$prefix/lib/cmake/Qt5Svg \
    -DQt5Widgets_DIR=$prefix/lib/cmake/Qt5Widgets \
    $macos_extra_flags \
    ../jlqml/
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
install_license $WORKSPACE/srcdir/jlqml*/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)

filter!(p -> libc(p) != "musl", platforms) # Qt_jll is currently not available for muslc

# The products that we will ensure are always built
products = [
    LibraryProduct("libjlqml", :libjlqml),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"),
    Dependency("Qt_jll"),
    BuildDependency("Libglvnd_jll"),
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"8",
    julia_compat = "$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
__precompile__()

module LibPNG

using ColorTypes,
      ImageCore,
      FixedPointNumbers,
      MicroLogging

export readimage, writeimage

const depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("LibPNG not properly installed. Please run Pkg.build(\"LibPNG\") then restart Julia.")
end

include("functions.jl")
include("io.jl")

end # module


text:
module ClimaCommsMPI

using ClimaComms

using KernelAbstractions
using MPI

struct MPICommsContext <: ClimaComms.AbstractCommsContext
    mpicomm::MPI.Comm
end
MPICommsContext() = MPICommsContext(MPI.COMM_WORLD)

function ClimaComms.init(ctx::MPICommsContext)
    if !MPI.Initialized()
        MPI.Init()
    end

    return ClimaComms.mypid(ctx), ClimaComms.nprocs(ctx)
end

ClimaComms.mypid(ctx::MPICommsContext) = MPI.Comm_rank(ctx.mpicomm) + 1
ClimaComms.iamroot(ctx::MPICommsContext) = ClimaComms.mypid(ctx) == 1
ClimaComms.nprocs(ctx::MPICommsContext) = MPI.Comm_size(ctx.mpicomm)

ClimaComms.barrier(ctx::MPICommsContext) = MPI.Barrier(ctx.mpicomm)

ClimaComms.reduce(ctx::MPICommsContext, val, op) =
    MPI.Reduce(val, op, 0, ctx.mpicomm)

function ClimaComms.gather(ctx::MPICommsContext, array)
    dims = size(array)
    lengths = MPI.Gather(dims[end], 0, ctx.mpicomm)
    if ClimaComms.iamroot(ctx)
        dimsout = (dims[1:(end - 1)]..., sum(lengths))
        arrayout = similar(array, dimsout)
        recvbuf = MPI.VBuffer(arrayout, lengths .* prod(dims[1:(end - 1)]))
    else
        recvbuf = nothing
    end
    MPI.Gatherv!(array, recvbuf, 0, ctx.mpicomm)
end

ClimaComms.abort(ctx::MPICommsContext, status::Int) =
    MPI.Abort(ctx.mpicomm, status)

# We could probably do something fancier here?
# Would need to be careful as there is no guarantee that all ranks will call
# finalizers at the same time.
const TAG = Ref(Cint(0))
function newtag(ctx::MPICommsContext)
    TAG[] = tag = mod(TAG[] + 1, 32767) # TODO: this should query MPI_TAG_UB attribute (https://github.com/JuliaParallel/MPI.jl/pull/551)
    if tag == 0
        @warn("MPICommsMPI: tag overflow")
    end
    return tag
end

"""
    MPISendRecvGraphContext

A simple ghost buffer implementation using MPI `Isend`/`Irecv` operations.
"""
mutable struct MPISendRecvGraphContext <: ClimaComms.AbstractGraphContext
    ctx::MPICommsContext
    tag::Cint
    send_bufs::Vector{MPI.Buffer}
    send_ranks::Vector{Cint}
    send_reqs::Vector{MPI.Request}
    recv_bufs::Vector{MPI.Buffer}
    recv_ranks::Vector{Cint}
    recv_reqs::Vector{MPI.Request}
end

function ClimaComms.graph_context(
    ctx::MPICommsContext,
    send_array,
    send_lengths,
    send_pids,
    recv_array,
    recv_lengths,
    recv_pids,
)
    @assert length(send_pids) == length(send_lengths)
    @assert length(recv_pids) == length(recv_lengths)

    tag = newtag(ctx)

    send_bufs = MPI.Buffer[]
    total_len = 0
    for len in send_lengths
        buf = MPI.Buffer(view(send_array, (total_len + 1):(total_len + len)))
        push!(send_bufs, buf)
        total_len += len
    end
    send_ranks = Cint[pid - 1 for pid in send_pids]
    send_reqs = MPI.Request[MPI.REQUEST_NULL for _ in send_ranks]

    recv_bufs = MPI.Buffer[]
    total_len = 0
    for len in recv_lengths
        buf = MPI.Buffer(view(recv_array, (total_len + 1):(total_len + len)))
        push!(recv_bufs, buf)
        total_len += len
    end
    recv_ranks = Cint[pid - 1 for pid in recv_pids]
    recv_reqs = MPI.Request[MPI.REQUEST_NULL for _ in recv_ranks]

    MPISendRecvGraphContext(
        ctx,
        tag,
        send_bufs,
        send_ranks,
        send_reqs,
        recv_bufs,
        recv_ranks,
        recv_reqs,
    )
end

function ClimaComms.start(
    ghost::MPISendRecvGraphContext;
    dependencies = nothing,
)
    if !all(r -> r == MPI.REQUEST_NULL, ghost.recv_reqs)
        error("Must finish() before next start()")
    end
    # post receives
    for n in 1:length(ghost.recv_bufs)
        ghost.recv_reqs[n] = MPI.Irecv!(
            ghost.recv_bufs[n],
            ghost.recv_ranks[n],
            ghost.tag,
            ghost.ctx.mpicomm,
        )
    end
    # post sends
    for n in 1:length(ghost.send_bufs)
        ghost.send_reqs[n] = MPI.Isend(
            ghost.send_bufs[n],
            ghost.send_ranks[n],
            ghost.tag,
            ghost.ctx.mpicomm,
        )
    end
end

function ClimaComms.progress(ghost::MPISendRecvGraphContext)
    MPI.Iprobe(MPI.MPI_ANY_SOURCE, ghost.tag, ghost.ctx.mpicomm)
end

function ClimaComms.finish(
    ghost::MPISendRecvGraphContext;
    dependencies = nothing,
)
    # wait on previous receives
    MPI.Waitall!(ghost.recv_reqs)
    # ensure that sends have completed
    # TODO: these could be moved to start()? but we would need to add a finalizer to make sure they complete.
    MPI.Waitall!(ghost.send_reqs)
end

end # module


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute()::BigInteger = 28_433 * powermod(BigInteger(2), 7_830_457, 10_000_000_000) % 10_000_000_000 + 1

compute()

@benchmark compute()

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
import PTChain
using PlotChain

function load_sample()
    filename = "si1032.tsp"
    f = open(filename)
    s = readstring(f)
    numbers = map(x -> parse(Int, x), split(s))
    a = 1032
    b = 1
    out = zeros(a, a)
    for i = 1:a
        out[1:i - 1] = zeros(i - 1)
        for j = i:a
            out[i, j] = numbers[b]
            b = b + 1
        end
    end
    close(f)
    return out + out'
end

function score_tour(graph, tour)
    val = 0
    assert(size(graph, 1) == length(tour))
    for i = 1:(length(tour) - 1)
        val = graph[tour[i], tour[i + 1]] + val
    end
    return val + graph[tour[end], tour[1]]
end

function swap_two(tour, mat)
    new_tour = copy(tour)
    randone = rand(1:length(tour))
    randtwo = rand(1:length(tour))
    new_tour[randone] = tour[randtwo]
    new_tour[randtwo] = tour[randone]
    len = length(tour)
    randone_n = randone == len ? 1 : randone + 1
    randone_p = randone == 1 ? len : randone - 1
    randtwo_n = randtwo == len ? 1 : randtwo + 1
    randtwo_p = randtwo == 1 ? len : randtwo - 1

    diff = 0
    diff += mat[new_tour[randone_p], new_tour[randone]] - mat[tour[randone_p], tour[randone]]
    diff += mat[new_tour[randone], new_tour[randone_n]] - mat[tour[randone], tour[randone_n]]
    diff += mat[new_tour[randtwo_p], new_tour[randtwo]] - mat[tour[randtwo_p], tour[randtwo]]
    diff += mat[new_tour[randtwo], new_tour[randtwo_n]] - mat[tour[randtwo], tour[randtwo_n]]

    new_tour, diff
end

function reverse_section(tour, mat)
    new_tour = copy(tour)
    randone, randtwo = minmax(rand(1:length(tour)), rand(1:length(tour)))
    new_tour[randone:randtwo] = tour[randtwo:-1:randone]
    len = length(tour)
    randone_n = randone == len ? 1 : randone + 1
    randone_p = randone == 1 ? len : randone - 1
    randtwo_n = randtwo == len ? 1 : randtwo + 1
    randtwo_p = randtwo == 1 ? len : randtwo - 1
    diff = 0
    diff += mat[new_tour[randone_p], new_tour[randone]] - mat[tour[randone_p], tour[randone]]
    diff += mat[new_tour[randone], new_tour[randone_n]] - mat[tour[randone], tour[randone_n]]
    diff += mat[new_tour[randtwo_p], new_tour[randtwo]] - mat[tour[randtwo_p], tour[randtwo]]
    diff += mat[new_tour[randtwo], new_tour[randtwo_n]] - mat[tour[randtwo], tour[randtwo_n]]
    new_tour, diff
end

function proposal_function(tour, mat)
    if rand() < 0.5
        return swap_two(tour, mat)
    else
       return reverse_section(tour, mat)
    end
end

function sim_tour()
    n = 1032
    start_x = randperm(n)
    mat = load_sample()
    reps = 1000
    temps = 100 * (0.95).^(0:0.5:200)
    f(x) = score_tour(mat, x)
    update(x) = proposal_function(x, mat)
    min_val, x, stats, results = PTChain.simulatedAnnealing(start_x, f, temps, update, fill(reps, length(temps)))
    return stats
end

function st_tour(micro_steps, macro_steps)
    n = 1032
    start_x = randperm(n)
    mat = load_sample()
    reps = 10
    temps = collect(0:2:100)
    f(x) = score_tour(mat, x)
    update(x) = proposal_function(x, mat)
    min_val_st, x_st, p_st = PTChain.sampleTempering(start_x, f, temps, update, micro_steps, macro_steps)
    min_val_gt, x_gt, p_gt = PTChain.greedyTempering(start_x, f, temps, update, micro_steps, macro_steps)
    #p_g = plot(p_gt', title="Greedy Tempering")
    #p_s = plot(p_st', title="Parallel Tempering")
    #plot(p_g, p_s, layout=(2,1))
    #gui()
    return min_val_st, min_val_gt
end

function diffs(N)
    x = randperm(1032)
    vals = Array{Float64}(N)
    mat = load_sample()
    for i = 1:N
        x = reverse_section(x)
        vals[i] = score_tour(mat, x)
    end
    histogram(diff(vals))
    gui()
    vals
end

println(sim_tour())

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
using Stage
using Base.Test

ckpts = Checkpoints(".ckpts")
testdata = Any[ 10, "test-2", Dict("test" => 1, "test2" => 2), [1 2; 3 4] ]

logger = Log(STDERR)
for i = 1:length(testdata)
  if !haskey(ckpts, "test-$i")
    ckpts["test-$i"] = testdata[i]
  end
  @info(logger, "running test $i")
  @test haskey(ckpts, "test-$i")
  @test fetch(ckpts["test-$i"]) == testdata[i]
end


text:
export accumulated

"""
    accumulated()

Creates an `accumulated` operator, which returns an Observable that emits the current item with all of the previous items emitted by the source Observable
in one single ordered array.

# Producing

Stream of type `<: Subscribable{Vector{L}}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> accumulated(), logger())
;

# output

[LogActor] Data: [1]
[LogActor] Data: [1, 2]
[LogActor] Data: [1, 2, 3]
[LogActor] Completed

```

```jldoctest
using Rocket

source = of(1)
subscribe!(source |> accumulated(), logger())
;

# output

[LogActor] Data: [1]
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
accumulated() = AccumulatedOperator()

struct AccumulatedOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{ Vector{L} }, operator::AccumulatedOperator, source) where { L }
    return proxy(Vector{L}, source, AccumulatedProxy())
end

operator_right(::AccumulatedOperator, ::Type{L}) where L = Vector{L}

struct AccumulatedProxy <: ActorProxy end

actor_proxy!(::Type{ Vector{L} }, proxy::AccumulatedProxy, actor::A) where { L, A } = AccumulatedActor{L, A}(Vector{L}(), actor)

struct AccumulatedActor{L, A} <: Actor{L}
    values :: Vector{L}
    actor  :: A
end

on_next!(actor::AccumulatedActor{L}, data::L) where L = begin push!(actor.values, data); next!(actor.actor, actor.values) end
on_error!(actor::AccumulatedActor, err)               = error!(actor.actor, err)
on_complete!(actor::AccumulatedActor)                 = complete!(actor.actor)

Base.show(io::IO, ::AccumulatedOperator)         = print(io, "AccumulatedOperator()")
Base.show(io::IO, ::AccumulatedProxy)            = print(io, "AccumulatedProxy()")
Base.show(io::IO, ::AccumulatedActor{L}) where L = print(io, "AccumulatedActor($L)")


text:
using AxisSets: flatten

@testset "utils" begin
    @testset "flatten" begin
        @testset "Pairs{Symbol}" begin
            data = [
                :val1 => [:a1 => 1, :a2 => 2],
                :val2 => [:b1 => 11, :b2 => 22],
                :val3 => [111, 222],
                :val4 => 4.3,
            ]
            expected = [
                (:val1, :a1) => 1,
                (:val1, :a2) => 2,
                (:val2, :b1) => 11,
                (:val2, :b2) => 22,
                (:val3,) => [111, 222],
                (:val4,) => 4.3,
            ]
            @test flatten(data) == expected

            delimited = [
                :val1_a1 => 1,
                :val1_a2 => 2,
                :val2_b1 => 11,
                :val2_b2 => 22,
                :val3 => [111, 222],
                :val4 => 4.3,
            ]
            @test flatten(data, :_) == delimited
        end

        @testset "Pairs{String}" begin
            data = [
                "val1" => ["a1" => 1, "a2" => 2],
                "val2" => ["b1" => 11, "b2" => 22],
                "val3" => [111, 222],
                "val4" => 4.3,
            ]
            expected = [
                ("val1", "a1") => 1,
                ("val1", "a2") => 2,
                ("val2", "b1") => 11,
                ("val2", "b2") => 22,
                ("val3",) => [111, 222],
                ("val4",) => 4.3,
            ]
            @test flatten(data) == expected

            delimited = [
                "val1_a1" => 1,
                "val1_a2" => 2,
                "val2_b1" => 11,
                "val2_b2" => 22,
                "val3" => [111, 222],
                "val4" => 4.3,
            ]
            @test flatten(data, "_") == delimited

            # Test a nested singleton pair
            @test flatten("val1" => "a1" => 1) == [("val1", "a1") => 1]
        end

        @testset "Dict" begin
            data = Dict(
                :val1 => Dict(:a1 => 1, :a2 => 2),
                :val2 => Dict(:b1 => 11, :b2 => 22),
                :val3 => [111, 222],
                :val4 => 4.3,
            )
            expected = Dict(
                (:val1, :a1) => 1,
                (:val1, :a2) => 2,
                (:val2, :b1) => 11,
                (:val2, :b2) => 22,
                (:val3,) => [111, 222],
                (:val4,) => 4.3,
            )
            @test flatten(data) == expected

            # Test that we can operate on keys that aren't symbols or strings
            data = Dict(
                DateTime(2021, 1, 1, 11) => Dict(1 => "a", 2 => "b"),
                DateTime(2021, 1, 1, 12) => Dict(1 => "x", 2 => "y"),
                DateTime(2021, 1, 1, 13) => [111, 222],
                DateTime(2021, 1, 1, 14) => 4.3,
            )

            expected = Dict(
                (DateTime(2021, 1, 1, 11), 1) => "a",
                (DateTime(2021, 1, 1, 11), 2) => "b",
                (DateTime(2021, 1, 1, 12), 1) => "x",
                (DateTime(2021, 1, 1, 12), 2) => "y",
                (DateTime(2021, 1, 1, 13),) => [111, 222],
                (DateTime(2021, 1, 1, 14),) => 4.3,
            )
            @test flatten(data) == expected
            @test flatten(expected) == expected

            expected = Dict(
                "2021-01-01T11:00:00_1" => "a",
                "2021-01-01T11:00:00_2" => "b",
                "2021-01-01T12:00:00_1" => "x",
                "2021-01-01T12:00:00_2" => "y",
                "2021-01-01T13:00:00" => [111, 222],
                "2021-01-01T14:00:00" => 4.3,
            )
            # Test that we can flatten that to a string
            @test flatten(data, "_") == expected
            @test flatten(expected, "_") == expected
        end

        @testset "NamedTuple" begin
            data = (
                val1 = (a1 = 1, a2 = 2),
                val2 = (b1 = 11, b2 = 22),
                val3 = [111, 222],
                val4 = 4.3,
            )
            # Expected names for NamedTuple must be symbols. Default delimiter is `:_`
            expected = (
                val1_a1 = 1,
                val1_a2 = 2,
                val2_b1 = 11,
                val2_b2 = 22,
                val3 = [111, 222],
                val4 = 4.3,
            )
            @test flatten(data, :_) == expected
            @test flatten(expected, :_) == expected

            # By default calling `flatten` on a named tuple without a delimiter will
            # return pairs
            expected = [
                (:val1, :a1) => 1,
                (:val1, :a2) => 2,
                (:val2, :b1) => 11,
                (:val2, :b2) => 22,
                (:val3,) => [111, 222],
                (:val4,) => 4.3,
            ]
            @test flatten(data) == expected
            @test flatten(expected) == expected
        end

        @testset "KeyedArray" begin
            dt = DateTime(2021, 1, 1, 11):Hour(1):DateTime(2021, 1, 1, 14)
            A = KeyedArray(
                reshape(1:24, (4, 3, 2));
                time=dt,
                loc=1:3,
                obj=[:a, :b],
            )

            @testset "tail dims" begin
                expected = KeyedArray(
                    reshape(1:24, (4, 6));
                    time=dt,
                    loc_obj=[Symbol(join((l, o), :_)) for l in 1:3, o in [:a, :b]][:],
                )
                @test flatten(A, (:loc, :obj), :_) == expected
            end

            @testset "head dims" begin
                expected = KeyedArray(
                    reshape(1:24, (12, 2));
                    time_loc=[Symbol(join((t, l), :_)) for t in dt, l in 1:3][:],
                    obj=[:a, :b],
                )
                @test flatten(A, (:time, :loc), :_) == expected
            end

            @test_throws ArgumentError flatten(A, (:time,), :_)
            @test_throws ArgumentError flatten(A, (:time, :obj), :_)
        end

        # TODO: Support group flatten operations
    end
end


text:
using Revise
using NNFEM
using PyPlot
using LinearAlgebra
using ADCME

NT = 200
Δt = 1/NT

n = 10
m = 2n
h = 1/n

# Create a very simple mesh
elements = SmallStrainContinuum[]
prop = Dict("name"=> "PlaneStrain", "rho"=> 1.0, "E"=>2.0, "nu"=>0.35)
coords = zeros((m+1)*(n+1), 2)
for j = 1:n
    for i = 1:m
        idx = (m+1)*(j-1)+i
        elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
        ngp = 3
        nodes = [
            (i-1)*h (j-1)*h
            i*h (j-1)*h
            i*h j*h
            (i-1)*h j*h
        ]
        coords[elnodes, :] = nodes
        push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
    end
end

Edge_Traction_Data = Array{Int64}[]
for i = 1:m
  elem = elements[i]
  for k = 1:4
    if elem.coords[k,2]<0.001 && elem.coords[k+1>4 ? 1 : k+1,2]<0.001
      push!(Edge_Traction_Data, [i, k, 1])
    end
  end
end

for i = 1:n
  elem = elements[(i-1)*m+1]
  for k = 1:4
    if elem.coords[k,1]<0.001 && elem.coords[k+1>4 ? 1 : k+1,1]<0.001
      push!(Edge_Traction_Data, [(i-1)*m+1, k, 0])
    end
  end
end

Edge_Traction_Data = hcat(Edge_Traction_Data...)'|>Array

# fixed on the bottom, push on the right
EBC = zeros(Int64, (m+1)*(n+1), 2)
FBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
f = zeros((m+1)*(n+1), 2)
for j = 1:n+1
    idx = (j-1)*(m+1) + m+1
    EBC[idx,:] .= -2 # time-dependent boundary, right
end
for i = 1:m+1
    idx = n*(m+1) + i
    EBC[idx,:] .= -1 # fixed boundary, bottoms
end

dimension = 2
domain = Domain(coords, elements, dimension, EBC, g, FBC, f, Edge_Traction_Data)

xy = domain.nodes[domain.dof_to_eq]
n_ = div(length(xy),2)
x = xy[1:n_,1]
y = xy[n_+1:end,1]
# Set initial condition
Dstate = zeros(domain.neqs) # d at last step
state = [(@. (1-y^2)*(x^2+y^2)); (@. (1-y^2)*(x^2-y^2))] * 0.1
velo = -[(@. (1-y^2)*(x^2+y^2)); (@. (1-y^2)*(x^2-y^2))] * 0.1
acce = [(@. (1-y^2)*(x^2+y^2)); (@. (1-y^2)*(x^2-y^2))] * 0.1
gt = nothing
ft = nothing

EBC_DOF = findall(EBC[:,1] .== -2)
x_EBC = domain.nodes[EBC_DOF,1]
y_EBC = domain.nodes[EBC_DOF,2]
function EBC_func(t)
  out = [(@. 0.1*(1-y_EBC^2)*(x_EBC^2+y_EBC^2)*exp(-t));
    (@. 0.1*(1-y_EBC^2)*(x_EBC^2-y_EBC^2)*exp(-t))]

  out, -out, out
end

function Body_func_linear_elasticity(x, y, t)
    f1 = @. 0.987654320987654*x*y*exp(-t) + 0.592592592592593*y^2*exp(-t) + (0.1 - 0.1*y^2)*(x^2 + y^2)*exp(-t) - (0.148148148148148 - 0.148148148148148*y^2)*exp(-t) - (0.641975308641975 - 0.641975308641975*y^2)*exp(-t) - (-0.148148148148148*x^2 - 0.148148148148148*y^2)*exp(-t)
    f2 = @. 0.987654320987654*x*y*exp(-t) - 2.5679012345679*y^2*exp(-t) + (0.1 - 0.1*y^2)*(x^2 - y^2)*exp(-t) - (0.148148148148148 - 0.148148148148148*y^2)*exp(-t) - (-0.641975308641975*x^2 + 0.641975308641975*y^2)*exp(-t) - (0.641975308641975*y^2 - 0.641975308641975)*exp(-t)
    [f1 f2]
end

FBC_func = nothing
Edge_func = nothing
globaldata = GlobalData(state, Dstate, velo, acce, domain.neqs, EBC_func, FBC_func,Body_func_linear_elasticity, Edge_func)

x = domain.nodes[:,1]
y = domain.nodes[:,2]
d0 = [(@. (1-y^2)*(x^2+y^2)); (@. (1-y^2)*(x^2-y^2))] * 0.1
v0 = -[(@. (1-y^2)*(x^2+y^2)); (@. (1-y^2)*(x^2-y^2))] * 0.1
a0 = [(@. (1-y^2)*(x^2+y^2)); (@. (1-y^2)*(x^2-y^2))] * 0.1
assembleMassMatrix!(globaldata, domain)

# linear elasticity matrix at each Gauss point
Hs = zeros(domain.neles*length(domain.elements[1].weights), 3, 3)
for i = 1:size(Hs,1)
    Hs[i,:,:] = elements[1].mat[1].H
end

# Construct Edge_func
function Edge_func_linear_elasticity(x, y, t, idx)
  if idx==0
      f1 = @. -6.41975308641975*x*(0.1 - 0.1*y^2)*exp(-t) + 3.45679012345679*y*(0.1 - 0.1*y^2)*exp(-t) + 0.345679012345679*y*(x^2 - y^2)*exp(-t)
      f2 = @. -1.48148148148148*x*(0.1 - 0.1*y^2)*exp(-t) - 1.48148148148148*y*(0.1 - 0.1*y^2)*exp(-t) + 0.148148148148148*y*(x^2 + y^2)*exp(-t)
    elseif idx==1
      f1 = @. -1.48148148148148*x*(0.1 - 0.1*y^2)*exp(-t) - 1.48148148148148*y*(0.1 - 0.1*y^2)*exp(-t) + 0.148148148148148*y*(x^2 + y^2)*exp(-t)
      f2 = @. -3.45679012345679*x*(0.1 - 0.1*y^2)*exp(-t) + 6.41975308641975*y*(0.1 - 0.1*y^2)*exp(-t) + 0.641975308641975*y*(x^2 - y^2)*exp(-t)
    end
    return [f1 f2]
end
globaldata.Edge_func = Edge_func_linear_elasticity

ts = ExplicitSolverTime(Δt, NT)
ubd, abd = compute_boundary_info(domain, globaldata, ts)
Fext = compute_external_force(domain, globaldata, ts)
d, v, a= ExplicitSolver(globaldata, domain, d0, v0, a0, Δt, NT, Hs, Fext, ubd, abd)

# # NOTE: You can also use the implicit alpha solvers
# ts = GeneralizedAlphaSolverTime(Δt, NT)
# ubd, abd = compute_boundary_info(domain, globaldata, ts)
# Fext = compute_external_force(domain, globaldata, ts)
# d, v, a= GeneralizedAlphaSolver(globaldata, domain, d0, v0, a0, Δt, NT, Hs, Fext, ubd, abd)

sess = Session(); init(sess)
d_, v_, a_ = run(sess, [d,v,a])

using Random; Random.seed!(233)
for k = 1:5
    i = rand(1:m+1)
    j = rand(1:n+1)
    if k==1
        plot(d_[:,(j-1)*(m+1)+i], color = "C$k", label="Computed")
    else
        plot(d_[:,(j-1)*(m+1)+i], color = "C$k")
    end
    x0 = (i-1)*h
    y0 = (j-1)*h

    if k==1
        plot((@. (1-y0^2)*(x0^2+y0^2)*exp(-ts))*0.1 ,"--", color="C$k", label="Reference")
    else
        plot((@. (1-y0^2)*(x0^2+y0^2)*exp(-ts))*0.1 ,"--", color="C$k")
    end
end
legend()


text:
function fdjac(y0::Vector{Float64},t0::Float64,tf::Float64,sparams::WK3.SolverParams,
    mparams::WK3.ModelParams,obs::Vector{Float64})
    sqrteps = sqrt(eps()); # ≈ square root of machine precision
    nparams = 9;
    nobs = 1;
    J = zeros(nobs,nparams);
    for j = 1:nparams
        if j == 1
            temp = mparams.Zc;
        elseif j == 2
            temp = mparams.R;
        elseif j == 3
            temp = mparams.C;
        elseif j == 4
            temp = mparams.V0;
        elseif j == 5
            temp = mparams.t1;
        elseif j == 6
            temp = mparams.t2;
        elseif j == 7
            temp = mparams.m1;
        elseif j == 8
            temp = mparams.m2;
        elseif j == 9
            temp = mparams.Emax;
        end
        h = sqrteps*abs.(temp); # scale step size by current value
        if h == 0
            h = sqrteps;
        end
        if j == 1
            mparams.Zc = temp+h; # trick to reduce finite precision error
            h = mparams.Zc - temp;
        elseif j == 2
            mparams.R = temp+h;
            h = mparams.R - temp;
        elseif j == 3
            mparams.C = temp+h;
            h = mparams.C - temp;
        elseif j == 4
            mparams.V0 = temp+h;
            h = mparams.V0 - temp;
        elseif j == 5
            mparams.t1 = temp+h;
            h = mparams.t1 - temp;
        elseif j == 6
            mparams.t2 = temp+h;
            h = mparams.t2 - temp;
        elseif j == 7
            mparams.m1 = temp+h;
            h = mparams.m1 - temp;
        elseif j == 8
            mparams.m2 = temp+h;
            h = mparams.m2 - temp;
        elseif j == 9
            mparams.Emax = temp+h;
            h = mparams.Emax - temp;
        end
        ~,yout = odeint(y0,t0,tf,sparams,mparams);
        if j == 1
            mparams.Zc = temp;
        elseif j == 2
            mparams.R = temp;
        elseif j == 3
            mparams.C = temp;
        elseif j == 4
            mparams.V0 = temp;
        elseif j == 5
            mparams.t1 = temp;
        elseif j == 6
            mparams.t2 = temp;
        elseif j == 7
            mparams.m1 = temp;
        elseif j == 8
            mparams.m2 = temp;
        elseif j == 9
            mparams.Emax = temp;
        end
        obsp = [yout[end][3]]; # dimensionless flowrate from perturbed parameter
        for i = 1:nobs # forward difference formula
            J[i,j] = (obsp[i]-obs[i])/h;
        end
    end
    return J
end


text:
"""
    This example shows how to use Term's introspection functionality.

Introspection refers to a piece of code analyzing another piece of code.
Say you want to know what a custom data type is like, or what methods does
a function have etc... What would you do? You can have a look at the docs, you
can go on the GitHub repository and look around, or yo can use Term!

Term provides an `inspect` method that provides a lot of useful information!
"""

import Term: inspect, typestree

# define some types
abstract type T1 end

abstract type T2 <: T1 end

"""
    MyType

It's just a useless type we've defined to provide an example of
Term's awesome `inspect` functionality!
"""
struct MyType <: T2
    x::Int
    y::String
    fn::Any
end

#  Now, what is MyType like?
inspect(MyType)

# Let's define some constructors and methods using MyType

"""
constructors!
"""
MyType(x::Int) = MyType(x, "no string", nothing)

MyType(x::Int, y) = MyType(x, y, nothing)

# methods
useless_method(m::MyType) = m

"""
these methods don't do much, just an example
"""
another_method(m1::MyType, m2::MyType) = print(m1, m2)

# let's inspect MyType again!
inspect(MyType)

"""
Not bad huh !?
It also works with instances of types
"""

inspect(MyType(2))

"""
And for abstract types of course1
"""

inspect(T2)

"""
finally, you can also use inspect to look at functions.
Let's get a bit meta
"""

inspect(inspect)

"""
as you can see inspect prints the docstrings of the methods it can find and lists some of
the methods available.
"""

"""
If you just want to see the type's hierarchy for a type, you can use typestree
"""

typestree(AbstractFloat)


text:
# this runs in your environment and assumes that all your scripts are up to date etc.
# do not run this if you have not previously run `serve()` and checked everything was fine.

Pkg.activate()

using JuDoc, Pkg, Logging, NodeJS, Literate

##################################################################
## PART 1: generate all scripts and notebooks and push everything
# on the gh-pages branch (yes we're abusing the system)
##################################################################

Logging.disable_logging(Logging.LogLevel(1500))

ifiles = joinpath.("scripts", readdir("scripts"))

nbpath = joinpath("generated", "notebooks")
scpath = joinpath("generated", "scripts")

isdir(nbpath) || mkpath(nbpath)
isdir(scpath) || mkpath(scpath)

ACTIVATE = """
    # Before running this, make sure to activate the environment corresponding to
    # [this `Project.toml`](https://raw.githubusercontent.com/tlienart/MLJTutorials/master/Project.toml)
    # and update it so that you get an environment which matches the one used to generate
    # the tutorials:
    #
    # ```julia
    # cd("MLJTutorials") # cd to folder with the Project.toml
    # using Pkg
    # Pkg.activate(".")
    # Pkg.update()
    # ```

    """

preproc(s) = ACTIVATE * s

Literate.notebook.(ifiles, nbpath, preprocess=preproc, execute=false, documenter=false)
Literate.script.(ifiles, scpath, postprocess=preproc, keep_comments=false, documenter=false)

JS_GHP = """
    var ghpages = require('gh-pages');
    ghpages.publish('generated/', function(err) {});
    """

run(`$(nodejs_cmd()) -e $JS_GHP`)

##################################################################
## PART 2
# optimize (compress and pre-render) the website and publish it
# to the master branch this time.
##################################################################

publish()

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
# Modelo simplificado gerador de vapor
type Gerador_Vapor_Simples
	Gerador_Vapor_Simples()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Potencia(),
			Eficiencia(),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}"
			)),
			[
				:(Fout.P = Fin.P),
				:([Fout.S,Fout.H] = PP2.propPTv(Fout.P,Fout.T)),
				:(Q_GV * EF_GV = Fin.F * (Fout.H - Fin.H)),
			],
			[
				"","","",
			],
			[:PP2,],
			[:Q_GV,:EF_GV,:Fin,:Fout,]
		)
	end
	PP2::DanaPlugin
	Q_GV::Potencia
	EF_GV::Eficiencia
	Fin::Corrente
	Fout::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
end
export Gerador_Vapor_Simples
function setEquationFlow(in::Gerador_Vapor_Simples)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	#	Fout.F = Fin.F;

end


text:
export BitFlippingEnv

"""
In Bit Flipping Environment we have n bits. The actions are 1 to n where executing i-th action flips the i-th bit of the state. For every episode we sample uniformly and initial state as well as the target state.

Refer [Hindsight Experience Replay paper](https://arxiv.org/pdf/1707.01495.pdf) for the motivation behind the environment.
"""

mutable struct BitFlippingEnv <: AbstractEnv
    N::Int
    rng::AbstractRNG
    state::BitArray{1}
    goal_state::BitArray{1}
    max_steps::Int
    t::Int
end

function BitFlippingEnv(; N = 8, T = N, rng = Random.GLOBAL_RNG)
    state = bitrand(rng, N)
    goal_state = bitrand(rng, N)
    max_steps = T
    BitFlippingEnv(N, rng, state, goal_state, max_steps, 0)
end

Random.seed!(env::BitFlippingEnv, s) = Random.seed!(env.rng, s)

RLBase.action_space(env::BitFlippingEnv) = Space(OneTo(env.N))

RLBase.legal_action_space(env::BitFlippingEnv) = Base.OneTo(env.N)

function (env::BitFlippingEnv)(action::Int)
    env.t += 1
    if 1 <= action <= env.N
        env.state[action] = !env.state[action]
        nothing
    else
        @error "Invalid Action"
    end
end

RLBase.state(env::BitFlippingEnv) = state(env::BitFlippingEnv, Observation{BitArray{1}}())
RLBase.state(env::BitFlippingEnv, ::Observation) = env.state
RLBase.state(env::BitFlippingEnv, ::GoalState) = env.goal_state
RLBase.state_space(env::BitFlippingEnv, ::Observation) = Space(Bool, env.N)
RLBase.state_space(env::BitFlippingEnv, ::GoalState) = Space(Bool, env.N)
RLBase.is_terminated(env::BitFlippingEnv) =
    (env.state == env.goal_state) || (env.t >= env.max_steps)

function RLBase.reset!(env::BitFlippingEnv)
    env.t = 0
    env.state .= bitrand(env.rng, env.N)
    env.goal_state .= bitrand(env.rng, env.N)
end

function RLBase.reward(env::BitFlippingEnv)
    if env.state == env.goal_state
        0.0
    else
        -1.0
    end
end

RLBase.NumAgentStyle(::BitFlippingEnv) = SINGLE_AGENT
RLBase.DynamicStyle(::BitFlippingEnv) = SEQUENTIAL
RLBase.ActionStyle(::BitFlippingEnv) = MINIMAL_ACTION_SET
RLBase.InformationStyle(::BitFlippingEnv) = PERFECT_INFORMATION
RLBase.StateStyle(::BitFlippingEnv) = (Observation{BitArray{1}}(), GoalState{BitArray{1}}())
RLBase.RewardStyle(::BitFlippingEnv) = STEP_REWARD
RLBase.UtilityStyle(::BitFlippingEnv) = GENERAL_SUM
RLBase.ChanceStyle(::BitFlippingEnv) = STOCHASTIC


text:
# Run the model for the ensemble of solutions -
include("SolveBalances.jl")
using PyPlot
using PyCall
@pyimport numpy as np

# Load the ensemble files from disk -
pc_array_full = readdlm("./data/pc_array_O1_O2.dat")
rank_array = readdlm("./data/rank_array_O1_O2.dat")

# Select the desired rank -
idx_rank = find(rank_array .<= 5.0)

# Setup time scale -
tStart = 0.0;
tStop = 25.0;
tStep = 0.1;
number_of_timesteps = 200
time_experimental = linspace(tStart,tStop,number_of_timesteps)

# initialize data array -
data_array = zeros(length(time_experimental),1)

data_dictionary = Dict()
initial_condition_array = zeros(19)
initial_condition_array[1] = 1.0      # 1 zymosan
initial_condition_array[2] = 1.9      # 2 C4
initial_condition_array[3] = 0.322    # 3 C2
initial_condition_array[8] = 7.57     # 4 C3
initial_condition_array[14] = 0.195   # 5 C5
initial_condition_array[17] = 2.23    # 6 Factor H
initial_condition_array[18] = 0.417   # 7 C4BP

initial_condition_array[19] = 0.6    # 19 C3a (ic initial)
initial_condition_array[15] = 0.00006   # 19 C5a (ic initial)

data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
bad_set = Set()

# Run the simulation for this rank -
number_of_samples = 1
for (index,rank_index_value) in enumerate(idx_rank)

  if (mod(index,1)==0)

    # Grab the parameter set from the cache -
    parameter_array = pc_array_full[:,rank_index_value]
    data_dictionary["PARAMETER_ARRAY"] = parameter_array

    # Run the model -
    (t,x) = SolveBalances(tStart,tStop,tStep,data_dictionary)

    # @show x[1,1]

    # Need to interpolate the simulation onto the experimental time scale -
    IC3a = np.interp(time_experimental[:],t,x[:,19])
    IC5a = np.interp(time_experimental[:],t,x[:,15])

    # grab -
    data_array = [data_array IC5a]

    # How many samples?
    number_of_samples = number_of_samples + 1

    @show index
  end
end

# calculate mean, and std
mean_value = mean(data_array,2)
std_value = std(data_array,2)

# make best fit -
SF = 2.58
UB = mean_value + (SF)*std_value
LB = mean_value - (SF)*std_value
idx_z = find(LB.<0)
LB[idx_z] = 0.0
fill_between(vec(time_experimental),vec(LB),vec(UB),color=[0.8,0.8,0.8],lw=2)

# 99% of mean -
SF = (2.58/sqrt(number_of_samples))
UB = mean_value + (SF)*std_value
LB = mean_value - (SF)*std_value
idx_z = find(LB.<0)
LB[idx_z] = 0.0

# Make the plot -
plot(time_experimental,mean_value,"k-")
fill_between(vec(time_experimental),vec(LB),vec(UB),color="gray",lw=2)

# Load the data -
# 1.0 mg/ml -
MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2e_C3a.txt",',')
MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3c_C5a_original.txt",',')

# 0.1 mg/ml -
# MEASURED_ARRAY_C3a = readdlm("./data/Shaw2015_Fig2d_C3a.txt",',')
# MEASURED_ARRAY_C5a = readdlm("./data/Shaw2015_Fig3b_C5a_original.txt",',')

MC3a = 0.11002*MEASURED_ARRAY_C3a[:,2]
MC5a = 0.00009615*MEASURED_ARRAY_C5a[:,2]
plot(MEASURED_ARRAY_C3a[:,1],MC5a,color="black",marker="o",markersize=8.0)


text:
using MatrixCompletion.Concepts
import LinearAlgebra

struct RelativeError end
struct AbsoluteError end

function within_radius(x;of=1e-5,at=0)
  return sum(Int.(abs.(x .- at) .> of))
end

function Concepts.provide(object::RelativeError,x,y;
                          metric::Any = x -> LinearAlgebra.norm(x,2),
                          base_metric::Any=metric)
  return metric(x-y)/base_metric(y)
end

function Concepts.provide(object::AbsoluteError,x,y;
                          metric::Any= x-> LinearAlgebra.norm(x,2))
  return metric(x - y)
end

function Concepts.provide(object::Diagnostics{Int64})
  return "int64"
end

function Concepts.provide(object::Diagnostics{<:Any};
                          reference::Optional{VecOrMatOfReals}=nothing,
                          input_data::Optional{VecOrMatOfReals}=nothing)
  if isnothing(input_data) && isnothing(reference)
    throw(DomainError("[provide(Diagnostics)]: input_data and/or reference variable missing))"))
  end
  return Dict("relative-error[#within-radius(1e-5)]" => Concepts.provide(RelativeError(),
                                                                         input_data,reference,
                                                                         metric      = x -> within_radius(x),
                                                                         base_metric = x -> LinearAlgebra.norm(abs.(x) .+ 1 ,0)),
              "absolute-error[#within-radius(1e-5)]" => Concepts.provide(AbsoluteError(),
                                                                         input_data, reference,
                                                                         metric      = x -> within_radius(x)),
              "relative-error[#within-radius(1)]" => Concepts.provide(RelativeError(),
                                                                      input_data,reference,
                                                                      metric      = x -> within_radius(x,of=1),
                                                                      base_metric = x -> LinearAlgebra.norm(abs.(x) .+ 1 ,0)),
              "absolute-error[#within-radius(1)]" => Concepts.provide(AbsoluteError(),
                                                                      input_data, reference,
                                                                      metric      = x -> within_radius(x,of=1)),
              "relative-error[L1]" => Concepts.provide(RelativeError(),
                                                       input_data,reference,
                                                       metric = x -> LinearAlgebra.norm(x,1)),
              "absolute-error[L1]" => Concepts.provide(AbsoluteError(),
                                                       input_data, reference,
                                                       metric = x -> LinearAlgebra.norm(x,1)),
              "relative-error[L2]" => Concepts.provide(RelativeError(),
                                                       input_data,reference,
                                                       metric = x -> LinearAlgebra.norm(x,2)^2),
              "absolute-error[L2]" => Concepts.provide(AbsoluteError(),
                                                       input_data,reference,
                                                       metric = x -> LinearAlgebra.norm(x,2)^2)
              )

end

# function print_diagnostics(object::Diagnostics{<:Any};
#                            reference::Optional{VecOrMatOfReals}=nothing,
#                            input_data::Optional{VecOrMatOfReals}=nothing)

# end


text:
using CSFML, CSFML.LibCSFML

import Base.Threads.Atomic
import Base.Threads.atomic_add!

mutable struct AudioStream
  header::Vector{UInt8}
  data::Channel{UInt8}
  position::Atomic{Int64}
end

AudioStream() = AudioStream(header(typemax(UInt32) - headerSize), Channel{UInt8}(Inf), Atomic{Int64}(1))

const headerSize = UInt32(44)

function header(datalength::UInt32)
  h = zeros(UInt8, headerSize)
  h[1:4] = Vector{UInt8}("RIFF")
  h[5:8] = reinterpret(UInt8, [headerSize + datalength])
  h[9:12] = Vector{UInt8}("WAVE")
  h[13:16] = Vector{UInt8}("fmt ")
  h[17:20] = reinterpret(UInt8, [UInt32(16)])
  h[21] = 0x01
  h[23] = 0x01
  h[25:28] = reinterpret(UInt8, [UInt32(44100)])
  h[29:32] = reinterpret(UInt8, [UInt32(44100 * 2)])
  h[33] = 0x02
  h[35] = 0x10
  h[37:40] = Vector{UInt8}("data")
  h[41:44] = reinterpret(UInt8, [datalength])
  h
end

function read!(dest::Vector{UInt8}, as::AudioStream, num::Int64)
  position = atomic_add!(as.position, num)
  for i = 1:num
    dest[i] = if position <= length(as.header)
      as.header[position]
    else
      take!(as.data)
    end
    position += 1
  end
  num
end

function readc!(data::Ptr{Cvoid}, size::Int64, userdata::Ptr{Cvoid})
  data = unsafe_wrap(Vector{UInt8}, Base.unsafe_convert(Ptr{UInt8}, data), size)
  as = unsafe_pointer_to_objref(userdata)
  read!(data, as, size)
end

function seek!(as::AudioStream, position::Int64)
  as.position[] = position
end

function seekc!(position::Clonglong, userdata::Ptr{Cvoid})
  as = unsafe_pointer_to_objref(userdata)
  seek!(as, position + 1) - 1
end

tell(as::AudioStream) = as.position[]

function tellc(userdata::Ptr{Cvoid})
  as = unsafe_pointer_to_objref(userdata)
  tell(as) - 1
end

getsize(as::AudioStream) = Int64(typemax(UInt32))

function getsizec(userdata::Ptr{Cvoid})
  as = unsafe_pointer_to_objref(userdata)
  getsize(as)
end

function inputstream(as::AudioStream)
  readptr = @cfunction(readc!, Clonglong, (Ptr{Cvoid}, Clonglong, Ptr{Cvoid}))
  seekptr = @cfunction(seekc!, Clonglong, (Clonglong, Ptr{Cvoid}))
  tellptr = @cfunction(tellc, Clonglong, (Ptr{Cvoid},))
  getsizeptr = @cfunction(getsizec, Clonglong, (Ptr{Cvoid},))
  sfInputStream(readptr, seekptr, tellptr, getsizeptr, pointer_from_objref(as))
end

music(as::AudioStream) = sfMusic_createFromStream(Ref(inputstream(as)))


text:
function doSort!(a,β = 1.0)
    n = length(a)
    k = min(Int(round(β*length(a))),n-1)
    @inbounds for i in 1:k-1
        change = false
        @inbounds for j in 1:n-i
            if a[j] > a[j+1]
                a[j],a[j+1] = a[j+1],a[j]
                change = true
            end
        end
        !change && break
    end
    return a
end

envelope(data) = accumulate(max,data)

B_curve(x, β, dist) = min(β/ccdf(dist,x),1) - min(β,ccdf(dist,x))
B_curve_normal(x,β) = B_curve(x,β,Normal())
B_curve_uniform(x,β) = B_curve(x,β,Uniform())
B_curve_wide_uniform(x,β) = B_curve(x,β,Uniform(-1,1))

function bs_stat(data, β, dist)
    n = length(data)
    sdata = copy(data)
    doSort!(sdata, β)
    envFunction = ecdf(envelope(sdata))
    running_max_data = envelope(sdata)

    m = max( abs(B_curve(running_max_data[1],β,dist)) , abs(envFunction(running_max_data[1]) - B_curve(running_max_data[1], β, dist)))
    for i in 2:n
        a_diff = abs(envFunction(running_max_data[i]) - B_curve(running_max_data[i], β, dist))
        (a_diff > m) && (m = a_diff)
        a_diff = abs(envFunction(running_max_data[i-1]) - B_curve(running_max_data[i], β, dist))
        (a_diff > m) && (m = a_diff)
    end
    (√n)*m
end

bs_stat_normal(data, β) = bs_stat(data, β, Normal())
bs_stat_uniform(data, β) = bs_stat(data, β, Uniform())


text:
export getCreateLastRow, Page, Scroller, Row, Element, TextElement,
       feedRow #Event,

#begin

# if this structure doesn't grow it should be elimenated. Pointless to have a structure for one member...
    mutable struct EventType
        pressed::Vector{Point}
        selectedNodes::Any
        EventType() = new([],[])
        EventType(x,y) = new([Point(x,y)],[])
        EventType(a) = new(a)
    end
#-==============================================================================
mutable struct Scroller
    x::Float64
    y::Float64
    contentWidth::Float64
    contentHeight::Float64
    Scroller() = new(0,0,0,0)
end
#-==============================================================================
mutable struct Row
    flags::BitArray{1} #Any
    nodes::Vector{Any} # cant't say Element because it's not yet defined
    space::Float64 # Space remaining 'til full
    height::Float64
    width::Float64 # because elements could be wider than the parent
    left::Float64  # TODO: change these to left/top for cmpatability
    top::Float64

    Finalize(n) = length(n.rows) > 0 ? FinalizeRow(n.rows[end]) : n.shape.top

    Row(left, top, width) = new(falses(32),[],width,0,width,left,top)
    function Row(node)
        top = Finalize(node)
        r = Row(node.shape.left, top, node.shape.width)
        push!(node.rows, r ) # There were no rows. Create one!
        return r
    end

    function Row(rows::Vector{Row}, left::Float64, top::Float64, width::Float64)
        r = Row(left, top, width)
        push!(rows, r)
        return r
    end
end
# maybe rename: GetNewRow or something similar.
getCreateLastRow(n) = length(n.rows) < 1 ? Row(n) : n.rows[end]

function getCreateLastRow(rows::Vector{Row}, l, t, w)
    if length(rows) < 1
        r = Row(l, t, w)
        push!(rows, r) # WAS: Row(rows, l, t, w)
    end
    return rows[end]
end

#-==============================================================================
#-==============================================================================
function feedRow(parent, node)

     width, height = getSize(node.shape)

    if parent.shape.width < 1 # set parent width
        parent.shape.width = width
    end

    row = getCreateLastRow(parent)

    if node.shape.flags[Marked]
        println("parrent.row.space: $(row.space) , node.width: $(width)")

    end

    if row.space < width # Not enough room for node!
        row = Row(parent) # This automatically finalises the row and returns a new one.
    end

    row.space -= width
    row.height = max(row.height, height)

    OffsetX, OffsetY = contentOffset( node.shape )
    node.shape.left = row.left + OffsetX
    row.left += width
    node.shape.top = row.top + OffsetY

    push!(row.nodes, node)
end
#-==============================================================================
mutable struct Element
    DOM::Dict       # Reference to dictionary counterpart of this node
    parent::Any               # This node's parent
    children::Vector{Any} # Children in order they appear in DOM
    rows::Vector{Row} # A layout property
    font::Any  # Font definition to be applied to all rows and children unless overriden.
    shape::Any # link to layout representation of node
    temp::Any  # Temporary copy of node to preserve data for when visual changes take place
    scroll::Scroller # Since shape get's destroyed we need the scroll-offsets here (prabably should be Nullable).
    Element(DOM) = new(DOM, nothing, [], [], nothing, nothing, nothing, Scroller())
    Element() = new(Dict(), nothing, [], [], nothing, nothing, nothing, Scroller())
end
#-==============================================================================
mutable struct TextElement
    DOM::Dict       # Reference to dictionary counterpart of this node
    parent::Any               # This node's parent
    shape::Any # link to layout representation of node
    TextElement(parent, shape) = new(Dict(), parent, shape)
    TextElement(shape) = new(Dict(), nothing, shape)
    TextElement() = new(Dict(), nothing, nothing)
end
#-==============================================================================
mutable struct AttachedEvents
    click::Vector{Element}
    mousedown::Vector{Any}
    mousemove::Vector{Any}
    mouseover::Vector{Any}
    mouseup::Vector{Any}
    mouseout::Vector{Any}

    drag::Vector{Any}
    drop::Vector{Any}

    keydown::Vector{Any}
    keypress::Vector{Any}
    keyup::Vector{Any}
    AttachedEvents() = new([], [], [], [], [], [], [], [], [], [], [])
    # onload::Array
end
# document.event
# mutable struct Event
#     node::Element
#     copy::Element # save a copy of the origional node when overwritten by new data
#     action::AbstractString # an event should generally result in an action!
#     Event(node, action) = new(node, nothing, action)
# end
#-==============================================================================
mutable struct Page
         parent::Any  # First node needs a Psudo-parent too ..maybe!
         children::Vector{Any} # First node in a tree-like data structure representing all elements on page
         fixed::Element # First node in a tree-like data structure representing all elements on page
         styles::Dict
         head::Dict
         width::Float64
         height::Float64

         win::Any
         url::Any         # URL of page
         flags::BitArray{1}        #  buttonPressed
         focusNode::Any
         canvas::Any
         event::EventType
         eventsList::AttachedEvents
             function Page(url::String)
                     children::Vector{Element} = [Element()]
                     parent = children[1]
                     children[1].parent = parent
                 new(parent, children, Element(), Dict(), Dict(), 0,0, nothing, url, falses(8),
                    nothing, 0, EventType(), AttachedEvents())
             end
end

# end


text:
"""
    Amplification(target, current, change_limit)

Apply amplification to the signal.

This modifier allows the user to specify a `target` linear amplification
value that will be applied to the signal.
The modifier will then change the amplification of the
signal until the desired amplification is achieved. The rate
at which the amplification can be changed per frame is parameterised by
the `change_limit` parameter.

To slowly ramp a signal to a desired value set the `target` amplification
to the desired value, and the `change_limit` to a small value.

To instantly change the signal set the `change_limit` to infinity and
modify the `target` value.

When initialising the modifier specify the desired starting point
using the `current` parameter.

You can access the exact amplification at any time by querying the
`current` parameter.

Inputs
------
* `target` desired linear amplification factor to be applied to signal.
* `current` linear amplification currently applied to signal.
  Also used to specify the intial value for the process.
* `change_limit` maximum change that can occur per frame.
* `enable` enable the modifier, if false the signal will be passed through without modification.

Example
-------
```julia
amplify = Amplification(0.1, 0.0, 0.05)
attenuated_sound = modify(amplify, original_sound)
```
"""
@with_kw mutable struct Amplification
    target::Float64=0             # Desired scaling factor
    current::Float64=0            # Current scaling factor
    change_limit::Float64=Inf     # Maximum change in scaling per frame
    enable::Bool=true             # Enable or pass through the signal

    Amplification(a, b, c, d) = new(a, b, c, d)
    Amplification(a, b, c) = new(a, b, c, true)
    Amplification(a, b) = new(a, b, Inf, true)
    Amplification(a) = new(a, 0, Inf, true)
end

function modify(sink::Amplification, buf)
    if sink.enable
        # Determine if the currect scaling needs to be updated
        if sink.target != sink.current
            error = sink.target - sink.current
            error_sign = sign(error)
            abs_error = abs(error)
            sink.current = sink.current + (min(sink.change_limit, abs_error) * error_sign)
        end
        buf.data = buf.data .* sink.current
    end
    return buf
end


text:
module Huggett890

using DocStringExtensions, Infiltrator, Interpolations, Lazy, Roots
using UtilityFunctions890

export init_demographics, validate_demographics
export init_test_endowments, validate_endowments
export init_test_capital_grid, validate_capital_grid, make_k_grid
export init_test_budget, validate_budget, income, kprime, consumption
export Model, init_model, validate_model, solve_worker, euler_dev
export WorkerSolution, init_test_worker_solution
export Solution, init_solution, init_test_solution

# The objects
include("types.jl");
include("demographics.jl");
include("endowments.jl");
include("capital_grid.jl");
include("budget.jl");
include("model.jl");
include("solution.jl");
# Solving the model
# include("solve.jl");
include("solve_worker.jl");
include("helpers.jl");

"""
	$(SIGNATURES)

Runs everything in sequence
"""
function run_all()
    m = init_model();
    sol = solve_worker(m);
    save_solution(sol);
    sim = simulate_model(m, sol);
    show_results(m, sol, sim);
end

# stub +++
function save_solution(sol) end

# stub +++
function simulate_model(m, sol) end

# stub +++
function show_results(m, sol, sim) end

end # module


text:
using FileIO
using DelimitedFiles

# Load the dataset and create a matrix with the inputs and a vector for the targets
dataset = DelimitedFiles.readdlm("datasets/561_cpu.tsv");
inputs  = Float64.(dataset[2:end, 1:end-1]);
targets = Float64.(dataset[2:end, end]);

# Load the DoME system
include("DoME.jl");
# Run DoME with these parameters
(trainingMSE, validationMSE, testMSE, bestTree) = dome(inputs, targets;
   minimumReductionMSE = 1e-6,
   maximumNodes = 30,
   strategy = StrategyExhaustive,
   showText = true
);
# Write the expression on screen
println("Best expression found: ", string(bestTree));
println("Best expression found (written in Latex): ", latexString(bestTree));
# If you want to rename the variable names, one of the easiest way is to do something like:
expression = string(bestTree);
expression = replace(expression, "X1" => "vendor");
expression = replace(expression, "X2" => "MYCT");
expression = replace(expression, "X3" => "MMIN");
expression = replace(expression, "X4" => "MMAX");
expression = replace(expression, "X5" => "CACH");
expression = replace(expression, "X6" => "CHMIN");
expression = replace(expression, "X7" => "CHMAX");
println("Best expression found (with the real names of the variables): ", expression);


text:
const STRING_GLOBALS = (
    mjDISABLESTRING = [
        "Constraint",
        "Equality",
        "Frictionloss",
        "Limit",
        "Contact",
        "Passive",
        "Gravity",
        "Clampctrl",
        "Warmstart",
        "Filterparent",
        "Actuation",
        "Refsafe",
    ],
    mjENABLESTRING = ["Override", "Energy", "Fwdinv", "Sensornoise"],
    mjTIMERSTRING = [
        "step",
        "forward",
        "inverse",
        "position",
        "velocity",
        "actuation",
        "acceleration",
        "constraint",
        "pos_kinematics",
        "pos_inertia",
        "pos_collision",
        "pos_make",
        "pos_project",
    ],
    mjLABELSTRING = [
        "None",
        "Body",
        "Joint",
        "Geom",
        "Site",
        "Camera",
        "Light",
        "Tendon",
        "Actuator",
        "Constraint",
        "Skin",
        "Selection",
        "SelPoint",
        "ContactForce",
    ],
    mjFRAMESTRING = ["None", "Body", "Geom", "Site", "Camera", "Light", "World"],
    mjVISSTRING = permutedims([
        "Convex Hull" "0" "H"
        "Texture" "1" "X"
        "Joint" "0" "J"
        "Actuator" "0" "U"
        "Camera" "0" "Q"
        "Light" "0" "Z"
        "Tendon" "1" "V"
        "Range Finder" "1" "Y"
        "Constraint" "0" "N"
        "Inertia" "0" "I"
        "Scale Inertia" "0" "'"
        "Perturb Force" "0" "B"
        "Perturb Object" "1" "O"
        "Contact Point" "0" "C"
        "Contact Force" "0" "F"
        "Contact Split" "0" "P"
        "Transparent" "0" "T"
        "Auto Connect" "0" "A"
        "Center of Mass" "0" "M"
        "Select Point" "0" "E"
        "Static Body" "1" "D"
        "Skin" "1" ";"
    ]),
    mjRNDSTRING = permutedims([
        "Shadow" "1" "S"
        "Wireframe" "0" "W"
        "Reflection" "1" "R"
        "Additive" "0" "L"
        "Skybox" "1" "K"
        "Fog" "0" "G"
        "Haze" "1" "/"
        "Segment" "0" ","
        "Id Color" "0" "."
    ]),
)

const CB_GLOBALS = (
    :mjcb_passive,
    :mjcb_control,
    :mjcb_contactfilter,
    :mjcb_sensor,
    :mjcb_time,
    :mjcb_act_dyn,
    :mjcb_act_gain,
    :mjcb_act_bias,
)

const mjCOLLISIONFUNC = [
    1  1  1  1  1  1  1  1
    1  1  1  1  1  1  1  1
    0  0  0  1  1  1  1  1
    0  0  0  0  1  1  1  1
    0  0  0  0  0  1  1  1
    0  0  0  0  0  0  1  1
    0  0  0  0  0  0  0  1
    0  0  0  0  0  0  0  0
]

@testset "$name" for (name, testval) in pairs(STRING_GLOBALS)
    @test getfield(MJCore, name) == testval
end

@testset "$name" for name in CB_GLOBALS
    ptr = Base.unsafe_convert(Ptr, getfield(MJCore, name))
    val = getfield(MJCore, name)[]
    @test ptr isa Ptr{Ptr{Cvoid}} && ptr != C_NULL
    @test val isa Ptr{Cvoid} && val == C_NULL
end

@testset "warn/error CB" begin
    for name in (:mju_user_error, :mju_user_warning)
        ptr = Base.unsafe_convert(Ptr, getfield(MJCore, name))
        val = getfield(MJCore, name)[]
        @test ptr isa Ptr{Ptr{Cvoid}} && ptr != C_NULL
        @test val isa Ptr{Cvoid} && val != C_NULL
    end
end

@test mjCOLLISIONFUNC == map(isequal(C_NULL), MJCore.mjCOLLISIONFUNC)


text:
struct DriftPath{T <: SSDFloat}
    path::Vector{<:AbstractCoordinatePoint{RealQuantity}}
    timestamps::Vector{RealQuantity}
end

struct EHDriftPath{T <: SSDFloat, TT <: RealQuantity}
    e_path::Vector{<:AbstractCoordinatePoint{T}}
    h_path::Vector{<:AbstractCoordinatePoint{T}}
    timestamps_e::Vector{TT}
    timestamps_h::Vector{TT}
end

_common_length(dp::EHDriftPath{T} where {T <: SSDFloat})::Int =
    max(length(dp.timestamps_e), length(dp.timestamps_h))
_common_length(dps::Vector{EHDriftPath{T}} where {T <: SSDFloat})::Int =
    maximum(_common_length.(dps))

function _common_time(dp::EHDriftPath{T, TT})::TT where {T <: SSDFloat, TT<:RealQuantity}
    max(last(dp.timestamps_e), last(dp.timestamps_h))
end
_common_time(dps::Vector{<:EHDriftPath}) =
maximum(_common_time.(dps))

function _common_timestamps(dp::Union{<:EHDriftPath{T}, Vector{<:EHDriftPath{T}}}, Δt) where {T}
    range(zero(Δt), step = Δt, stop = typeof(Δt)(_common_time(dp)) + Δt)
end

function get_velocity_vector(interpolation_field::Interpolations.Extrapolation{<:StaticVector{3}, 3}, point::CartesianPoint{T})::CartesianVector{T} where {T <: SSDFloat}
    return CartesianVector{T}(interpolation_field(point.x, point.y, point.z))
end

@inline function get_velocity_vector(interpolated_vectorfield, point::CylindricalPoint{T}) where {T <: SSDFloat}
    return CartesianVector{T}(interpolated_vectorfield(point.r, point.φ, point.z))
end

function _drift_charges(detector::SolidStateDetector{T}, grid::Grid{T, 3}, point_types::PointTypes{T, 3},
                        starting_points::Vector{CartesianPoint{T}},
                        velocity_field_e::Interpolations.Extrapolation{<:SVector{3}, 3},
                        velocity_field_h::Interpolations.Extrapolation{<:SVector{3}, 3},
                        Δt::RQ; max_nsteps::Int = 2000, verbose::Bool = true)::Vector{EHDriftPath{T}} where {T <: SSDFloat, RQ <: RealQuantity}

    drift_paths::Vector{EHDriftPath{T}} = Vector{EHDriftPath{T}}(undef, length(starting_points))

    dt::T = T(to_internal_units(internal_time_unit, Δt))

    for i in eachindex(starting_points)
        drift_path_e::Vector{CartesianPoint{T}} = zeros(CartesianPoint{T}, max_nsteps )#Vector{CartesianPoint{T}}(undef, max_nsteps)
        drift_path_h::Vector{CartesianPoint{T}} = zeros(CartesianPoint{T}, max_nsteps )#Vector{CartesianPoint{T}}(undef, max_nsteps)
        timestamps_e::Vector{T} = Vector{T}(undef, max_nsteps)
        timestamps_h::Vector{T} = Vector{T}(undef, max_nsteps)
        n_e::Int = _drift_charge!(drift_path_e, timestamps_e, detector, point_types, grid, starting_points[i], dt, velocity_field_e, verbose = verbose)
        n_h::Int = _drift_charge!(drift_path_h, timestamps_h, detector, point_types, grid, starting_points[i], dt, velocity_field_h, verbose = verbose)
        drift_paths[i] = EHDriftPath{T, T}( drift_path_e[1:n_e], drift_path_h[1:n_h], timestamps_e[1:n_e], timestamps_h[1:n_h] )
    end

    return drift_paths
end

function _drift_charge( detector::SolidStateDetector{T}, grid::Grid{T, 3}, point_types::PointTypes{T, 3},
                       starting_point::CartesianPoint{<:SSDFloat},
                       velocity_field_e::Interpolations.Extrapolation{SVector{3, T}, 3},
                       velocity_field_h::Interpolations.Extrapolation{SVector{3, T}, 3},
                       Δt::RealQuantity, max_nsteps::Int = 2000, verbose::Bool = true)::Vector{EHDriftPath{T}} where {T <: SSDFloat}
    return _drift_charges(detector, grid, CartesianPoint{T}.(point_types), [starting_point], velocity_field_e, velocity_field_h, T(Δt.val) * unit(Δt), max_nsteps = max_nsteps, verbose = verbose)
end

@inline _convert_vector(pt::CartesianPoint, ::Val{:cylindrical}) = CylindricalPoint(pt)
@inline _convert_vector(pt::CartesianPoint, ::Val{:cartesian}) = pt

function modulate_surface_drift(p::CartesianVector{T})::CartesianVector{T} where {T <: SSDFloat}
    return p
end

"""
    _drift_charge(...)

Before calling this function one should check that `startpos` is inside `det`: `in(startpos, det`
"""
function _drift_charge!(
                            drift_path::Vector{CartesianPoint{T}},
                            timestamps::Vector{T},
                            det::SolidStateDetector{T, S},
                            point_types::PointTypes{T, 3, S},
                            grid::Grid{T, 3, S},
                            startpos::CartesianPoint{T},
                            Δt::T,
                            velocity_field::Interpolations.Extrapolation{<:StaticVector{3}, 3};
                            verbose::Bool = true
                        )::Int where {T <: SSDFloat, S}
    drifttime::T = zero(T)
    done::Bool = false
    drift_path[1] = startpos
    timestamps[1] = zero(T)
    null_step::CartesianVector{T} = CartesianVector{T}(0, 0, 0)
    last_real_step_index::Int = 1
    @inbounds for istep in eachindex(drift_path)[2:end]
        if done == false
            last_real_step_index += 1
            current_pos::CartesianPoint{T} = drift_path[istep - 1]
            stepvector::CartesianVector{T} = get_velocity_vector(velocity_field, _convert_vector(current_pos, Val(S))) * Δt
            if geom_round.(stepvector) == null_step
                done = true
            end
            next_pos::CartesianPoint{T} = current_pos + stepvector
            next_pos_cyl::CylindricalPoint{T} = _convert_vector(next_pos, Val(S))
            if next_pos_cyl in point_types || next_pos_cyl in det
                drift_path[istep] = next_pos
                drifttime += Δt
                timestamps[istep] = drifttime
            else
                crossing_pos::CartesianPoint{T}, cd_point_type::UInt8, boundary_index::Int, surface_normal::CartesianVector{T} = get_crossing_pos(det, grid, current_pos, next_pos)
                if cd_point_type == CD_ELECTRODE
                    drift_path[istep] = crossing_pos
                    drifttime += Δt
                    timestamps[istep] = drifttime
                    done = true
                elseif cd_point_type == CD_FLOATING_BOUNDARY
                    projected_vector::CartesianVector{T} = CartesianVector{T}(project_to_plane(stepvector, surface_normal))
                    projected_vector = modulate_surface_drift(projected_vector)
                    next_pos = current_pos + projected_vector
                    # ToDo: We actually need a time array as well to do this properly...
                    small_projected_vector = projected_vector * T(0.001)
                    i::Int = 0
                    while i < 1000 && !(next_pos in det)
                        next_pos -= small_projected_vector
                        i += 1
                    end
                    if i == 1000 && verbose @warn("Handling of charge at floating boundary did not work as intended. Start Position (Cart): $startpos") end
                    drift_path[istep] = next_pos
                    drifttime += Δt * (1 - i * T(0.001))
                    timestamps[istep] = drifttime
                    if geom_round.(next_pos - current_pos) == null_step
                        done = true
                    end
                elseif cd_point_type == CD_BULK
                    if verbose @warn ("Internal error for charge starting at $startpos") end
                    drift_path[istep] = current_pos
                    drifttime += Δt
                    timestamps[istep] = drifttime
                    done = true
                else # elseif cd_point_type == CD_OUTSIDE
                    if verbose @warn ("Internal error for charge starting at $startpos") end
                    drift_path[istep] = current_pos
                    drifttime += Δt
                    timestamps[istep] = drifttime
                    done = true
                end
            end
        end
    end
    return last_real_step_index
end

function is_real(pt::AbstractCoordinatePoint{T})::Bool where {T <: SSDFloat}
    any(v -> isnan(v) || isinf(v), pt[:] )
end
function is_real(pt::AbstractCoordinateVector{T})::Bool where {T <: SSDFloat}
    any(v -> isnan(v) || isinf(v), pt[:] )
end

# Point types for charge drift: Defined in DetectorGeometries/DetectorGeometries.jl
# const CD_ELECTRODE = 0x00
# const CD_OUTSIDE = 0x01
# const CD_BULK = 0x02
# const CD_FLOATING_BOUNDARY = 0x04 # not 0x03, so that one could use bit operations here...

function get_crossing_pos(  detector::SolidStateDetector{T, S}, grid::Grid{T, 3}, point_in::CartesianPoint{T}, point_out::CartesianPoint{T};
                            max_n_iter::Int = 500)::Tuple{CartesianPoint{T}, UInt8, Int, CartesianVector{T}} where {T <: SSDFloat, S}
    point_mid::CartesianPoint{T} = T(0.5) * (point_in + point_out)
    cd_point_type::UInt8, contact_idx::Int, surface_normal::CartesianVector{T} = point_type(detector, grid, _convert_vector(point_mid, Val(S)))
    for i in 1:max_n_iter
        if cd_point_type == CD_BULK
            point_in = point_mid
        elseif cd_point_type == CD_OUTSIDE
            point_out = point_mid
        elseif cd_point_type == CD_ELECTRODE
            break
        else #elseif cd_point_type == CD_FLOATING_BOUNDARY
            break
        end
        point_mid = T(0.5) * (point_in + point_out)
        cd_point_type, contact_idx, surface_normal = point_type(detector, grid, _convert_vector(point_mid, Val(S)))
    end
    return point_mid, cd_point_type, contact_idx, surface_normal
end

include("plot_recipes.jl")