text:
# This file contains equilibria for specific quadratures

function _equilibrium(q::D2Q17, ρ, weight, u_dot_xi, u_squared, T, xi_squared)
    # Truncated upto order 2
    cs = q.speed_of_sound_squared
    D = dimension(q)
    # H_2_temperature = cs * (T .- 1) .* (cs * xi_squared - D)
    # H_3_temperature = 3.0 * cs * (T .- 1) * (cs * xi_squared - 2 - D)
    H_2_temperature = 0.0
    H_3_temperature = 0.0

    a_H_0 = 1.0
    a_H_1 = cs * u_dot_xi
    a_H_2 = cs^2 * (u_dot_xi .* u_dot_xi) .+ H_2_temperature .+ -cs * u_squared
    a_H_3 =
        cs * u_dot_xi .*
        (cs^2 * (u_dot_xi .* u_dot_xi) - 3 * cs * u_squared .+ H_3_temperature)
    return ρ .* weight .* (a_H_0 .+ a_H_1 .+ (1 / 2) * a_H_2 .+ (1 / 6) * a_H_3)
end

# function _equilibrium(q::D2Q17, ρ, weight, u_dot_xi, u_squared, T, xi_squared)
#     # Truncated upto order 2
#     cs = q.speed_of_sound_squared
#     D = dimension(q)
#     # H_2_temperature = (cs * T .- 1) .* (cs * xi_squared - D)
#     # H_3_temperature = 3.0 * (cs * T .- 1) * (cs * xi_squared - 2 - D)
#     H_2_temperature = 0.0
#     H_3_temperature = 0.0

#     a_H_0 = 1.0
#     a_H_1 = cs * u_dot_xi

#     a_H_2 = cs^2 * (u_dot_xi .* u_dot_xi) .+ H_2_temperature .+ - cs * u_squared
#     a_H_3 = cs * u_dot_xi .* (
#         cs^2 * (u_dot_xi .* u_dot_xi) -  3 * cs * u_squared .+ H_3_temperature
#     )
#     return ρ .* weight .* (
#         a_H_0 .+
#         a_H_1 .+
#         (1 / 2) * a_H_2 .+
#         (1 / 6) * a_H_3
#     )
# end

function _equilibrium(q::D2Q21, ρ, weight, u_dot_xi, u_squared, T, xi_squared)
    # Truncated upto order 2
    cs = q.speed_of_sound_squared
    D = dimension(q)
    # H_2_temperature = (cs * T .- 1) .* (cs * xi_squared - D)
    # H_3_temperature = 3.0 * (cs * T .- 1) * (cs * xi_squared - 2 - D)
    H_2_temperature = 0.0
    H_3_temperature = 0.0

    a_H_0 = 1.0
    a_H_1 = cs * u_dot_xi
    a_H_2 = cs^2 * (u_dot_xi .* u_dot_xi) .+ H_2_temperature .+ -cs * u_squared
    a_H_3 =
        cs * u_dot_xi .*
        (cs^2 * (u_dot_xi .* u_dot_xi) - 3 * cs * u_squared .+ H_3_temperature)
    return ρ .* weight .* (a_H_0 .+ a_H_1 .+ (1 / 2) * a_H_2 .+ (1 / 6) * a_H_3)
end

function equilibrium!(
    q::D2Q37,
    ρ::T,
    u::VT,
    temperature::T,
    f::VT,
) where {T <: Real, VT <: AbstractVector{T}}
    u_squared = zero(T)
    u_fourth = zero(T)
    for d in 1:dimension(q)
        u_squared += u[d] .^ 2
        u_fourth += u[d] .^ 4
    end

    for idx in 1:length(q.weights)
        u_dot_xi = q.abscissae[1, idx] .* u[1] .+ q.abscissae[2, idx] .* u[2]

        f[idx] = _equilibrium(
            q,
            ρ,
            q.weights[idx],
            u_dot_xi,
            u_squared,
            u_fourth,
            temperature,
            q.abscissae[1, idx]^2 + q.abscissae[2, idx]^2,
        )
    end

    return
end

function _equilibrium(
    q::D2Q37,
    ρ::T,
    weight::T,
    u_dot_xi::T,
    u_squared::T,
    u_fourth::T,
    temperature::T,
    xi_squared,
) where {T <: Real}
    cs = q.speed_of_sound_squared
    D = dimension(q)
    # H_2_temperature = (cs * temperature .- 1) .* (cs * xi_squared - D)
    # H_3_temperature = 3.0 * (cs * temperature .- 1) * (cs * xi_squared - 2 - D)
    H_2_temperature = zero(T)
    H_3_temperature = zero(T)

    a_H_0 = 1.0
    a_H_1 = cs * u_dot_xi
    a_H_2 = cs^2 * (u_dot_xi .* u_dot_xi) .+ H_2_temperature .+ -cs * u_squared
    a_H_3 =
        cs * u_dot_xi .*
        (cs^2 * (u_dot_xi .* u_dot_xi) - 3 * cs * u_squared .+ H_3_temperature)

    # H_4_temperature = 0.0
    a_H_4 = (cs^4 * u_dot_xi^4 - 6 * cs^3 * u_squared * u_dot_xi^2 + 3 * cs^2 * u_squared^2)

    return ρ .* weight .*
           (a_H_0 .+ a_H_1 .+ (1 / 2) * a_H_2 .+ (1 / 6) * a_H_3 .+ (1 / 24) * a_H_4)
end

# The D2Q4 lattice is unable to include temperature
pressure(q::D2Q4, f::VT, ρ::T, u::VT) where {T <: Real, VT <: AbstractVector{T}} = 1.0

function _equilibrium(
    q::D2Q4,
    ρ::T,
    weight::T,
    u_dot_xi::T,
    u_squared::T,
    temperature::T,
    xi_squared,
) where {T <: Real}
    cs = q.speed_of_sound_squared
    a_H_0 = one(T)
    a_H_1 = cs * u_dot_xi

    return ρ .* weight .* (a_H_0 .+ a_H_1)
end

function _equilibrium(
    q::D2Q5,
    ρ::T,
    weight::T,
    u_dot_xi::T,
    u_squared::T,
    temperature::T,
    xi_squared,
) where {T <: Real}
    cs = q.speed_of_sound_squared
    a_H_0 = one(T)
    a_H_1 = cs * u_dot_xi

    return ρ .* weight .* (a_H_0 .+ a_H_1)
end


text:
function remove_constraint!(data, idx::Int)
    constraint_idx = data.working_set[idx]
    ignored_constraints_add_row!(data, data.A[constraint_idx, :], data.b[constraint_idx])
    prepend!(data.ignored_set, constraint_idx)
    deleteat!(data.working_set, idx)
    remove_constraint!(data.F, idx)
    update_views!(data) # Updates A_ignored, b_ignored views
end

function add_constraint!(data, idx::Int)
    constraint_idx = data.ignored_set[idx]
    add_constraint!(data.F, data.A[constraint_idx, :])
    ignored_constraints_remove_row!(data, idx)
    deleteat!(data.ignored_set, idx)
    append!(data.working_set, constraint_idx)
    update_views!(data)
end

function ignored_constraints_add_row!(data, a::AbstractVector, b::AbstractFloat)
    l = length(data.ignored_set)
    data.A_shuffled[data.m - l, :] = a
    data.b_shuffled[data.m - l] = b
end

function ignored_constraints_remove_row!(data, idx::Int)
    l = length(data.ignored_set)
    start = data.m - l + 1
    @inbounds for j in 1:data.n, i in start+idx-1:-1:start+1
        data.A_shuffled[i, j] = data.A_shuffled[i - 1, j]
    end
    @inbounds for i in start+idx-1:-1:start+1
        data.b_shuffled[i] = data.b_shuffled[i - 1]
    end
end

function update_views!(data)
    range = data.m-length(data.ignored_set)+1:data.m
    data.A_ignored = view(data.A_shuffled, range, :)
    data.b_ignored = view(data.b_shuffled, range)
end

ext: jl


text:
using SetProg
using Documenter, Literate

const EXAMPLES_DIR = joinpath(@__DIR__, "src", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

const EXAMPLES = readdir(EXAMPLES_DIR)

for example in EXAMPLES
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR)
    Literate.notebook(example_filepath, OUTPUT_DIR)
    Literate.script(example_filepath, OUTPUT_DIR)
end

makedocs(
    sitename = "SetProg",
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    # See https://github.com/jump-dev/JuMP.jl/issues/1576
    strict = true,
    pages = [
        "Index" => "index.md",
        "Tutorials" => map(
            file -> joinpath("generated", file),
            filter(
                file -> endswith(file, ".md"),
                sort(readdir(OUTPUT_DIR)),
            ),
        ),
    ],
)

deploydocs(
    repo   = "github.com/blegat/SetProg.jl.git",
    push_preview = true,
)


text:
using Parameters

## Create a type which has default values:
@with_kw immutable PhysicalPara{R}
    rw::R = 1000.
    ri::R = 900.
    L::R = 3.34e5
    g::R = 9.81
    cw::R = 4220.
    day::R = 24*3600.
end

# Create an instance with the defaults
pp = PhysicalPara{Float64}()
# Make one with some non-defaults
pp2 = PhysicalPara{Float32}(cw=77, day= 987)
# Make another one based on the previous one with some modifications
pp3 = PhysicalPara(pp2; cw=.11e-7, rw=100.)

## It's possible to use @asserts straight in the type-def.  (Note, as
# usual, that for mutables, these asserts can be violated by updating
# the fields.)
@with_kw immutable PhysicalPara2{R}
    rw::R = 1000.; @assert rw>0
    ri::R = 900.
    @assert rw>ri # Note that the placement of assertions is not
                  # relevant. (They are moved to the constructor.
end

## Parameter interdependence
@with_kw immutable Para{R<:Real}
    a::R = 5
    b::R
    c::R = a+b
end
pa = Para{Int}(b=7)

## Setting a default type annotation, as often the bulk of fields will
# have the same type.  The last example more compactly:
@with_kw immutable Para2{R<:Real} @deftype R
    a = 5
    b
    c = a+b
end
pa2 = Para2{Int}(b=7)
# or more pedestrian
@with_kw immutable Para3 @deftype Float64
    a = 5
    b
    c = a+b
end
pa3 = Para3(b=7)

## Custom inner constructors:
@with_kw immutable MyS{R}
    a::R = 5
    b = 4

    # Can define inner constructors as long as:
    #  - one defining all positional arguments is given
    #  - no zero-positional arguments constructor is defined (as that
    #    would clash with the keyword constructor)
    #  - no @assert are used in the type-def (above example)
    #
    # Note that the keyword constructor goes through the positional
    # constructor, thus invariants defined there will be honored.

    MyS(a,b) = (@assert a>b; new(a,b))  # The keyword constructor
                                        # calls this constructor, so
                                        # the invariant is satisfied.
                                        # Note that invariants can be
                                        # done with @asserts as in
                                        # above example.
    MyS(a) = MyS{R}(a, a-1) # For this provide your own outer constructor:
end
MyS{R}(a::R) = MyS{R}(a)

MyS{Int}() # MyS(5,4)
ms = MyS(3) # MyS(3,2)
MyS(ms, b=-1) # MyS(3,-1)
try
    MyS(ms, b=6) # this will fail the assertion
end

## (Un)pack macros
#
# When working with parameters it is often convenient to unpack (and
# pack then):
function fn1(var, pa::Para)
    @unpack_Para pa # the macro is constructed during the @with_kw
                    # and called @unpack_*
    out = var + a + b
    b = 77
    @pack_Para pa # now pa.b==77
    return out, pa
end

out, pa = fn1(7, pa)

# If only a few parameters are needed, or possibly in general, it is
# more prudent to be explicit which variables are introduced into the
# local scope:

function fn2(var, pa::Para)
    @unpack a, b = pa
    out = var + a + b
    b = 77
    #@pack pa = b # doesn't work with immutable type
    return out, pa
end

out, pa = fn2(7, pa)


text:
Base.@kwdef struct Skeleton
    # Data
    y0::Vector{Float64}
    X0::Matrix{Float64}
    y1::Vector{Float64}
    X1::Matrix{Float64}
    mapping::Vector{Vector{Int}} = [[i] for i in 1:size(X0, 2)]
    update_g::Vector{Bool} = ones(Bool, length(mapping))
    # Transformed data
    N0::Int = size(X0, 1)
    N1::Int = size(X1, 1)
    D0::Int = size(X0, 2)
    D1::Int = size(X1, 2)
    X0vec::Vector{Vector{Float64}} = [X0[i, :] for i in 1:N0]
    X1vec::Vector{Vector{Float64}} = [X1[i, :] for i in 1:N1]
    # HyperParameters
    ζ0g::Float64 = 1.0
    m0β::Vector{Float64} = zeros(D0)
    Σ0β::Matrix{Float64} = Matrix(9.0 * I(D0))
    # Parameters
    rmodel::BNB.Sampler = BNB.Sampler(ones(Int, N0), X0; μ0β = m0β, Σ0β, ζ0g, mapping, update_g = zeros(Bool, length(mapping)), s = [2])
    r::Vector{Int} = ones(Int, N0)
    d::Vector{Int} = ones(Int, N0)
    β::Vector{Float64} = rmodel.β
    g::Vector{Bool} = rmodel.g
    s::Vector{Int} = rmodel.s
    # Transformed parameters
    rmax::Base.RefValue{Int} = Ref(maximum(r))
    ϕ0::Vector{Float64} = ones(N0) / 2
    ϕ1::Vector{Float64} = ones(N1) / 2
    f::Vector{Float64} = zeros(N1)
    n::Vector{Int} = [N0]
end

ygrid(skl::Skeleton) = skl.y1
Xgrid(skl::Skeleton) = skl.X1
Ngrid(skl::Skeleton) = skl.N1
y(skl::Skeleton) = skl.y0
X(skl::Skeleton) = skl.X0
N(skl::Skeleton) = skl.N0
cluster_sizes(skl::Skeleton) = skl.n
cluster_labels(skl::Skeleton) = skl.d
r(skl::Skeleton) = skl.r
f(skl::Skeleton) = skl.f
rmax(skl::Skeleton) = skl.rmax[]
# ok: 2022-03-22


text:
using VLFeat
using Images
using Distances

# Return [keypoints, descriptor]
# keypoints: Each column is a feature frame and has the format [X;Y;S;TH], where X,Y is the (fractional) center of the frame, S is the scale and TH is the orientation (in radians). [from VLFeat documentation]
# descriptors: 128xn array of x,y points
sift(im::Image) = VLFeat.vl_sift(im)

# Find matches (returns 2xn array of matches)
# Returns ([2xm array of matches], [euclidean distance measures])
# TODO: Replace with FLANN
# match A to B iff dist(A, B) * threshold < dist(A, [all other matches])
function find_matches(siftA, siftB, threshold::Float64)
    # Yes, this is very brute force right now
    nA = size(siftA, 2)
    nB = size(siftB, 2)
    bestIndex = zeros(Float32, nA)
    bestValue = fill(inf(Float32), nA)
    secondValue = fill(inf(Float32), nA)
    matchFlag = fill(false, nA)
    nMatches = 0
    for a = 1:nA
        for b = 1:nB
            dist = euclidean(siftA[:, a], siftB[:, b])
            if dist < bestValue[a]
                if bestValue[a] < secondValue[a]
                    secondValue[a] = bestValue[a]
                end
                bestValue[a] = dist
                bestIndex[a] = b
            end
        end
        if bestValue[a] < secondValue[a] * threshold
            matchFlag[a] = true
            nMatches += 1
        end
    end
    matches = zeros(Int32, 2, nMatches)
    dists = zeros(Float32, nMatches)
    cMatch = 1
    for a = 1:nA
        if matchFlag[a]
            matches[1, cMatch] = a
            matches[2, cMatch] = bestIndex[a]
            dists[cMatch] = bestValue[a]
            cMatch += 1
        end
    end
    return (matches, dists)
end

function indices_to_coordinates(coordsA, coordsB, matches)
    coords = zeros(Float32, 4, size(matches, 2))
    for match = 1:size(matches, 2)
        coords[1:2, match] = coordsA[:, matches[1, match]]
        coords[3:4, match] = coordsB[:, matches[2, match]]
    end
    return coords
end


text:
import Statistics:mean

function ParetoSub(X, perm)
    Z= X.-X[:,1]
    nondominated=reshape(any(Z.<0,dims=1),:)
    ispareto = all(any(Z[:,nondominated].>0, dims=1))
    X=X[:, nondominated]
    nondominated = perm[nondominated]
    return ispareto, nondominated, X
end

"""
    index = ParetoSet(X[, sense])
Get the Pareto set from a given set of oberservations `X`.
# Argument
- `X::AbstractMatrix{Real}` size(X)=(n_objectives, n_samples)
- `sense::AbstractVector{Union{Symbol, EGOSense}}=repeat(:Min, n_objectives)` sense of each objective, size(S)=(n_objectives,)
"""
function ParetoSet(X::AbstractMatrix)
    n,m = size(X)
    Xmin = minimum(X, dims=2)
    Xt = X.-Xmin
    Xm = mean(Xt, dims=2)
    Xn = maximum(Xt ./ (Xm.+maximum(Xm)), dims=1)
    perm = sortperm(reshape(Xn,:))
    Y = X[:, perm]
    membership = falses(m)
    while length(perm)>1
        k = perm[1]
        membership[k], perm, Y = ParetoSub(Y, perm)
    end
    membership[perm] .= true
    membership
end

ParetoSet(X, sense::AbstractVector{EGOSense}) = ParetoSet(X.*Int.(sense))

ParetoSet(X, sense::AbstractVector{Symbol}) = ParetoSet(X, sym2sense.(sense))


text:
# # A fit with BAT.jl
#
# In this example we show how to bring the PDF parametrisation and
# forward model together with `BAT.jl` to perform a fit of simulated data.
# This fit is a work in progress and just a starting point for verification
# of the method.

using BAT, DensityInterface
using PartonDensity
using QCDNUM
using Plots, Random, Distributions, ValueShapes, ParallelProcessingTools
using StatsBase, LinearAlgebra

gr(fmt=:png);

# ## Simulate some data

seed = 42
Random.seed!(seed) # for reproducibility

# We can start off by simulating some fake data for us to fit. This way,
# we know exactly what initial conditions we have specified and can check
# the validity of our inference, assuming the generative model is the one that is producing our data.
#
# This is a good first check to work with.
#
# ### Specify the input PDFs
#
# See the *Input PDF parametrisation and priors* example for more information on the
# definition of the input PDFs.

pdf_params = ValencePDFParams(λ_u=0.64, K_u=3.38, λ_d=0.67, K_d=4.73,
                              λ_g1=0.50, λ_g2=-0.63, K_g=4.23, λ_q=-0.23, weights=[5., 5., 1., 1., 1., 0.5, 0.5]);

plot_input_pdfs(pdf_params)

# ### Go from PDFs to counts in ZEUS detector bins
#
# Given the input PDFs, we can then evolve, calculate the cross sections, and fold through
# the ZEUS transfer matrix to get counts in bins. Here, we make use of some simple helper
# functions to do so. For more details, see the *Forward model* example.

# first specify QCDNUM inputs
qcdnum_grid = QCDNUMGrid(x_min=[1.0e-3], x_weights=[1], nx=100,
    qq_bounds=[1.0e2, 3.0e4], qq_weights=[1.0, 1.0], nq=50, spline_interp=3)
qcdnum_params = QCDNUMParameters(order=2, α_S=0.118, q0=100.0, grid=qcdnum_grid,
    n_fixed_flav=5, iqc=1, iqb=1, iqt=1, weight_type=1);

# now SPLINT and quark coefficients
splint_params = SPLINTParameters();
quark_coeffs = QuarkCoefficients();

# initialise QCDNUM
forward_model_init(qcdnum_grid, qcdnum_params, splint_params)

# run forward model
counts_pred_ep, counts_pred_em = forward_model(pdf_params, qcdnum_params,
    splint_params, quark_coeffs);

#
# take a poisson sample
nbins = size(counts_pred_ep)[1]
counts_obs_ep = zeros(UInt64, nbins)
counts_obs_em = zeros(UInt64, nbins)

for i in 1:nbins
    counts_obs_ep[i] = rand(Poisson(counts_pred_ep[i]))
    counts_obs_em[i] = rand(Poisson(counts_pred_em[i]))
end

#

plot(1:nbins, counts_pred_ep, label="Expected counts (eP)", color="blue")
plot!(1:nbins, counts_pred_em, label="Expected counts (eM)", color="red")
scatter!(1:nbins, counts_obs_ep, label="Detected counts (eP)", color="blue")
scatter!(1:nbins, counts_obs_em, label="Detected counts (eM)", color="red")
plot!(xlabel="Bin number")

# store
sim_data = Dict{String, Any}()
sim_data["nbins"] = nbins;
sim_data["counts_obs_ep"] = counts_obs_ep;
sim_data["counts_obs_em"] = counts_obs_em;

# write to file
pd_write_sim("output/simulation.h5", pdf_params, sim_data)

# ## Fit the simulated data
#
# Now we can try to fit this simulated data using `Bat.jl`.
# The first step is to define the prior and likelihood.
# For now, let's try relatively narrow priors centred on the true values.

prior = NamedTupleDist(
    θ_tmp = Dirichlet(pdf_params.weights),
    λ_u = Truncated(Normal(pdf_params.λ_u, 1), 0, 1),
    K_u = Truncated(Normal(pdf_params.K_u, 1), 2, 10),
    λ_d = Truncated(Normal(pdf_params.λ_d, 1), 0, 1),
    K_d = Truncated(Normal(pdf_params.K_d, 1), 2, 10),
    λ_g1 = Truncated(Normal(pdf_params.λ_g1, 1), -1, 0),
    λ_g2 = Truncated(Normal(pdf_params.λ_g2, 1), -1, 0),
    K_g =  Truncated(Normal(pdf_params.K_g, 1), 2, 10),
    λ_q = Truncated(Normal(pdf_params.λ_q, 0.1), -1, 0),
);

# The likelihood is similar to that used in the *input PDF parametrisation* example.
# We start by accessing the current parameter set of the sampler's iteration,
# then running the forward model to get the predicted counts and comparing to
# the observed counts using a simple Poisson likelihood.
#
# The `@critical` macro is used because `forward_model()` is currently not thread safe, so
# this protects it from being run in parallel.

likelihood = let d = sim_data

    counts_obs_ep = d["counts_obs_ep"]
    counts_obs_em = d["counts_obs_em"]
    nbins = d["nbins"]

    logfuncdensity(function (params)

            θ = get_scaled_θ(params.λ_u, params.K_u, params.λ_d,
                             params.K_d, Vector(params.θ_tmp))

            pdf_params = ValencePDFParams(λ_u=params.λ_u, K_u=params.K_u, λ_d=params.λ_d,
                                          K_d=params.K_d, λ_g1=params.λ_g1, λ_g2=params.λ_g2,
                                          K_g=params.K_g, λ_q=params.λ_q, θ=θ)

            counts_pred_ep, counts_pred_em = @critical forward_model(pdf_params,
                qcdnum_params, splint_params, quark_coeffs);

            ll_value = 0.0
            for i in 1:nbins

                if counts_pred_ep[i] < 0
                   #@warn "counts_pred_ep[i] < 0, setting to 0" i counts_pred_ep[i]
                   counts_pred_ep[i] = 0
                end

                if counts_pred_em[i] < 0
                   #@warn "counts_pred_em[i] < 0, setting to 0" i counts_pred_em[i]
                   counts_pred_em[i] = 0
                end

                ll_value += logpdf(Poisson(counts_pred_ep[i]), counts_obs_ep[i])
                ll_value += logpdf(Poisson(counts_pred_em[i]), counts_obs_em[i])
            end

            return ll_value
    end)
end

# We can now run the MCMC sampler. We will start by using the
# Metropolis-Hastings algorithm as implemented in `BAT.jl`.
# To get reasonable results, we need to run the sampler for a
# long time (several hours). To save time in this demo, we will
# work with a ready-made results file. To actually run the sampler,
# simply uncomment the code below.

#posterior = PosteriorDensity(likelihood, prior);
#convergence = BrooksGelmanConvergence(threshold=1.3);
#burnin = MCMCMultiCycleBurnin(max_ncycles=100);

#samples = bat_sample(posterior, MCMCSampling(mcalg=MetropolisHastings(), nsteps=10^5, nchains=2)).result;

# Alternatively, we could also try a nested sampling approach
# here for comparison. This is easily done thanks to the
# interface of `BAT.jl`, you will just need to add the
# `NestedSamplers.jl` package.

#import NestedSamplers
#samples = bat_sample(posterior, EllipsoidalNestedSampling()).result

# If you run the sampler, be sure to save
# the results for further analysis

#import HDF5
#bat_write("output/results.h5", samples)

# ## Analysing the results
#
# First, let's load our simulation inputs and results

pdf_params, sim_data = pd_read_sim("output/demo_simulation.h5");
samples = bat_read("output/demo_results.h5").result;

# We can check some diagnostics using built in `BAT.jl`, such as the
# effective sample size shown below

bat_eff_sample_size(unshaped.(samples))[1]

# We see a value for each of our 15 total parameters. As the
# Metropolis-Hastings algorithm's default implementation
# isn't very efficient, we see that the effective sample size
# is only a small percentage of the input `nsteps`. We should try
# to improve this if possible, or use a much larger `nsteps` value.
#
# For demonstration purposes, we will continue to show how we can
# visualise the results in this case. For robust inference, we need
# to improve the sampling stage above.

# We can use `BAT.jl`'s built in plotting recipes to show the marginals,
# for example, consider `λ_u`, and compare to the known truth.

plot(
    samples, :(λ_u),
    nbins=50,
    colors=[:skyblue4, :skyblue3, :skyblue1],
    alpha=0.7,
    marginalmode=false,
    legend=:topleft
)
vline!([pdf_params.λ_u], color="black", label="truth", lw=3)

# If we want to compare the momentum weights, we must transform
# from `θ_tmp` to `θ`, as shown below. Here, we transform using
# a helper function, convert the result to a matrix, and
# access the ith weight with the integer `i`.

θ = [get_scaled_θ(λ_u, K_u, λ_d, K_d, Vector(θ_tmp)) for (λ_u, K_u, λ_d, K_d, θ_tmp)
     in zip(samples.v.λ_u, samples.v.K_u, samples.v.λ_d, samples.v.K_d, samples.v.θ_tmp)];

θ = transpose(reduce(vcat,transpose.(θ)))

i = 1
hist = append!(Histogram(0:0.02:1), θ[i,:])
plot(
    normalize(hist, mode=:density),
    st = :steps, label="Marginal posterior"
)
vline!([pdf_params.θ[i]], color="black", label="truth", lw=3)

# Rather than making a large plot 15 different marginals,
# it can be more useful to visualise the posterior distribution
# in differently, such as the shape of the distributions
# we are trying to fit, or the *model space*. Helper functions
# exist for doing just this.

# Using BAT recipe
function wrap_xtotx(p::NamedTuple{(:K_d, :K_g, :K_u, :θ_tmp, :λ_d, :λ_g1,
                                   :λ_g2, :λ_q, :λ_u)}, x::Real)
    θ = get_scaled_θ(p.λ_u, p.K_u, p.λ_d, p.K_d, Vector(p.θ_tmp))
    pdf_params = ValencePDFParams(λ_u=p.λ_u, K_u=p.K_u, λ_d=p.λ_d, K_d=p.K_d, λ_g1=p.λ_g1,
        λ_g2=p.λ_g2, K_g=p.K_g, λ_q=p.λ_q, θ=θ)
    return log(xtotx(x, pdf_params))
end

x_grid = range(1e-3, stop=1, length=50)
plot(x_grid, wrap_xtotx, samples, colors=[:skyblue4, :skyblue3, :skyblue1],
     legend=:topright)
plot!(x_grid, [log(xtotx(x, pdf_params)) for x in x_grid], color="black", lw=3,
      label="Truth", linestyle=:dash)
plot!(ylabel="log(xtotx)")

# Using `PartonDensity.jl`
plot_model_space(pdf_params, samples, nsamples=500)

# Alternatively, we can also visualise the implications of the fit
# in the *data space*, as shown below.

plot_data_space(pdf_params, sim_data, samples, qcdnum_grid, qcdnum_params,
                splint_params, quark_coeffs, nsamples=500)

# The first results seem promising, but these are really just first checks
# and more work will have to be done to verify the method.


text:
using SALTBase
using AbbreviatedTypes
using Test

# Define a concrete subtype of AbstractMaxwellSolver for the test purspose.
struct DefaultMXS <: AbstractMaxwellSolver end
Base.similar(::DefaultMXS) = DefaultMXS()
Base.size(::DefaultMXS) = (0,0)

function SALTBase.init_mxsolver!(mxsolver::DefaultMXS, ω::Number, ε::AbsVecNumber) end
function SALTBase.mxsolve!(x::AbsVecComplexF, mxsolver::DefaultMXS, b::AbsVecNumber) end
function SALTBase.mxsolve_transpose!(x::AbsVecComplexF, mxsolver::DefaultMXS, b::AbsVecNumber) end
function SALTBase.mxapply!(b::AbsVecComplexF, mxsolver::DefaultMXS, x::AbsVecNumber) end
function SALTBase.abs2_interp!(abs2ψ::AbsVecReal, mxsolver::DefaultMXS, ψ::AbsVecNumber) end

include("util.jl")
include("gain.jl")
include("lasing.jl")


text:
function solve(case, N, D, A)
    dict = Dict{Rational,Int}()
    for a in A
        dict[a] = get(dict, a, 0) + 1
    end

    cuts = 0
    while max(values(dict)) < D

    end

    println("Case #$(case): $(cuts)")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        N, D = parse.(Int, stdin |> readline |> split)
        A = parse.(Int, stdin |> readline |> split)
        solve(case, N, D, A)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# Wrong Answer

const offsets = (N = (0, -1), S = (0, 1), E = (1, 0), W = (-1, 0))

function solveRec(X, i)
    multiplier = X[i]
    offset = [0, 0]
    i += 2 # skip number and open paren

    while true
        if X[i] == ')'
            offset .*= parse(Int, multiplier)
            return offset, i
        end

        if isdigit(X[i])
            suboffset, close = solveRec(X, i)
            i = close + 1
            offset .+= suboffset
        else
            offset .+= offsets[Symbol(X[i])]
            i += 1
        end
    end
end

function solve(case, X)
    coords = [1, 1]

    i = firstindex(X)
    while i <= lastindex(X)
        if isdigit(X[i])
            offset, close = solveRec(X, i)
            coords .+= offset
            i = close + 1
        else
            coords .+= offsets[Symbol(X[i])]
            i += 1
        end
    end

    x, y = mod1.(coords, Int(1e9))
    println("Case #$(case): $x $y")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        path = readline(stdin)
        solve(case, path)
    end
end

isinteractive() || main()


text:
using Test
using DataDeps: try_determine_load_path, determine_save_path, try_determine_package_datadeps_dir

@testset "package data deps dir" begin
    target = joinpath(realpath(joinpath(dirname(@__FILE__),"..")),"deps","data")
    @test try_determine_package_datadeps_dir(@__FILE__) == target

    mktemp() do fn, fh
        @test try_determine_package_datadeps_dir(fn) == nothing
    end

end


text:
importall Bukdu

type WelcomeController <: ApplicationController
end

function plugins(c::WelcomeController)
    plug(Logger.log_message, c)
end

function index(::WelcomeController)
    Logger.info("hi")
    "hello"
end

Router() do
    get("/", WelcomeController, index)
end

Endpoint() do
    plug(Plug.Logger, level=:debug)
    plug(Router)
end

using Base.Test

Logger.have_color(false)
let oldout = STDERR
   rdout, wrout = redirect_stdout()

conn = (Router)(get, "/")

   reader = @async readstring(rdout)
   redirect_stdout(oldout)
   close(wrout)

   lines = split(wait(reader), '\n')
@test startswith(lines[1], "DEBUG GET /")
@test lines[2] == "INFO WelcomeController.index hi"
end


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
sz = MPI.Comm_size(comm)
filename = MPI.bcast(tempname(), 0, comm)

MPI.Barrier(comm)

# Collective write
fh = MPI.File.open(comm, filename, read=true, write=true, create=true)
@test MPI.File.get_position_shared(fh) == 0

if !MPI.File.get_atomicity(fh)
    MPI.File.set_atomicity(fh, true)
end

@test MPI.File.get_atomicity(fh)

MPI.Barrier(comm)
MPI.File.sync(fh)

header = "my header"

if rank == 0
    MPI.File.write_shared(fh, header)
end

# TODO: is there a better way to synchronise shared pointers?
MPI.Barrier(comm)
MPI.File.sync(fh)

offset = MPI.File.get_position_shared(fh)
@test offset == sizeof(header)
byte_offset = MPI.File.get_byte_offset(fh, offset)
@test byte_offset == offset

MPI.File.set_view!(fh, byte_offset, MPI.Datatype(Int64), MPI.Datatype(Int64))
MPI.Barrier(comm)
MPI.File.sync(fh)
@test MPI.File.get_position_shared(fh) == 0

MPI.Barrier(comm)
MPI.File.sync(fh)

MPI.File.write_ordered(fh, fill(Int64(rank), rank+1))
@test MPI.File.get_position_shared(fh) == sum(1:sz)

MPI.File.seek_shared(fh, 0)
@test MPI.File.get_position_shared(fh) == 0

MPI.Barrier(comm)
MPI.File.sync(fh)

buf = zeros(Int64, rank+1)
MPI.File.read_ordered!(fh, buf)
@test buf == fill(Int64(rank), rank+1)

MPI.Barrier(comm)
MPI.File.sync(fh)
@test MPI.File.get_position_shared(fh) == sum(1:sz)

MPI.File.set_view!(fh, 0, MPI.Datatype(UInt8), MPI.Datatype(UInt8))
MPI.Barrier(comm)
MPI.File.sync(fh)
MPI.File.seek_shared(fh, 0)
@test MPI.File.get_position_shared(fh) == 0

MPI.Barrier(comm)
MPI.File.sync(fh)

if rank == sz-1
    buf = Array{UInt8}(undef, sizeof(header))
    MPI.File.read_shared!(fh, buf)
    @test String(buf) == header
end

MPI.Barrier(comm)
MPI.File.sync(fh)

@test MPI.File.get_position_shared(fh) == sizeof(header)

close(fh)


text:
# These tests require an FT device which supports D2XX to be connected
#
# By Reuben Hill 2019, Gowerlabs Ltd, reuben@gowerlabs.co.uk
#
# Copyright (c) Gowerlabs Ltd.

module TestLibFTD2XX

using Test
using LibFTD2XX
import LibFTD2XX.Wrapper

@testset "high level" begin

  # libversion
  if Sys.iswindows()
    ver = libversion()
    @test ver isa VersionNumber
  else
    @test_throws UndefVarError libversion()
  end

  # createdeviceinfolist
  numdevs = LibFTD2XX.createdeviceinfolist()
  @test numdevs > 0

  # LibFTD2XX.getdeviceinfodetail
  @test_throws D2XXException LibFTD2XX.getdeviceinfodetail(numdevs)
  for deviceidx = 0:(numdevs-1)
    idx, flgs, typ, devid, locid, serialn, descr, fthand = LibFTD2XX.getdeviceinfodetail(deviceidx)
    @test idx == deviceidx
    if Sys.iswindows() # should not have a locid on windows
      @test locid == 0
    end
    @test serialn isa String
    @test descr isa String
    @test fthand isa FT_HANDLE
  end
  idx, flgs, typ, devid, locid, serialn, descr, fthand = LibFTD2XX.getdeviceinfodetail(0)
  @test_throws DomainError LibFTD2XX.getdeviceinfodetail(-1)

  # FT_HANDLE functions...
  @testset "FT_HANDLE" begin

    # open by description
    if !isempty(descr)
      handle = open(descr, OPEN_BY_DESCRIPTION)
      @test handle isa FT_HANDLE
      @test isopen(handle)
      @test_throws Wrapper.FT_DEVICE_NOT_FOUND open(descr, OPEN_BY_DESCRIPTION) # can't open twice
      close(handle)
      @test !isopen(handle)
  end

    # open by serialnumber
    if !isempty(serialn)
      handle = open(serialn, OPEN_BY_SERIAL_NUMBER)
      @test handle isa FT_HANDLE
      @test isopen(handle)
      @test_throws Wrapper.FT_DEVICE_NOT_FOUND open(serialn, OPEN_BY_SERIAL_NUMBER) # can't open twice
      close(handle)
      @test !isopen(handle)
  end

    # bytesavailable
    handle = open(descr, OPEN_BY_DESCRIPTION)
    nb = bytesavailable(handle)
    @test nb >= 0
    close(handle) # can't use on closed device
    @test_throws D2XXException bytesavailable(handle)

    # read
    handle = open(descr, OPEN_BY_DESCRIPTION)
    rxbuf = read(handle, nb)
    @test length(rxbuf) == nb
    @test_throws ArgumentError read(handle, -1) # exception type set by Base/io.jl
    close(handle) # can't use on closed device
    @test_throws D2XXException read(handle, nb)

    # write
    handle = open(descr, OPEN_BY_DESCRIPTION)
    txbuf = ones(UInt8, 10)
    nwr = write(handle, txbuf)
    @test nwr == length(txbuf)
    @test txbuf == ones(UInt8, 10)
    @test_throws ErrorException write(handle, Int.(txbuf)) # No byte I/O...
    close(handle) # can't use on closed device
    @test_throws D2XXException read(handle, nb)

    # readavailable
    handle = open(descr, OPEN_BY_DESCRIPTION)
    rxbuf = readavailable(handle)
    @test rxbuf isa AbstractVector{UInt8}
    close(handle) # can't use on closed device
    @test_throws D2XXException readavailable(handle)

    # baudrate
    handle = open(descr, OPEN_BY_DESCRIPTION)
    retval = baudrate(handle, 2000000)
    @test retval == nothing
    txbuf = ones(UInt8, 10)
    nwr = write(handle, txbuf)
    @test nwr == length(txbuf)
    @test txbuf == ones(UInt8, 10)
    @test_throws DomainError baudrate(handle, 0)
    @test_throws DomainError baudrate(handle, -1)
    close(handle) # can't use on closed device
    @test_throws D2XXException baudrate(handle, 2000000)

    # flush and eof
    handle = open(descr, OPEN_BY_DESCRIPTION)
    retval = flush(handle)
    @test eof(handle)
    @test retval == nothing
    @test isopen(handle)
    close(handle) # can't use on closed device
    @test_throws D2XXException flush(handle)
    @test_throws D2XXException eof(handle)

    # driverversion
    if Sys.iswindows()
      handle = open(descr, OPEN_BY_DESCRIPTION)
      ver = driverversion(handle)
      @test ver isa VersionNumber
      close(handle) # can't use on closed device
      @test_throws D2XXException driverversion(handle)
    else
      handle = open(descr, OPEN_BY_DESCRIPTION)
      @test_throws UndefVarError driverversion(handle)
      close(handle) # can't use on closed device
    end

    # datacharacteristics
    handle = open(descr, OPEN_BY_DESCRIPTION)
    retval = datacharacteristics(handle, wordlength = BITS_8, stopbits = STOP_BITS_1, parity = PARITY_NONE)
    @test retval == nothing
    close(handle) # can't use on closed device
    @test_throws D2XXException datacharacteristics(handle, wordlength = BITS_8, stopbits = STOP_BITS_1, parity = PARITY_NONE)

    # timeouts tests...
    handle = open(descr, OPEN_BY_DESCRIPTION)
    baudrate(handle, 9600)
    timeout_read, timeout_wr = 200, 100 # milliseconds
    timeouts(handle, timeout_read, timeout_wr)
    tread = 1000 * @elapsed read(handle, 5000)
    buffer = zeros(UInt8, 5000);
    twr = 1000 * @elapsed write(handle, buffer)
    @test timeout_read < 1.5*tread
    @test timeout_wr < 1.5*twr
    @test_throws DomainError timeouts(handle, timeout_read, -1)
    @test_throws DomainError timeouts(handle, -1, timeout_wr)
    close(handle) # can't use on closed device
    @test_throws D2XXException timeouts(handle, timeout_read, timeout_wr)

    # status
    handle = open(descr, OPEN_BY_DESCRIPTION)
    mflaglist, lflaglist = status(handle)
    @test mflaglist isa Dict{String, Bool}
    @test lflaglist isa Dict{String, Bool}
    @test haskey(mflaglist, "CTS")
    @test haskey(mflaglist, "DSR")
    @test haskey(mflaglist, "RI")
    @test haskey(mflaglist, "DCD")
    @test haskey(lflaglist, "OE")
    @test haskey(lflaglist, "PE")
    @test haskey(lflaglist, "FE")
    @test haskey(lflaglist, "BI")
    close(handle) # can't use on closed device
    @test_throws D2XXException status(handle)

    # close and isopen
    handle = open(descr, OPEN_BY_DESCRIPTION)
    retval = close(handle)
    @test retval == nothing
    @test !isopen(handle)
    @test LibFTD2XX.Wrapper._ptr(handle) == C_NULL
    retval = close(handle) # check can close more than once without issue...
    @test !isopen(handle)
  end

  # D2XXDevice
  @testset "D2XXDevice" begin

    # Constructor
    @test_throws DomainError D2XXDevice(-1)
    for i = 0:(numdevs-1)
      idx, flgs, typ, devid, locid, serialn, descr, fthand = LibFTD2XX.getdeviceinfodetail(i)
      dev = D2XXDevice(i)
      @test deviceidx(dev) == idx == i
      @test deviceflags(dev) == flgs
      @test devicetype(dev) == typ
      @test deviceid(dev) == devid
      if Sys.iswindows()
        @test locationid(dev) == locid == 0
      else
        @test locationid(dev) == locid
      end
      @test serialnumber(dev) == serialn
      @test description(dev) == descr
      @test LibFTD2XX.Wrapper._ptr(fthandle(dev)) == LibFTD2XX.Wrapper._ptr(fthand)
      @test !isopen(fthandle(dev))
    end

    # D2XXDevices
    devices = D2XXDevices()
    @test length(devices) == numdevs
    @test all(deviceidx(devices[d]) == deviceidx(D2XXDevice(d-1)) for d = 1:numdevs)

    # isopen
    @test all(.!isopen.(devices))

    # open
    retval = open.(devices)
    @test all(retval .== nothing)
    @test all(isopen.(devices))
    @test_throws D2XXException open.(devices) # can't open twice

    # bytesavailable
    nbs = bytesavailable.(devices)
    @test all(nbs .>= 0)
    close.(devices) # can't use on closed device
    @test_throws D2XXException bytesavailable.(devices)

    device = devices[1] # choose device 1...
    nb = nbs[1]

    # read
    open(device)
    rxbuf = read(device, nb)
    @test length(rxbuf) == nb
    @test_throws ArgumentError read(device, -1) # exception type set by Base/io.jl
    close(device) # can't use on closed device
    @test_throws D2XXException read(device, nb)

    # write
    open(device)
    txbuf = ones(UInt8, 10)
    nwr = write(device, txbuf)
    @test nwr == length(txbuf)
    @test txbuf == ones(UInt8, 10)
    @test_throws ErrorException write(device, Int.(txbuf)) # No byte I/O...
    close(device) # can't use on closed device
    @test_throws D2XXException write(device, txbuf)

    # readavailable
    open(device)
    rxbuf = readavailable(device)
    @test rxbuf isa AbstractVector{UInt8}
    close(device) # can't use on closed device
    @test_throws D2XXException readavailable(device)

    # baudrate
    open(device)
    retval = baudrate(device, 2000000)
    @test retval == nothing
    txbuf = ones(UInt8, 10)
    nwr = write(device, txbuf)
    @test nwr == length(txbuf)
    @test txbuf == ones(UInt8, 10)
    @test_throws DomainError baudrate(device, 0)
    @test_throws DomainError baudrate(device, -1)
    close(device) # can't use on closed device
    @test_throws D2XXException baudrate(device, 2000000)

    # flush and eof
    open(device)
    retval = flush(device)
    @test eof(device)
    @test retval == nothing
    @test isopen(device)
    close(device) # can't use on closed device
    @test_throws D2XXException flush(device)

    # driverversion
    if Sys.iswindows()
      open(device)
      ver = driverversion(device)
      @test ver isa VersionNumber
      close(device) # can't use on closed device
      @test_throws D2XXException driverversion(device)
    else
      @test_throws UndefVarError driverversion(device)
    end

    # datacharacteristics
    open(device)
    retval = datacharacteristics(device, wordlength = BITS_8, stopbits = STOP_BITS_1, parity = PARITY_NONE)
    @test retval == nothing
    close(device) # can't use on closed device
    @test_throws D2XXException datacharacteristics(device, wordlength = BITS_8, stopbits = STOP_BITS_1, parity = PARITY_NONE)

    # timeouts tests...
    open(device)
    baudrate(device, 9600)
    timeout_read, timeout_wr = 200, 100 # milliseconds
    timeouts(device, timeout_read, timeout_wr)
    tread = 1000 * @elapsed read(device, 5000)
    buffer = zeros(UInt8, 5000);
    twr = 1000 * @elapsed write(device, buffer)
    @test timeout_read < 1.5*tread
    @test timeout_wr < 1.5*twr
    @test_throws DomainError timeouts(device, timeout_read, -1)
    @test_throws DomainError timeouts(device, -1, timeout_wr)
    close(device) # can't use on closed device
    @test_throws D2XXException timeouts(device, timeout_read, timeout_wr)

    # status
    open(device)
    mflaglist, lflaglist = status(device)
    @test mflaglist isa Dict{String, Bool}
    @test lflaglist isa Dict{String, Bool}
    @test haskey(mflaglist, "CTS")
    @test haskey(mflaglist, "DSR")
    @test haskey(mflaglist, "RI")
    @test haskey(mflaglist, "DCD")
    @test haskey(lflaglist, "OE")
    @test haskey(lflaglist, "PE")
    @test haskey(lflaglist, "FE")
    @test haskey(lflaglist, "BI")
    close(device) # can't use on closed device
    @test_throws D2XXException status(device)

    # close and isopen (all devices)
    retval = close.(devices)
    @test all(retval .== nothing)
    @test all(.!isopen.(devices))
    @test all(LibFTD2XX.Wrapper._ptr.(fthandle.(devices)) .== C_NULL)
    close.(devices) # check can close more than once without issue...
    @test all(.!isopen.(devices))

  end

end

end # module TestLibFTD2XX


text:
"""
    Output

Abstract supertype for simulation outputs.

Outputs are store or display simulation results, usually
as a vector of grids, one for each timestep - but they may also
sum, combine or otherwise manipulate the simulation grids to improve
performance, reduce memory overheads or similar.

Simulation outputs are decoupled from simulation behaviour,
and in many cases can be used interchangeably.
"""
abstract type Output{T,A} <: AbstractDimArray{T,1,Tuple{Ti},A} end
# Generic ImageOutput constructor. Converts an init array to vector of arrays.
function (::Type{T})(
    init::Union{NamedTuple,AbstractMatrix}; extent=nothing, kw...
) where T <: Output
    extent = extent isa Nothing ? Extent(; init=init, kw...) : extent
    T(; frames=[deepcopy(init)], running=false, extent=extent, kw...)
end

# Forward base methods to the frames array
Base.parent(o::Output) = frames(o)
Base.length(o::Output) = length(parent(o))
Base.size(o::Output) = size(parent(o))
Base.firstindex(o::Output) = firstindex(parent(o))
Base.lastindex(o::Output) = lastindex(parent(o))
Base.push!(o::Output, x) = push!(parent(o), x)
Base.step(o::Output) = step(tspan(o))

# DimensionalData interface
function DimensionalData.dims(o::Output)
    ts = tspan(o)
    val = isstored(o) ? ts : ts[end]:step(ts):ts[end]
    (Ti(val; mode=Sampled(Ordered(), Regular(step(ts)), Intervals(Start()))),)
end
DimensionalData.refdims(o::Output) = ()
DimensionalData.name(o::Output) = NoName()
DimensionalData.metadata(o::Output) = NoMetadata()
# Output bebuild just returns a DimArray
DimensionalData.rebuild(o::Output, data, dims::Tuple, refdims, name, metadata) =
    DimArray(data, dims, refdims, name, metadata)

# Getters and setters
frames(o::Output) = o.frames
isrunning(o::Output) = o.running
extent(o::Output) = o.extent
init(o::Output) = init(extent(o))
mask(o::Output) = mask(extent(o))
aux(o::Output, key...) = aux(extent(o), key...)
tspan(o::Output) = tspan(extent(o))
timestep(o::Output) = step(tspan(o))

ruleset(o::Output) =
    throw(ArgumentError("No ruleset on the output. Pass one to `sim!` as the second argument"))
fps(o::Output) = nothing
stoppedframe(o::Output) = lastindex(o)

setrunning!(o::Output, val) = o.running = val
settspan!(o::Output, tspan) = settspan!(extent(o), tspan)
setfps!(o::Output, x) = nothing
settimestamp!(o::Output, f) = nothing
setstoppedframe!(o::Output, f) = nothing

isasync(o::Output) = false
isstored(o::Output) = true
isshowable(o::Output, frame) = false
initialise!(o::Output, data) = nothing
finalise!(o::Output, data) = nothing
initialisegraphics(o::Output, data) = nothing
finalisegraphics(o::Output, data) = nothing
maybesleep(o::Output, frame) = nothing
showframe(o::Output, data) = nothing

frameindex(o::Output, data::AbstractSimData) = frameindex(o, currentframe(data))
frameindex(o::Output, f::Int) = isstored(o) ? f : oneunit(f)

function storeframe!(output::Output, data)
    checkbounds(output, frameindex(output, data))
    _storeframe!(eltype(output), output, data)
end

function _storeframe!(::Type{<:NamedTuple}, output::Output, data)
    map(values(grids(data)), keys(data)) do grid, key
        _copyto_output!(output[frameindex(output, data)][key], grid, proc(data))
    end
end
function _storeframe!(::Type{<:AbstractArray}, output::Output, data)
    grid = first(grids(data))
    _copyto_output!(output[frameindex(output, data)], grid, proc(data))
end

function _copyto_output!(outgrid, grid::GridData, proc)
    copyto!(outgrid, CartesianIndices(outgrid), source(grid), CartesianIndices(outgrid))
end
function _copyto_output!(outgrid, grid::GridData, proc::ThreadedCPU)
    Threads.@threads for j in axes(outgrid, 2)
        for i in axes(outgrid, 1)
            @inbounds outgrid[i, j] = grid[i, j]
        end
    end
end

# Grids are preallocated and reused.
init_output_grids!(o::Output, init) = init_output_grids!(o[1], o::Output, init)
# Array grids are copied
function init_output_grids!(grid::AbstractArray, o::Output, init::AbstractArray)
    copyto!(grid, init)
    return o
end
# All arrays are copied if both are named tuples
function init_output_grids!(grids::NamedTuple, o::Output, inits::NamedTuple)
    map(grids, inits) do grid, init
        copyto!(grid, init)
    end
    return o
end


text:
@enum TerminalConfig::Cuint RSE=Val_RSE NRSE=Val_NRSE Differential=Val_Diff PseudoDifferential=Val_PseudoDiff
@enum InputType::Cuint Voltage=Val_ChannelVoltage Current=Val_ChannelCurrent

analog_input_configs = Dict{AbstractString,TerminalConfig}(  # deprecate
    "referenced single-ended"     => RSE,
    "non-referenced single-ended" => NRSE,
    "pseudo-differential"         => PseudoDifferential,
    "differential"                => Differential)

"""
`analog_input(channel; terminal_config, range, type=Voltage) -> task`

`analog_input(task, channel; terminal_config, range, type=Voltage)`

create an analog input channel, and a new task if one is not specified.

terminal_config can be Default, RSE, NRSE, Differential, PseudoDifferential
range specifies the minimum and maximum value to measure in volt or ampere
type indicates whether it is a voltage or current input. For current inputs, default it is
 the internal shunt resistor.

The measurements are returned in volt or ampere.

For more information see the NI documentation:
https://zone.ni.com/reference/en-XX/help/370471AM-01/daqmxcfunc/daqmxcreateaivoltagechan/
https://zone.ni.com/reference/en-XX/help/370471AM-01/daqmxcfunc/daqmxcreateaicurrentchan/

"""
function analog_input(channel::String;
                      terminal_config::Union{String,TerminalConfig} = Differential,
                      range = nothing,
                      type  = Voltage
                      )
    if typeof(terminal_config) == String  # deprecate
      Base.depwarn("specifying terminal configurations with Strings is deprecated.  Use the TerminalConfig Enum instead.", :analog_input)
      terminal_config = analog_input_configs[terminal_config]
    end
    t = AITask()
    analog_input(t, channel; terminal_config, range, type)
    t
end

function analog_input(t::AITask,
                      channel::String;
                      terminal_config::Union{String,TerminalConfig}=Differential,
                      range=nothing,
                      type=Voltage)
    if typeof(terminal_config) == String  # deprecate
      Base.depwarn("specifying terminal configurations with Strings is deprecated.  Use the TerminalConfig Enum instead.", :analog_input)
      terminal_config = analog_input_configs[terminal_config]
    end
    if isnothing(range)
        device::String = split(channel,'/')[1]
        if type == Voltage
            range = float(analog_voltage_input_ranges(device)[end,:])
        else
            range = float(analog_current_input_ranges(device)[end,:])
        end
    end
    if type == Voltage # https://zone.ni.com/reference/en-XX/help/370471AM-01/daqmxcfunc/daqmxcreateaivoltagechan/
        catch_error( CreateAIVoltageChan(t.th,
                Ref(codeunits(channel),1),
                Ref(codeunits(""), 1),
                terminal_config,
                range[1], range[2],
                Val_Volts,
                convert(Ptr{UInt8},C_NULL)), "see https://www.ni.com/documentation/en/ni-daqmx/latest/devconsid/defaulttermconfig/" )
    elseif type == Current # https://zone.ni.com/reference/en-XX/help/370471AM-01/daqmxcfunc/daqmxcreateaicurrentchan/
        catch_error( CreateAICurrentChan(t.th,
                Ref(codeunits(channel),1),
                Ref(codeunits(""), 1),
                terminal_config,
                range[1], range[2],
                Val_Amps,
                Val_Default, # shuntResistorLoc
                0 ,
                convert(Ptr{UInt8},C_NULL)), "see https://www.ni.com/documentation/en/ni-daqmx/latest/devconsid/defaulttermconfig/" )
    end

    nothing
end

"""
`acceleration_input(channel, config, range) -> task`

`acceleration_input(task, channel, config, range)`

create an acceleration input channel, and a new task if one is not specified

terminal_config can be Default, RSE, NRSE, Differential, PseudoDifferential
range specifies the minimum and maximum value to measure in g (9.81 m/s^2)
sensitivity specifies the sensor´s sensitivity in mV/g (default is 100mV/g)
currentexcitval is the excitation current in ampere (default 2mA)

The measurements are returned in g.

For more information see the NI documentation:
https://zone.ni.com/reference/en-XX/help/370471AA-01/daqmxcfunc/daqmxcreateaiaccelchan/
"""

function acceleration_input(channel::String;
                      terminal_config::TerminalConfig = Differential,
                      range = nothing, # in g
                      sensitivity::Real = 100. ,
                      currentexcitval::Real = 0.002)

    t = AITask()
    acceleration_input(t, channel; terminal_config, range, sensitivity, currentexcitval)
    t
end

function acceleration_input(t::AITask, channel::String;
                      terminal_config::TerminalConfig = Differential,
                      range = nothing, # in g
                      sensitivity::Real = 100. , # mV / g
                      currentexcitval::Real = 0.002) # Ampere
    if isnothing(range)
        @error "specifiy input ranges"
    end
    # https://zone.ni.com/reference/en-XX/help/370471AA-01/daqmxcfunc/daqmxcreateaiaccelchan/
    catch_error( CreateAIAccelChan(t.th,
            Ref(codeunits(channel),1),
            Ref(codeunits(""), 1),
            terminal_config,
            range[1], range[2],
            Val_AccelUnit_g, # or Val_MetersPerSecondSquared
            sensitivity,
            Val_mVoltsPerG, # or Val_VoltsPerG
            Val_Internal,
            currentexcitval,
            convert(Ptr{UInt8},C_NULL)), "see https://www.ni.com/documentation/en/ni-daqmx/latest/devconsid/defaulttermconfig/" )
    @warn("Attention, this channel has input/output values based in g, not SI (m/s2)")
    nothing
end
"""
`analog_output(channel, range) -> task`

`analog_output(task, channel, range)`

create an analog output channel, and a new NIDAQ task if one is not specified
"""
function analog_output(channel::String; range=nothing)
    t = AOTask()
    analog_output(t, channel, range=range)
    t
end

function analog_output(t::AOTask, channel::String; range=nothing)
    if range == nothing
        device::String = split(channel,'/')[1]
        range=float(analog_output_ranges(device)[end,:])
    end
    catch_error( CreateAOVoltageChan(t.th,
            Ref(codeunits(channel),1),
            Ref(codeunits(""),1),
            range[1], range[2],
            Val_Volts,
            convert(Ptr{UInt8},C_NULL)) )
    nothing
end

read_analog_cfunctions = Dict{Type,Function}(
    Float64 => ReadAnalogF64,
    Int16   => ReadBinaryI16,
    Int32   => ReadBinaryI32,
    UInt16  => ReadBinaryU16,
    UInt32  => ReadBinaryU32 )

function read(t::AITask, num_samples_per_chan::Integer = -1, precision::DataType = Float64)
    outdata_ref = Ref{Cuint}()
    DAQmxGetTaskNumChans(t.th, outdata_ref)
    num_channels = outdata_ref.x
    num_samples_per_chan_read = Int32[0]
    buffer_size = num_samples_per_chan==-1 ? 1024 : num_samples_per_chan
    data = Array{precision}(undef, buffer_size*num_channels)
    catch_error( read_analog_cfunctions[precision](t.th,
        convert(Int32,num_samples_per_chan),
        1.0,
        reinterpret(Bool32,Val_GroupByChannel),
        Ref(data,1),
        convert(UInt32,buffer_size*num_channels),
        Ref(num_samples_per_chan_read,1),
        reinterpret(Ptr{Bool32},C_NULL)) )
    data = data[1:num_samples_per_chan_read[1]*num_channels]
    num_channels==1 ? data : reshape(data, (div(length(data),num_channels), convert(Int64,num_channels)))
end
function Base.read!(data,t::AITask, num_samples_per_chan::Integer = -1, precision::DataType = Float64)
    outdata_ref = Ref{Cuint}()
    DAQmxGetTaskNumChans(t.th, outdata_ref)
    num_channels = outdata_ref.x
    num_samples_per_chan_read = Int32[0]
    buffer_size = num_samples_per_chan==-1 ? 1024 : num_samples_per_chan
    # data = Array{precision}(undef, buffer_size*num_channels)
    catch_error( read_analog_cfunctions[precision](t.th,
        convert(Int32,num_samples_per_chan),
        1.0,
        reinterpret(Bool32,Val_GroupByChannel),
        Ref(data,1),
        convert(UInt32,buffer_size*num_channels),
        Ref(num_samples_per_chan_read,1),
        reinterpret(Ptr{Bool32},C_NULL)) )
end

for (cfunction, types) in (
        (WriteAnalogF64, Float64),
        (WriteBinaryI16, Int16),
        (WriteBinaryI32, Int32),
        (WriteBinaryU16, UInt16),
        (WriteBinaryU32, UInt32))
    @eval function write(t::AOTask, data::Matrix{$types})
        num_samples_per_chan::Int32 = size(data, 1)
        data = reshape(data, length(data))
        num_samples_per_chan_written = Int32[0]
        catch_error( $cfunction(t.th,
            num_samples_per_chan,
            reinterpret(Bool32, UInt32(false)),
            1.0,
            reinterpret(Bool32,Val_GroupByChannel),
            Ref(data,1),
            Ref(num_samples_per_chan_written,1),
            reinterpret(Ptr{Bool32},C_NULL)) )
        num_samples_per_chan_written[1]
    end
    @eval write(t::AOTask, data::Vector{$types}) =
        write(t, reshape(data,(length(data),1)))
end


text:
function insertionsort!(data::Array{Int64})
	for i in 2 : length(data)
		current_value = data[i]
		for j in i - 1 : -1 : 1
			if current_value < data[j]
				data[j + 1] = data[j]
				if j == 1
					data[j] = current_value
				end
			else
				data[j + 1] = current_value
				break
			end
		end
	end
end

function builtinsort!(data::Array{Int64})
	sort!(data)
end

function quicksort!(data::Array, start=-1, stop=-1, lowest=-1, highest=-1)
	if start == -1
		start = 1
		stop = length(data)
		pivot = 0
		lowest = data[1]
		highest = data[1]
		for i in 2 : stop
			if data[i] < lowest
				lowest = data[i]
			end
			if data[i] > highest
				highest = data[i]
			end
		end
	end

	pivot = (highest + lowest) / 2

#=	print(start)
	print(", ")
	print(stop)
	print(", ")
	print(pivot)
	print(", ")
	print(lowest)
	print(", ")
	println(highest)
=#
	if highest == lowest
		return
	end

	if stop - start <= 2
		if stop - start == 2
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
			if data[start + 1] > data[start + 2]
				buffer = data[start + 1]
				data[start + 1] = data[start + 2]
				data[start + 2] = buffer
			end
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
#			println(data)
			return
		end

		if stop - start == 1
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
#			println(data)
			return
		end

		if stop - start == 0
			return
		end
	end

	leftPointer = start
	rightPointer = stop

	leftHighest = lowest
	rightLowest = highest

	done = false

	while !done
#=		print(leftPointer)
		print(", ")
		println(rightPointer)
=#
		checkForLeftExtreme = false
		checkForRightExtreme = false

		left = data[leftPointer] <= pivot
		right = data[rightPointer] > pivot

		if rightPointer - leftPointer <= 0
			done = true
		elseif !left && !right
			buffer = data[leftPointer]
			data[leftPointer] = data[rightPointer]
			data[rightPointer] = buffer
			if rightPointer - leftPointer > 1
				checkForLeftExtreme = true
				checkForRightExtreme = true
				leftPointer += 1
				rightPointer -= 1
			else
				done = true
				if data[leftPointer] > leftHighest
					leftHighest = data[leftPointer]
				end
				if data[rightPointer] < rightLowest
					rightLowest = data[rightPointer]
				end
			end
		elseif !left && right
			rightPointer -= 1
			checkForRightExtreme = true
		elseif left && !right
			leftPointer += 1
			checkForLeftExtreme = true
		elseif rightPointer - leftPointer > 1
			leftPointer += 1
			rightPointer -= 1
			checkForLeftExtreme = true
			checkForRightExtreme = true
		else
			done = true
		end

		if checkForLeftExtreme
			if data[leftPointer - 1] > leftHighest
				leftHighest = data[leftPointer - 1]
			end
		end
		if checkForRightExtreme
			if data[rightPointer + 1] < rightLowest
				rightLowest = data[rightPointer + 1]
			end
		end

#		println(data)
	end

	pivotPointer = leftPointer

	if leftPointer == rightPointer
		if data[leftPointer] > pivot
			pivotPointer -= 1
			if data[leftPointer] < rightLowest
				rightLowest = data[leftPointer]
			end
		else
			if data[leftPointer] > leftHighest
				leftHighest = data[leftPointer]
			end
		end
	else
		if data[leftPointer] > leftHighest
			leftHighest = data[leftPointer]
		end
		if data[rightPointer] < rightLowest
			rightLowest = data[rightPointer]
		end
	end

	quicksort!(data, start, pivotPointer, lowest, leftHighest)
	quicksort!(data, pivotPointer + 1, stop, rightLowest, highest)
end


text:
using RLESUtils, RunUtils, IFTTTUtils

const COMP = gethostname()
const I_START = 1
const I_END = 50
const NOTIFY = false
const VER = "sin"

#uses default config

template_gp(i) =
"""
using ExprSearch,Symbolic_GP,RLESUtils,IFTTTUtils
symbolic_gp(;seed=$i, ver=:$VER, outdir=Pkg.dir("ExprSearch/results/SymbolicCompare_$VER/Symbolic_GP/seed$i"))
$NOTIFY && sendifttt(;value1="gp,$i", value2="$COMP")
"""

template_ge(i) =
"""
using ExprSearch,Symbolic_GE,RLESUtils,IFTTTUtils
symbolic_ge(;seed=$i, ver=:$VER, outdir=Pkg.dir("ExprSearch/results/SymbolicCompare_$VER/Symbolic_GE/seed$i"))
$NOTIFY && sendifttt(;value1="ge,$i", value2="$COMP")
"""

template_mc(i) =
"""
using ExprSearch,Symbolic_MC,RLESUtils,IFTTTUtils
symbolic_mc(;seed=$i, ver=:$VER, outdir=Pkg.dir("ExprSearch/results/SymbolicCompare_$VER/Symbolic_MC/seed$i"))
$NOTIFY && sendifttt(;value1="mc,$i", value2="$COMP")
"""

template_mcts(i) =
"""
using ExprSearch,Symbolic_MCTS,RLESUtils,IFTTTUtils
symbolic_mcts(;seed=$i, ver=:$VER, outdir=Pkg.dir("ExprSearch/results/SymbolicCompare_$VER/Symbolic_MCTS/seed$i"))
$NOTIFY && sendifttt(;value1="mcts,$i", value2="$COMP")
"""

template_ce(i) =
"""
using ExprSearch,Symbolic_CE,RLESUtils,IFTTTUtils
symbolic_ce(;seed=$i, ver=:$VER, outdir=Pkg.dir("ExprSearch/results/SymbolicCompare_$VER/Symbolic_CE/seed$i"))
$NOTIFY && sendifttt(;value1="ce,$i", value2="$COMP")
"""

notifydone() = sendifttt(;value1="$COMP done")

A = JuliaSource[]
append!(A, [JuliaSource(template_ce(i)) for i=I_START:I_END])
append!(A, [JuliaSource(template_ge(i)) for i=I_START:I_END])
append!(A, [JuliaSource(template_gp(i)) for i=I_START:I_END])
append!(A, [JuliaSource(template_mc(i)) for i=I_START:I_END])
#append!(A, [JuliaSource(template_mcts(i)) for i=I_START:I_END])

#include("script.jl")
#fs=julia_process(A,np); notifydone()

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
using Distances
type KMeans
    centers::Matrix{Float64}
    clusters::Vector{Int64}
end

function fit(::Type{KMeans}, X::Matrix{Float64}, k=5, n_iter=20)
    centers = X[:, randcycle(size(X)[2])[1:k]]
    closest = zeros(Int64, size(X)[2])
    dists = zeros(size(X)[2], k)

    for i=1:n_iter
        pairwise!(dists,SqEuclidean(), X, centers)
        closest .= vec(mapslices(indmin,dists,2))
        for j=1:k
            centers[:,j] = mean(X[:, closest .== j], 2)
        end
    end
    return(KMeans(centers, closest))
end


text:
"""
    parse_uints(f, io)

Call function `f(i::UInt, nl::Bool)` for every number `i` in the stream `io`
whereby `nl` tells whether it was the last number in that line
"""
@inline function parse_uints(f::Function, io::IO)
    local i::UInt = 0
    local first::Bool = true
    local finished::Bool = false
    local newline::Bool = false
    while !eof(io)
        local c::UInt8 = read(io, UInt8)
        if c == UInt8(' ') || c == UInt8(',') || c == UInt8('\t')
            finished = true
        elseif c == UInt8('\n') || c == UInt8('\r')
            finished = true
            newline = true
        elseif c == UInt('_')
            nothing
        elseif UInt8('0') <= c <= UInt8('9')
            if finished
                if !first
                    f(i, newline)
                end
                finished = false
                newline = false
                i = 0
            end
            first = false
            i *= 10
            i += c - UInt8('0')
        else
            error("Unknown character: $c")
        end
    end
    if !first
        f(i, true)
    end
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
""" AST and pretty printer for Graphviz's DOT language.

References:

- DOT grammar: http://www.graphviz.org/doc/info/lang.html
- DOT language guide: http://www.graphviz.org/pdf/dotguide.pdf
"""
module Graphviz
export Expression, Statement, Attributes, Graph, Digraph, Subgraph,
  Node, NodeID, Edge, pprint, run_graphviz

using Compat
using DataStructures: OrderedDict
using Parameters: @with_kw_noshow
using Requires: @require

const USE_GV_JLL = Ref(false)

function __init__()
  @require Graphviz_jll="3c863552-8265-54e4-a6dc-903eb78fde85" begin
    USE_GV_JLL[] = true
    let cfg = joinpath(Graphviz_jll.artifact_dir, "lib", "graphviz", "config6")
      if !isfile(cfg)
        Graphviz_jll.dot(path -> run(`$path -c`))
      end
    end
  end
end

# AST
#####

abstract type Expression end
abstract type Statement <: Expression end

""" AST type for Graphviz's "HTML-like" node labels.

For now, the HTML content is just a string.
"""
struct Html
  content::String
end
Base.print(io::IO, html::Html) = print(io, html.content)

const AttributeValue = Union{String,Html}
const Attributes = OrderedDict{Symbol,AttributeValue}

as_attributes(attrs::Attributes) = attrs
as_attributes(d::OrderedDict) = Attributes(Symbol(k) => d[k] for k in keys(d))
as_attributes(d::AbstractDict) =
  Attributes(Symbol(k) => d[k] for k in sort!(collect(keys(d))))

@with_kw_noshow struct Graph <: Expression
  name::String
  directed::Bool
  prog::String="dot"
  stmts::Vector{Statement}=Statement[]
  graph_attrs::Attributes=Attributes()
  node_attrs::Attributes=Attributes()
  edge_attrs::Attributes=Attributes()
end

Graph(name::String, stmts::Vector{Statement}; kw...) =
  Graph(; name=name, directed=false, stmts=stmts, kw...)
Graph(name::String, stmts::Vararg{Statement}; kw...) =
  Graph(; name=name, directed=false, stmts=collect(stmts), kw...)
Digraph(name::String, stmts::Vector{Statement}; kw...) =
  Graph(; name=name, directed=true, stmts=stmts, kw...)
Digraph(name::String, stmts::Vararg{Statement}; kw...) =
  Graph(; name=name, directed=true, stmts=collect(stmts), kw...)

@with_kw_noshow struct Subgraph <: Statement
  name::String="" # Subgraphs can be anonymous
  stmts::Vector{Statement}=Statement[]
  graph_attrs::Attributes=Attributes()
  node_attrs::Attributes=Attributes()
  edge_attrs::Attributes=Attributes()
end

Subgraph(stmts::Vector{Statement}; kw...) = Subgraph(; stmts=stmts, kw...)
Subgraph(stmts::Vararg{Statement}; kw...) = Subgraph(; stmts=collect(stmts), kw...)
Subgraph(name::String, stmts::Vector{Statement}; kw...) =
  Subgraph(; name=name, stmts=stmts, kw...)
Subgraph(name::String, stmts::Vararg{Statement}; kw...) =
  Subgraph(; name=name, stmts=collect(stmts), kw...)

struct Node <: Statement
  name::String
  attrs::Attributes
end
Node(name::String, attrs::AbstractDict) = Node(name, as_attributes(attrs))
Node(name::String; attrs...) = Node(name, attrs)

struct NodeID <: Expression
  name::String
  port::String
  anchor::String
  NodeID(name::String, port::String="", anchor::String="") = new(name, port, anchor)
end

struct Edge <: Statement
  path::Vector{NodeID}
  attrs::Attributes
end
Edge(path::Vector{NodeID}, attrs::AbstractDict) = Edge(path, as_attributes(attrs))
Edge(path::Vector{NodeID}; attrs...) = Edge(path, attrs)
Edge(path::Vararg{NodeID}; attrs...) = Edge(collect(path), attrs)
Edge(path::Vector{String}, attrs::AbstractDict) = Edge(map(NodeID, path), attrs)
Edge(path::Vector{String}; attrs...) = Edge(map(NodeID, path), attrs)
Edge(path::Vararg{String}; attrs...) = Edge(map(NodeID, collect(path)), attrs)

# Bindings
##########

""" Run a Graphviz program.

Invokes Graphviz through its command-line interface. If the `Graphviz_jll`
package is installed and loaded, it is used; otherwise, Graphviz must be
installed on the local system.

For bindings to the Graphviz C API, see the the package
[GraphViz.jl](https://github.com/Keno/GraphViz.jl). At the time of this writing,
GraphViz.jl is unmaintained.
"""
function run_graphviz(io::IO, graph::Graph; prog::Union{String,Nothing}=nothing,
                      format::String="json0")
  if isnothing(prog)
    prog = graph.prog
  end
  @assert prog in ("dot","neato","fdp","sfdp","twopi","circo")
  if USE_GV_JLL[]
    fun = getfield(Graphviz_jll, Symbol(prog))
    prog = fun(identity)
  end
  open(`$prog -T$format`, io, write=true) do gv
    pprint(gv, graph)
  end
end
function run_graphviz(graph::Graph; kw...)
  io = IOBuffer()
  run_graphviz(io, graph; kw...)
  seekstart(io)
end

function Base.show(io::IO, ::MIME"image/svg+xml", graph::Graph)
  run_graphviz(io, graph, format="svg")
end

# Pretty-print
##############

""" Pretty-print the Graphviz expression.
"""
pprint(expr::Expression) = pprint(stdout, expr)
pprint(io::IO, expr::Expression) = pprint(io, expr, 0)

function pprint(io::IO, graph::Graph, n::Int)
  indent(io, n)
  print(io, graph.directed ? "digraph " : "graph ")
  print(io, graph.name)
  println(io, " {")
  pprint_attrs(io, graph.graph_attrs, n+2; pre="graph", post=";\n")
  pprint_attrs(io, graph.node_attrs, n+2; pre="node", post=";\n")
  pprint_attrs(io, graph.edge_attrs, n+2; pre="edge", post=";\n")
  for stmt in graph.stmts
    pprint(io, stmt, n+2, directed=graph.directed)
    println(io)
  end
  indent(io, n)
  println(io, "}")
end

function pprint(io::IO, subgraph::Subgraph, n::Int; directed::Bool=false)
  indent(io, n)
  if isempty(subgraph.name)
    println(io, "{")
  else
    print(io, "subgraph ")
    print(io, subgraph.name)
    println(io, " {")
  end
  pprint_attrs(io, subgraph.graph_attrs, n+2; pre="graph", post=";\n")
  pprint_attrs(io, subgraph.node_attrs, n+2; pre="node", post=";\n")
  pprint_attrs(io, subgraph.edge_attrs, n+2; pre="edge", post=";\n")
  for stmt in subgraph.stmts
    pprint(io, stmt, n+2, directed=directed)
    println(io)
  end
  indent(io, n)
  print(io, "}")
end

function pprint(io::IO, node::Node, n::Int; directed::Bool=false)
  indent(io, n)
  print(io, node.name)
  pprint_attrs(io, node.attrs)
  print(io, ";")
end

function pprint(io::IO, node::NodeID, n::Int)
  print(io, node.name)
  if !isempty(node.port)
    print(io, ":")
    print(io, node.port)
  end
  if !isempty(node.anchor)
    print(io, ":")
    print(io, node.anchor)
  end
end

function pprint(io::IO, edge::Edge, n::Int; directed::Bool=false)
  indent(io, n)
  for (i, node) in enumerate(edge.path)
    if i > 1
      print(io, directed ? " -> " : " -- ")
    end
    pprint(io, node, n)
  end
  pprint_attrs(io, edge.attrs)
  print(io, ";")
end

function pprint_attrs(io::IO, attrs::Attributes, n::Int=0;
                      pre::String="", post::String="")
  if !isempty(attrs)
    indent(io, n)
    print(io, pre)
    print(io, " [")
    for (i, (key, value)) in enumerate(attrs)
      if (i > 1) print(io, ",") end
      print(io, key)
      print(io, "=")
      print(io, value isa Html ? "<" : "\"")
      print(io, value)
      print(io, value isa Html ? ">" : "\"")
    end
    print(io, "]")
    print(io, post)
  end
end

indent(io::IO, n::Int) = print(io, " "^n)

end


text:
module TestQuantizedArrays

using Test
using InteractiveUtils
using LinearAlgebra
using Distances
using QuantizedArrays

include("codebook.jl")
include("quantizer.jl")
include("interface.jl")
include("utils.jl")

end


text:
# The viewer_geometry_data_t has a field named `type`, which is a reserved
# keyword in Julia v0.7 and below. To get around that, we can just
# construct the expression directly instead of relying on the parser
eval(Expr(:struct, true,
    Expr(:(<:), :viewer_geometry_data_t, :LCMType),
    Expr(:block,
        Expr(:(::), :type, :Int8),
        Expr(:(::), :position, :(SVector{3, Float32})),
        Expr(:(::), :quaternion, :(SVector{4, Float32})),
        Expr(:(::), :color, :(SVector{4, Float32})),
        Expr(:(::), :string_data, :String),
        Expr(:(::), :num_float_data, :Int32),
        Expr(:(::), :float_data, :(Vector{Float32})),
)))

@lcmtypesetup(viewer_geometry_data_t,
    float_data => (num_float_data,)
)


text:
function find_minimum(f::Function, itr)
    min_f_value = Inf
    min_element = 0
    min_idx = 0

    for (idx, element) in enumerate(itr)
        f_value = f(element)
        if f_value < min_f_value
            min_f_value = f_value
            min_element = element
            min_idx = idx
        end
    end

    return min_f_value, min_element, min_idx
end

function find_maximum(f::Function, itr)
    max_f_value = -Inf
    max_element = 0
    max_idx = 0

    for (idx, element) in enumerate(itr)
        f_value = f(element)
        if f_value > max_f_value
            max_f_value = f_value
            max_element = element
            max_idx = idx
        end
    end

    return max_f_value, max_element, max_idx
end


text:
module TEST_LINKED_NODE
using Samsara, Conception, Test

@testset "LinkNode initiation" begin
    Conception.__purge_everything!!()

    caseNode = Samsara.LinkedCardinalNode()
    @test isa(caseNode, Samsara.LinkedCardinalNode)
    @test isa(caseNode, Samsara.EuclideanNode)
    @test isa(caseNode, Conception.Time.TemporalType)
    " LinkedCardinalNode lager en struct av type Conception.AbstractConcept --- som også SAT gjør... "

    @test get_id(LinkedCardinalNode("test")) == "test"

    state_set = Conception.MuExS()
    caseNode = Samsara.LinkedCardinalNode(:asdffdsa, in_MuEx=state_set)
    @test Conception.exists_in_MuExS(caseNode, state_set)
    " Lage en LinkedCardinalNode => registreres i arg: [in_MuEx] "

    @test isnothing(caseNode._node_W)
    @test isnothing(caseNode._node_E)
    " Default neste node (i retning E og W) er nothing "

    node1 = LinkedCardinalNode("1")
    node3 = LinkedCardinalNode("3")
    caseNode = LinkedCardinalNode(node_to_E = node3, node_to_W = node1)
    @test isa(caseNode._node_E, Conception.AbstractConcept)
    @test caseNode._node_E == node3
    @test isa(caseNode._node_W, Conception.AbstractConcept)
    @test caseNode._node_W == node1
    " LinkedCardinalNode can have node to West and node to East, set by constructor "

    @test Samsara.east_of(caseNode) == node3
    @test Samsara.west_of(caseNode) == node1
    " convence functions: east of, west of, ... "

    other_node = LinkedCardinalNode()
    @test get_id(LinkedCardinalNode()) != get_id(other_node)
    " Constructed LinkedCardinalNode without id-arg are assigned unique id values "

    node1 = LinkedCardinalNode(:1)
    @assert node1._node_W === nothing
    @assert node1._node_E === nothing
    node2 = LinkedCardinalNode(:2)

    Samsara._set_node_to_W!(node1, node2)
    @test node1._node_W == node2
    " set node to W with function _set_node_to_W!(nodeA, nodeB) "

    Samsara._set_node_to_W!(node1, nothing)
    @test node1._node_W === nothing
    " set node to W with function _set_node_to_W!(nodeA, nothing) "

    Samsara._set_node_to_E!(node2, node1)
    @test node2._node_E == node1
    Samsara._set_node_to_E!(node2, nothing)
    @test node2._node_E === nothing
    " set node to W with function _set_node_to_W!(nodeA, nodeB) "
end

@testset "Conception.TemporalType funksjonalitet" begin
    """ LinkedCardinalNode arber fra Conception.TemporalType, og er meint for å
    kunne bruke funksjonane i Conception.jl/srt/trait.jl
    (sjå også time.jl i Conception.jl om TemporalType)
    """
    muex = Conception.MuExS()
    n1 = LinkedCardinalNode(in_MuEx=muex)
    @test length(muex._elements) == 1
    activate!(n1)
    @test is_active(n1)
    @test muex._active_element == n1._node
    " LinkedCardinalNode fungerer som TemporalType, dvs. likestilt med SAT! "

    muex_stateset = Conception.MuExS()

    #
    # TODO Må være medlem av MuExS! Hopper dit no..
end

@testset "linked_list_factory(N)" begin
    caseList = Samsara.linked_list_factory(1)
    @test length(caseList) == 1
    @test isa(caseList[1], LinkedCardinalNode)
    " factory with N=1 creates list of 1 linked_node "

    Conception.__purge__all_SAT!!()
    caseList = Samsara.linked_list_factory(3)
    @test isa(caseList, Vector{LinkedCardinalNode})
    @test length(caseList) == 3
    " With N = 3    =>     factory makes list of 3 nodes "

    @test caseList[1]._node_E == caseList[2]
    @test caseList[2]._node_W == caseList[1]
    " Element 1 er East for node 2, og element 2 er West for node 1 "

    Conception.__purge__all_SAT!!()
    caseList = Samsara.linked_list_factory(5)
    for i in 2:length(caseList)-1
        @test caseList[i]._node_E == caseList[i+1]
        @test caseList[i]._node_W == caseList[i-1]
    end
    @test caseList[1]._node_E == caseList[2]
    @test caseList[end]._node_W == caseList[end-1]
    " linked_list_factory(N) integration test: creates a doubly-linked list of N elements "

    Conception.__purge_everything!!()
    muex = Conception.MuExS()
    @assert !isnothing(muex)
    caseList = Samsara.linked_list_factory(3, in_MuExS=muex)
    for item in caseList
        @test item._node ∈ muex._elements
    end
    """ LL_factory(N, in_MuExS) lager elementer i muex [in_MuExS] --
    det er enkelt å sette MuExS ved arg i LL_factory """

    for item in caseList
        activate!(item)
        @test Conception.the_active_event_of(muex) == item._node
        " The activating an item causes that item to be active in that muex.. (verifying  "
    end
    for i in 2:length(caseList)
        @test length(caseList[i]._node._incoming_inferences) == 0
    end
    """ .. but the lenght of incoming asscons are 0 (since no concurrent event happened
    while traversing the linked state-list """

    for item in caseList
        activate!(item)
        activate!(SAT(:some_action))
        @test Conception.the_active_event_of(muex) == item._node
    end
    for i in 2:length(caseList)
        @test caseList[i]._node._incoming_inferences[1].nodeA == caseList[i-1]._node
        @test caseList[i]._node._incoming_inferences[1].nodeR == SAT(:some_action)
    end
    """ incoming asscons is sorted, such that the first element is the most recent in history;
    This test verifies that each element has
        * the most recent item in the linked list as nodeL
        * the concurrent event with that node is SAT(:some_action)  (as nodeR)
    """

    for item in caseList
        @show item, item._node._incoming_inferences
    end
    @show caseList
    # TODO lag bra automatiserte testar for forrige opplegget! TODO
end

end #module TEST_LINKED_NODE

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
export gettrajectory

function gettrajectory(Collection,Meta,PArray,ncollisions,fullhistory)
    # Tracks the trajectory of a particle at
    # (x0,y0) with initial (unit speed) velocity (vx,vy)
    # through its next ncollisions collisions.

    edgesx,edgesy  = Meta.edgesx, Meta.edgesy
    nx,ny = length(edgesx)-1, length(edgesy)-1
    lx,ly  = edgesx[end]-edgesx[1], edgesy[end]-edgesy[1]
	nparticles = size(PArray)[1]
	allEvents = Array{Event,1}(undef, 0)
	allPathlengths = zeros(nparticles*ncollisions)
	allCollisionsCounter = 1

	for particleid = 1:nparticles
		P = PArray[particleid]
		collisioncount = 0
		Events = []
		Pathlengths = zeros(ncollisions)

		x0,y0,vx0,vy0,particlerad =  P.x,P.y,P.vx,P.vy,P.r
		# Normalize velocity
		normv = norm([vx0,vy0])
		vx0,vy0 = vx0/normv,vy0/normv

		# Find initial cell of the particle
		i,j  = findindex(x0,edgesx), findindex(y0,edgesy)

		# Make sure the particle is not already inside an obstacle
		if length(Collection[i,j].x)>0
			dists =  sqrt.((x0.-Collection[i,j].x).^2 .+(y0.-Collection[i,j].y).^2) .- Collection[i,j].r
			if minimum(dists) < 0.0  error("A particle is placed inside an obstacle") end
		end

		while collisioncount < ncollisions
			if (collisioncount > 0 )&( mod(collisioncount+1,1000000)==0 )
				println("Collision $(1+collisioncount) of $ncollisions.")
			end
			iscollided = false
			distancesincelastreset = 0.0
			while !iscollided
				# Check if the particle collides in the current cell
				iscollided,distanceincell,newx,newy,newvx,newvy = checkallcollision(
					x0,y0,vx0,vy0,particlerad, Collection[i,j])

				# find out where the particle where to move if it did not collide
				ishift, jshift,distanceincelltrav,newxtrav,newytrav,newvxtrav,newvytrav =traverse(
							x0,y0,vx0,vy0,edgesx[i],edgesx[i+1],edgesy[j],edgesy[j+1])

				if !iscollided
					distanceincell = distanceincelltrav
					newx,newy,newvx,newvy = newxtrav,newytrav, newvxtrav, newvytrav
				else
					# Even if we collided, there is this weird case where a scatterer's center
					# is in the neighboring cell and a particle collides with that scatterer
					# but outside of its current cell
					# This would then not actually be a collisions
					if distanceincell>distanceincelltrav # That is the case!
						iscollided = false
						distanceincell = distanceincelltrav
						newx,newy,newvx,newvy = newxtrav,newytrav, newvxtrav, newvytrav
					else
						collisioncount += 1 # A collision takes place inside the cell
					end
				end

				# Store an event. If a particle "makes use of periodicity", create two events
				# First event without checking for periodicity
				if fullhistory
					e = Event(x0,y0,vx0,vy0,i,j,newx,newy,newvx,newvy,i,j,distanceincell,iscollided,false)
					push!(Events,e)
				end
				# Always store pathlengths if this is the end
				distancesincelastreset += distanceincell

				x0,y0,vx0,vy0 = newx,newy,newvx,newvy

				# Second event: check for periodicity, this is only possible in the absence of a collision
				# actually this is not correct. it can also happy in the case of a collision with a particle
				if !iscollided
				   isperiodic,newx,newy, newi,newj = getperiodic(newx,newy,i,j,ishift,jshift,edgesx,edgesy)
					if isperiodic
						if fullhistory
							e = Event(x0,y0,vx0,vy0,i,j,newx,newy,newvx,newvy,newi,newj,0.0,false,true)
							#println(e)
							push!(Events,e)
						end
					end
					x0,y0,vx0,vy0,i,j = newx,newy,newvx,newvy,newi,newj
				end
			end
			allPathlengths[allCollisionsCounter] = distancesincelastreset
			allCollisionsCounter += 1
			if (x0>edgesx[end]) | (x0<edgesx[1]) | (y0>edgesy[end]) | (y0<edgesy[1])
				error("Out of bounds problem. Log: $x0, $y0, $vx0, $vy0,  $i, $j, $iscollided")
			end
		end
		if fullhistory
			allEvents = [allEvents;Events]
		end
	end
    return allEvents, allPathlengths
end


text:
function modeltest(mdl::ModelInit,method_cube::String="CPU")

    mdl = init_model(mdl)
    mdl.Utils = execute_cube(mdl, method_cube)
    iterate_decission_rules(mdl)
    Ψ_search(mdl)
    new_prices = copy(mdl.MCcondition(mdl))
    push!( mdl.History, tuple(mdl.Prices..., new_prices...) )

    return(mdl)

end
function modeltest(mdl::Model,method_cube::String="CPU")

    mdl = init_model(mdl)
    mdl.Utils = execute_cube(mdl, method_cube)
    iterate_decission_rules(mdl)
    Ψ_search(mdl)
    new_prices = copy(mdl.MCcondition(mdl))
    push!( mdl.History, tuple(mdl.Prices..., new_prices...) )

    return(mdl)

end


text:
"""
    _optimise(approximator::ParametrisedConvexApproximator, x::AbstractVector;
        u_min=nothing, u_max=nothing,
    )

Find a minimiser and optimal value (optval) of `approximator::ParametrisedConvexApproximator` for given
data point `x::AbstractVector`.
Solve with Convex.jl [1].
Available solvers include SCS [2], COSMO [3], Mosek [4], etc.
# Refs.
[1] https://github.com/jump-dev/Convex.jl
[2] https://github.com/jump-dev/SCS.jl
[3] https://github.com/oxfordcontrol/COSMO.jl
[4] https://github.com/MOSEK/Mosek.jl
"""
function _optimise(approximator::ParametrisedConvexApproximator, x::AbstractVector, u_min, u_max;
        solver=SCS,
    )
    @unpack m = approximator
    u = Convex.Variable(m)
    if u_min != nothing
        @assert length(u) == length(u_min)
    end
    if u_max != nothing
        @assert length(u) == length(u_max)
    end
    problem = minimize(approximator(x, u)[1])
    if u_min != nothing
        problem.constraints += [u >= u_min]
    end
    if u_max != nothing
        problem.constraints += [u <= u_max]
    end
    solve!(problem, solver.Optimizer(); verbose=false, silent_solver=true)
    minimiser = typeof(u.value) <: Number ? [u.value] : u.value[:]  # to make it a vector
    optval = [problem.optval]  # to make it a vector
    minimiser, optval
end

"""
    _optimise(approximator::AbstractApproximator, x::AbstractVector;
        u_min=nothing, u_max=nothing,
    )

Find a minimiser and optimal value (optval) of `approximator::AbstractApproximator` (non-parametrised-convex) for given
data point `x::AbstractVector`.
Default solver is `IPNewton` in Optim.jl for box constraints [1].
# Refs.
[1] https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/#box-minimzation
"""
function _optimise(approximator::AbstractApproximator, x::AbstractVector, u_min, u_max)
    @unpack m = approximator
    obj(u) = approximator(x, u)[1]
    if u_min == nothing
        u_min = Float64[]  # no constraint
    end
    if u_max == nothing
        u_max = Float64[]  # no constraint
    end
    # TODO: initial guess; to make sure that `u_min + eps*ones(m) <= u_guess <= u_max - eps*ones(m)`; this algorithm needs an initial guess in the interior of given box constraints.
    u_guess = randn(m)
    if u_min != Float64[]
        u_guess = maximum(hcat(u_min, u_guess); dims=2)[:] + 2*eps()*ones(m)  # make it an interior point
    end
    if u_max != Float64[]
        u_guess = minimum(hcat(u_max, u_guess); dims=2)[:] - eps()*ones(m)  # make it an interior point
    end
    dfc = TwiceDifferentiableConstraints(u_min, u_max)
    res = Optim.optimize(obj, dfc, u_guess, IPNewton())
    minimiser = prod(size(u_guess)) == 1 ? res.minimizer[1] : res.minimizer
    optval = res.minimum
    minimiser, optval
end

function optimise(approximator::AbstractApproximator, x::AbstractVector;
        u_min=nothing, u_max=nothing,
    )
    minimiser, optval = _optimise(approximator, x, u_min, u_max)
    if minimiser == nothing
        @unpack m = approximator
        minimiser = repeat([nothing], m)
    end
    result = (; minimiser=minimiser, optval=optval)
end

"""
    optimise(approximator::AbstractApproximator, x::AbstractMatrix;
        u_min=nothing, u_max=nothing,
    )

Find a minimiser and optimal value (optval) of `approximator::AbstractApproximator` for given
data point `x::AbstractMatrix` using multi-thread computing (powered by Transducers.jl).
"""
function optimise(approximator::AbstractApproximator, x::AbstractMatrix;
        u_min=nothing, u_max=nothing,
        multithreading=true,
    )
    collector = multithreading ? Transducers.tcollect : collect
    d = size(x)[2]
    ress = 1:d |> Map(i -> optimise(approximator, x[:, i]; u_min=u_min, u_max=u_max)) |> collector
    minimiser_matrix = hcat((ress |> Map(res -> res.minimiser) |> collector)...)
    optval_matrix = hcat((ress |> Map(res -> res.optval) |> collector)...)
    (; minimiser = minimiser_matrix, optval = optval_matrix)
end