text:
module Juleva

#include("Gopts.jl")
using ..Gopts

export optimize

# Add the C library to PATH
global const mylib=joinpath(pwd(),"C","libgen.so")
# Arguments
global argv=["-h"]

function f_callback(p::Ptr{Cvoid}, f_::Ptr{Cvoid})
    f = unsafe_pointer_to_objref(f_)::Function
    f(p)::Ptr{Cvoid}
end

function optimize(f::Function)
    return Symbol(f)
end

function optimize_t(f::Function, o::GOpts, x::Vector{Cdouble})
    if (length(x) != o.dim)
        throw(BoundsError())
    end
   f_ = unsafe_pointer_to_objref(f)::Function
    o.x_s = pointer(x)
    f_ptr = @cfunction( $f_
        ,Cdouble,(Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))
    lk=ReentrantLock()
    begin
        lock(lk)
        try

            GC.@preserve f_ptr ccall((:g_optimize,mylib),Int32,(Int32,Ptr{Ptr{UInt8}}
                                                                ,Ptr{Cvoid},Ref{GOpts}),length(argv), argv
                                     ,Base.unsafe_convert(Ptr{Cvoid},f_ptr), Ref(o))
        finally
            unlock(lk)
        end
    end

    return (o.retRes)
end

function optimize(o::GOpts, x::Vector{Cdouble})
    if (length(x) != o.dim)
        throw(BoundsError())
    end
    o.x_s = pointer(x)
    ccall((:g_optimize2,mylib)
          ,Int32,(Int32,Ptr{Ptr{UInt8}}, Ref{GOpts})
          ,length(argv), argv, Ref(o))
    return (unsafe_load(o.retRes), unsafe_load(o.retErr))
end

end

max_stars_repo_name: denisbertini/juleva

max_issues_repo_name: denisbertini/juleva

max_forks_repo_name: denisbertini/juleva


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "iss_model.jl")
include(@relpath "iss_specifications.jl")

S = iss_model()
X0, options = iss_specification()
options = Dict(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# partition with mixed 1D and 135D blocks
partition = vcat([[i] for i in 1:135], [136:270])

# reachability algorithms
algorithm_reach = BFFPSV18(:δ => 1e-3, :vars => [182], :partition => partition,
                           :assume_sparse => true)
algorithm_check = BFFPSV18(:δ => 6e-4, :vars => 136:270, :partition => partition,
                           :assume_sparse => true, :lazy_inputs_interval => -1)

# check property
options[:mode] = "check"
solution = solve(problem, options; op=algorithm_check)
@assert solution.satisfied

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options; op=algorithm_reach)

# project flowpipe to time and x₁₈₂
solution.options[:plot_vars] = [0, 182]
solution_proj = project(solution)

# plot projection
plot(solution_proj)
savefig("iss")

ext: jl

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
export FFJORD

"""
Implementation of FFJORD from

[Grathwohl, Will, Ricky TQ Chen, Jesse Bettencourt, Ilya Sutskever, and David Duvenaud. "Ffjord: Free-form continuous dynamics for scalable reversible generative models." arXiv preprint arXiv:1810.01367 (2018).](https://arxiv.org/abs/1810.01367)
"""
struct FFJORD{T <: AbstractFloat} <: AbstractICNF{T}
    re::Optimisers.Restructure
    p::AbstractVector{T}

    nvars::Integer
    basedist::Distribution
    tspan::Tuple{T, T}

    solver_test::SciMLBase.AbstractODEAlgorithm
    solver_train::SciMLBase.AbstractODEAlgorithm

    sensealg_test::SciMLBase.AbstractSensitivityAlgorithm
    sensealg_train::SciMLBase.AbstractSensitivityAlgorithm

    acceleration::AbstractResource
    array_mover::Function

    # trace_test
    # trace_train
end

function FFJORD{T}(
        nn,
        nvars::Integer,
        ;
        basedist::Distribution=MvNormal(zeros(T, nvars), Diagonal(ones(T, nvars))),
        tspan::Tuple{T, T}=convert(Tuple{T, T}, (0, 1)),

        solver_test::SciMLBase.AbstractODEAlgorithm=default_solver_test,
        solver_train::SciMLBase.AbstractODEAlgorithm=default_solver_train,

        sensealg_test::SciMLBase.AbstractSensitivityAlgorithm=default_sensealg,
        sensealg_train::SciMLBase.AbstractSensitivityAlgorithm=default_sensealg,

        acceleration::AbstractResource=default_acceleration,
        ) where {T <: AbstractFloat}
    array_mover = make_mover(acceleration, T)
    p, re = destructure(nn)
    FFJORD{T}(
        re, p |> array_mover, nvars, basedist, tspan,
        solver_test, solver_train,
        sensealg_test, sensealg_train,
        acceleration, array_mover,
    )
end

function augmented_f(icnf::FFJORD{T}, mode::TestMode)::Function where {T <: AbstractFloat}

    function f_aug(u, p, t)
        m = icnf.re(p)
        z = u[1:end - 1, :]
        mz, J = jacobian_batched(m, z, icnf.array_mover)
        trace_J = transpose(tr.(eachslice(J; dims=3)))
        vcat(mz, -trace_J)
    end
    f_aug
end

function augmented_f(icnf::FFJORD{T}, mode::TrainMode, sz::Tuple{T2, T2}; rng::Union{AbstractRNG, Nothing}=nothing)::Function where {T <: AbstractFloat, T2 <: Integer}
    ϵ = isnothing(rng) ? randn(T, sz) : randn(rng, T, sz)
    ϵ = ϵ |> icnf.array_mover

    function f_aug(u, p, t)
        m = icnf.re(p)
        z = u[1:end - 1, :]
        mz, back = Zygote.pullback(m, z)
        ϵJ = only(back(ϵ))
        trace_J = sum(ϵJ .* ϵ; dims=1)
        vcat(mz, -trace_J)
    end
    f_aug
end

function inference(icnf::FFJORD{T}, mode::TestMode, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractVector where {T <: AbstractFloat}
    xs = xs |> icnf.array_mover
    zrs = zeros(T, 1, size(xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, mode)
    prob = ODEProblem{false}(f_aug, vcat(xs, zrs), icnf.tspan, p)
    sol = solve(prob, icnf.solver_test; sensealg=icnf.sensealg_test)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    Δlogp = fsol[end, :]
    logp̂x = logpdf(icnf.basedist, z) - Δlogp
    logp̂x
end

function inference(icnf::FFJORD{T}, mode::TrainMode, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractVector where {T <: AbstractFloat}
    xs = xs |> icnf.array_mover
    zrs = zeros(T, 1, size(xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, mode, size(xs); rng)
    prob = ODEProblem{false}(f_aug, vcat(xs, zrs), icnf.tspan, p)
    sol = solve(prob, icnf.solver_train; sensealg=icnf.sensealg_train)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    Δlogp = fsol[end, :]
    logp̂x = logpdf(icnf.basedist, z) - Δlogp
    logp̂x
end

function generate(icnf::FFJORD{T}, mode::TestMode, n::Integer, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractMatrix{T} where {T <: AbstractFloat}
    new_xs = isnothing(rng) ? rand(icnf.basedist, n) : rand(rng, icnf.basedist, n)
    new_xs = new_xs |> icnf.array_mover
    zrs = zeros(T, 1, size(new_xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, mode)
    prob = ODEProblem{false}(f_aug, vcat(new_xs, zrs), reverse(icnf.tspan), p)
    sol = solve(prob, icnf.solver_test; sensealg=icnf.sensealg_test)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    z
end

function generate(icnf::FFJORD{T}, mode::TrainMode, n::Integer, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractMatrix{T} where {T <: AbstractFloat}
    new_xs = isnothing(rng) ? rand(icnf.basedist, n) : rand(rng, icnf.basedist, n)
    new_xs = new_xs |> icnf.array_mover
    zrs = zeros(T, 1, size(new_xs, 2)) |> icnf.array_mover
    f_aug = augmented_f(icnf, mode, size(new_xs))
    prob = ODEProblem{false}(f_aug, vcat(new_xs, zrs), reverse(icnf.tspan), p)
    sol = solve(prob, icnf.solver_train; sensealg=icnf.sensealg_train)
    fsol = sol[:, :, end]
    z = fsol[1:end - 1, :]
    z
end

Flux.@functor FFJORD (p,)

function loss(icnf::FFJORD{T}, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; agg::Function=mean) where {T <: AbstractFloat}
    logp̂x = inference(icnf, TrainMode(), xs, p)
    agg(-logp̂x)
end


text:
"""
    dfplsr_cg(X, y; nlv, reorth = true)
Compute the model complexity (df) of PLSR models with the CGLS algorithm.
* `X` : X-data.
* `y` : Univariate Y-data.
* `nlv` : Nb. latent variables (LVs).
* `reorth` : If `true`, a Gram-Schmidt reorthogonalization of the normal equation
    residual vectors is done.

The number of degrees of freedom (df) of the model is returned for 0, 1, ..., nlv LVs.

`X` and `y` are internally centered. The models are computed with an intercept.

## References

Hansen, P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems, Mathematical Modeling and Computation.
Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719697

Hansen, P.C., 2008. Regularization Tools version 4.0 for Matlab 7.3.
Numer Algor 46, 189–194. https://doi.org/10.1007/s11075-007-9136-9

Lesnoff, M., Roger, J.M., Rutledge, D.N., Submitted. Monte Carlo methods for estimating
Mallows’s Cp and AIC criteria for PLSR models. Illustration on agronomic spectroscopic NIR data.
Journal of Chemometrics.
"""
function dfplsr_cg(X, y; nlv, reorth = true)
    F = cglsr(X, y; nlv = nlv, reorth = reorth, filt = true).F
    df = [1 ; vec(1 .+ sum(F, dims = 1))]
    (df = df,)
end

"""
    aicplsr(X, y; nlv, reorth = true, filt = false)
Compute the Akaike's (AIC) and Mallows's (Cp) criteria for univariate PLSR models.
* `X` : X-data.
* `y` : Univariate Y-data.
* `nlv` : Nb. latent variables (LVs).
* `reorth` : If `true`, a Gram-Schmidt reorthogonalization of the normal equation
    residual vectors is done.

For the methodolgy, see Hansen 1998, 2008 and Lesnoff et al. Submitted.

`X` and `y` are internally centered. The models are computed with an intercept.

## References

Hansen, P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems, Mathematical Modeling and Computation.
Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719697

Hansen, P.C., 2008. Regularization Tools version 4.0 for Matlab 7.3.
Numer Algor 46, 189–194. https://doi.org/10.1007/s11075-007-9136-9

Lesnoff, M., Roger, J.M., Rutledge, D.N., Submitted. Monte Carlo methods for estimating
Mallows’s Cp and AIC criteria for PLSR models. Illustration on agronomic spectroscopic NIR data.
Journal of Chemometrics
"""
function aicplsr(X, y; nlv, correct = true)
    X = ensure_mat(X)
    n = size(X, 1)
    p = size(X, 2)
    nlv = min(nlv, n, p)
    res = gridscorelv(X, y, X, y;
        fun = plskern, score = ssr, nlv = 0:nlv)
    zssr = res.y1
    df = dfplsr_cg(X, y, nlv = nlv, reorth = true).df
    df_ssr = n .- df
    # For Cp, unbiased estimate of sigma2
    # ----- Cp1: From a low biased model
    # Not stable with dfcov and nlv too large compared to best model !!
    # If df stays below .95 * n, this corresponds
    # to the maximal model (nlv)
    # Option 2 gives in general results
    # very close to those of option 1,
    # but can give poor results with dfcov
    # when nlv is set too large to the best model
    k = maximum(findall(df .<= .5 * n))
    s2_1 = zssr[k] / df_ssr[k]
    # ----- Cp2: FPE-like
    # s2 is estimated from the model under evaluation
    # Used in Kraemer & Sugiyama 2011 Eq.5-6
    s2_2 = zssr ./ df_ssr
    ct = ones(nlv + 1)
    correct ? ct .= n ./ (n .- df .- 2) : nothing
    ct[(df .> n) .| (ct .<= 0)] .= NaN
    # For safe predictions when df stabilizes and fluctuates
    ct[df .> .8 * n] .= NaN
    # End
    u = findall(isnan.(ct))
    if length(u) > 0
        ct[minimum(u):(nlv + 1)] .= NaN
    end
    aic = n * log.(zssr) + 2 * (df .+ 1) .* ct
    cp1 = zssr .+ 2 * s2_1 * df .* ct
    cp2 = zssr .+ 2 * s2_2 .* df .* ct
    cp1 = cp1 / n
    cp2 = cp2 / n
    res = (aic = aic, cp1 = cp1, cp2 = cp2)
    znlv = 0:nlv
    tab = DataFrame(nlv = znlv, n = fill(n, nlv + 1), df = df, ct = ct, ssr = zssr)
    crit = hcat(tab, DataFrame(res))
    opt = map(x -> findmin(x[isnan.(x) .== 0])[2] - 1, res)
    delta = map(x -> x .- findmin(x[isnan.(x) .== 0])[1], res)                  # Differences "Delta"
    w = map(x -> exp.(-x / 2) / sum(exp.(-x[isnan.(x) .== 0] / 2)), delta)      # AIC odel weights
    delta = reduce(hcat, delta)
    w = reduce(hcat, w)
    nam = [:aic, :cp1, :cp2]
    delta = DataFrame(delta, nam)
    insertcols!(delta, 1, :nlv => znlv)
    w = DataFrame(w, nam)
    insertcols!(w, 1, :nlv => znlv)
    (crit = crit, opt = opt, delta = delta, w = w)
end


text:
"""
	mutable struct GLMesh

Container of (a) the `Int32` code of primitives in `GLMesh`; (b) a matrix transform in local space, the `GLVertexArray`, and three objects `GLVertexBuffer` for `vertices`, `normals`, and `colors`. Inner constructors for initialization.

# Example
Generation of the object `GLMesh` modeling the unit 3D cube with one vertex on the origin:
```
julia> GL.GLCuboid(GL.Box3d(GL.Point3d(0,0,0),GL.Point3d(1,1,1)))
ViewerGL.GLMesh(4, [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0], ViewerGL.GLVertexArray(-1), ViewerGL.GLVertexBuffer(-1, Float32[0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0  …  1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0]), ViewerGL.GLVertexBuffer(-1, Float32[0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0  …  0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0]), ViewerGL.GLVertexBuffer(-1, Float32[]))
```
"""
mutable struct GLMesh

	primitive::Int32
	T::Matrix4
	vertex_array::GLVertexArray

	vertices::GLVertexBuffer
	normals::GLVertexBuffer
	colors::GLVertexBuffer

	# constructor
	function GLMesh()
		ret=new(POINTS,Matrix4(),GLVertexArray(),GLVertexBuffer(),GLVertexBuffer(),GLVertexBuffer())
		finalizer(releaseGpuResources, ret)
		return ret
	end

	# constructor
	function GLMesh(primitive)
		ret=new(primitive,Matrix4(),GLVertexArray(),GLVertexBuffer(),GLVertexBuffer(),GLVertexBuffer())
		finalizer(releaseGpuResources, ret)
		return ret
	end

end

"""
	releaseGpuResources(mesh::GLMesh)

Release the memory allocated by the `vertex_array`, `vertices`, `normals`, and `colors` on the GPU.
"""
function releaseGpuResources(mesh::GLMesh)
	releaseGpuResources(mesh.vertex_array)
	releaseGpuResources(mesh.vertices)
	releaseGpuResources(mesh.normals)
	releaseGpuResources(mesh.colors)
end

"""
	computeNormal(p1::Point2d, p2::Point2d)

Compute in 2D the ccw normal vector to a line segment oriented from `p1` to `p2`.

# Example
```
julia> GL.computeNormal(GL.Point2d(1,0), GL.Point2d(0.8,0.5))
2-element StaticArrays.MArray{Tuple{2},Float64,1,2}:
 0.9284766908852594
 0.37139067635410367

julia> GL.computeNormal(GL.Point2d(-0.8,0.5), GL.Point2d(-1,0))
2-element StaticArrays.MArray{Tuple{2},Float64,1,2}:
 -0.9284766908852594
  0.37139067635410367
```
"""
function computeNormal(p1::Point2d, p2::Point2d)
	t = p2-p1
	n = Point2d(+t[2], -t[1])
	return normalized(n)
end

"""
	computeNormal(p0::Point3d,p1::Point3d,p2::Point3d)

Compute the ccw normal vector to a 3D triangle with orientation `(p0,p1,p2)`.
# Example
```
julia> GL.computeNormal(GL.Point3d(0,0,0), GL.Point3d(1,0,0),GL.Point3d(0,1,0))
3-element StaticArrays.MArray{Tuple{3},Float64,1,3}:
 0.0
 0.0
 1.0
```
"""
function computeNormal(p0::Point3d,p1::Point3d,p2::Point3d)
	return normalized(cross(p1-p0,p2-p0))
end

"""
	getBoundingBox(mesh::GLMesh)

Return an object of ViewerGL.Box3d
# Example
```
julia> mesh = GL.GLCuboid(GL.Box3d(GL.Point3d(0,0,0),GL.Point3d(1,1,1)));

julia> GL.getBoundingBox(mesh)
ViewerGL.Box3d([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
```
"""
function getBoundingBox(mesh::GLMesh)
	box=invalidBox()
	vertices=mesh.vertices.vector
	for I in 1:3:length(vertices)
		point=Point3d(vertices[I+0],vertices[I+1],vertices[I+2])
		addPoint(box,point)
	end
	return box
end
function getBoundingBox(meshes::Array)
	box=invalidBox()
	for mesh in meshes
		vertices=mesh.vertices.vector
		for I in 1:3:length(vertices)
			point=Point3d(vertices[I+0],vertices[I+1],vertices[I+2])
			addPoint(box,point)
		end
	end
	return box
end

"""
	GLCuboid(box::Box3d)

Return the `GLMesh` modeling a 3-cuboid, i.e. the 3D parallelepiped parallel to the axes of reference frame).
# Example
```
julia> mesh = GL.GLAxis(GL.Point3d(-2,-2,-2),GL.Point3d(+2,+2,+2));

julia> box = GL.getBoundingBox(mesh)
ViewerGL.Box3d([-2.0, -2.0, -2.0], [2.0, 2.0, 2.0])

julia> GL.GLCuboid(box)
ViewerGL.GLMesh(4, [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0], ViewerGL.GLVertexArray(-1), ViewerGL.GLVertexBuffer(-1, Float32[-2.0, 2.0, -2.0, 2.0, 2.0, -2.0, 2.0, -2.0, -2.0, -2.0  …  2.0, -2.0, 2.0, -2.0, -2.0, -2.0, 2.0, -2.0, 2.0, 2.0]), ViewerGL.GLVertexBuffer(-1, Float32[0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0  …  0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0, -1.0, 0.0, 0.0]), ViewerGL.GLVertexBuffer(-1, Float32[]))
```
"""
function GLCuboid(box::Box3d)

	points=getPoints(box)
	faces=[[1, 2, 3, 4],[4, 3, 7, 8],[8, 7, 6, 5],[5, 6, 2, 1],[6, 7, 3, 2],[8, 5, 1, 4]]

	vertices=Vector{Float32}()
	normals =Vector{Float32}()
	for face in faces

		p3,p2,p1,p0 = points[face[1]],points[face[2]],points[face[3]],points[face[4]] # reverse order
		n=0.5*(computeNormal(p0,p1,p2) + computeNormal(p0,p2,p3))

		append!(vertices,p0); append!(normals,n)
		append!(vertices,p1); append!(normals,n)
		append!(vertices,p2); append!(normals,n)
		append!(vertices,p0); append!(normals,n)
		append!(vertices,p2); append!(normals,n)
		append!(vertices,p3); append!(normals,n)
	end

	ret=GLMesh(GL_TRIANGLES)
	ret.vertices = GLVertexBuffer(vertices)
	ret.normals  = GLVertexBuffer(normals)
	return ret
end

"""
	GLAxis(p0::Point3d,p1::Point3d)

Generate the GLMesh of a reference frame in local space.
Axes *position* and *length* is determined by the bounding box of `p1` and `p2`.
The reference frame is colored according to the *`RGB` colors*, respectively.

# Example
Note the size and position of reference frame in locol coordinates ...
```
GL.VIEW([
	GL.GLCuboid(GL.Box3d(GL.Point3d(0,0,0),GL.Point3d(1,1,1)))
	GL.GLAxis(GL.Point3d(-2,-2,-2),GL.Point3d(+2,+2,+2))
	])
```
"""
function GLAxis(p0::Point3d,p1::Point3d)

	vertices=Vector{Float32}()
	colors  =Vector{Float32}()

	R=Point4d(1,0,0,1); append!(vertices,p0); append!(vertices,Point3d(p1[1],p0[2],p0[3])); append!(colors,R); append!(colors,R)
	G=Point4d(0,1,0,1); append!(vertices,p0); append!(vertices,Point3d(p0[1],p1[2],p0[3])); append!(colors,G); append!(colors,G)
	B=Point4d(0,0,1,1); append!(vertices,p0); append!(vertices,Point3d(p0[1],p0[2],p1[3])); append!(colors,B); append!(colors,B)

	ret=GLMesh(GL_LINES)
	ret.vertices=GLVertexBuffer(vertices)
	ret.colors  =GLVertexBuffer(colors)
	return ret
end


text:
module Psycho

export
    standardize,
    standardize!,
    perfectNormal,

    r2_tjur,
    datagrid,

    Rules,
    interpret,
    interpret_p,
    format_p,
    interpret_d,
    interpret_odds,

    simulate_groupnames,
    simulate_coefs_correlation,
    simulate_data_correlation,
    simulate_data_logistic,

    sdt_indices,

    Report,
    report,
    table,
    values

include("core/standardize.jl")
include("core/perfectNormal.jl")

include("modelling/R2.jl")
include("modelling/datagrid.jl")

include("interpret/interpret.jl")
include("interpret/pvalue.jl")
include("interpret/cohen_d.jl")
include("interpret/odds.jl")

include("simulate/data_groups.jl")
include("simulate/data_correlation.jl")
include("simulate/data_logistic.jl")

include("sdt/GLM.jl")
include("sdt/basics.jl")

include("report/report.jl")
include("report/Vector.jl")
include("report/DataFrame.jl")
include("report/LM.jl")
include("report/GLM.jl")

end # module


text:
using ReinforcementLearning
using StableRNGs
using OpenSpiel
using Flux

mutable struct KuhnOpenEDHook <: AbstractHook
    results::Vector{Float64}
end

function (hook::KuhnOpenEDHook)(::PreEpisodeStage, policy, env)
    # get nash_conv of the current policy.
    push!(hook.results, RLZoo.nash_conv(policy, env))

    # update agents' learning rate.
    for (_, agent) in policy.agents
        agent.learner.optimizer[2].eta = 1.0 / sqrt(length(hook.results))
    end
end

function (hook::KuhnOpenEDHook)(::PostExperimentStage, policy, env)
    reset!(env)

    # get nash_conv of the latest model.
    push!(hook.results, RLZoo.nash_conv(policy, env))
end

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:ED},
    ::Val{:OpenSpiel},
    game;
    seed = 123,
)
    rng = StableRNG(seed)

    env = OpenSpielEnv(game)
    wrapped_env = ActionTransformedEnv(
        env,
        action_mapping = a -> RLBase.current_player(env) == chance_player(env) ? a : Int(a - 1),
        action_space_mapping = as -> RLBase.current_player(env) == chance_player(env) ?
            as : Base.OneTo(num_distinct_actions(env.game)),
    )
    wrapped_env = DefaultStateStyleEnv{InformationSet{Array}()}(wrapped_env)
    player = 0 # or 1
    ns, na = length(state(wrapped_env, player)), length(action_space(wrapped_env, player))

    create_network() = Chain(
        Dense(ns, 64, relu;init = glorot_uniform(rng)),
        Dense(64, na;init = glorot_uniform(rng))
    )

    create_learner() = NeuralNetworkApproximator(
        model = create_network(),
        optimizer = Flux.Optimise.Optimiser(WeightDecay(0.001), Descent())
    )

    EDmanager = EDManager(
        Dict(
            player => EDPolicy(
                1 - player, # opponent
                create_learner(), # neural network learner
                WeightedSoftmaxExplorer(), # explorer
            ) for player in players(env) if player != chance_player(env)
        )
    )

    stop_condition = StopAfterEpisode(500, is_show_progress=!haskey(ENV, "CI"))
    hook = KuhnOpenEDHook([])

    Experiment(EDmanager, wrapped_env, stop_condition, hook, "# play OpenSpiel $game with ED algorithm")
end

using Plots
ex = E`JuliaRL_ED_OpenSpiel(kuhn_poker)`
run(ex)
plot(ex.hook.results, xlabel="episode", ylabel="nash_conv")

savefig("assets/JuliaRL_ED_OpenSpiel(kuhn_poker).png")#hide

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
#=
Popis
  Funkcia prepise obsah datovej struktury grafu do suboru.

Parametre
  g : datova struktura grafu :: GenericGraph
  path : absolutna cesta k suboru kam sa ulozi graf :: String
  dist_key : nazov atributu, v ktorom je ulozene ohodnotenie hrany :: String
  index : parameter pre zmenu sposobu ulozenia vrcholov :: Bool

Navratova hodnota
  Funkcia nema ziadnu navratovu hodnotu.

Poznamky
  Funkcia vytvori subor alebo ak uz v danom adresari subor
  existuje tak je jeho obsah prepisany. Parameter index umoznuje pouzivatelovi
  prepnut medzi zapisom vrcholov do suboru podla indexu alebo ich oznacenia tzv. labelu.
  Parameter dist_key oznacuje nazov atributu kde je ulozene ohodnotenie hrany.
=#

function save_graph(g::GenericGraph, path::AbstractString; dist_key::AbstractString="weight", index::Bool=true)

  edge_array = edges(g)
  vertex_array = vertices(g)
  file = nothing

  try
    file = open(path, "w")
  catch e
    println(e)
  end

  if index

    for edge in edge_array

      write(file, dec(edge.source.index), " ",
            dec(edge.target.index), " ",
            string(edge.attributes[dist_key]), "\r\n")

	  end

  else

    for edge in edge_array

      write(file, edge.source.label, " ",
            edge.target.label, " ",
            string(edge.attributes[dist_key]), "\r\n")

	  end

  end

  try
    close(file)
  catch e
    println(e)
  end

  return

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
module MangalTestConversions
    using EcologicalNetworks
    using Mangal
    using Test

    # random network
    N = first(networks("name" => "hadfield_2014_20140201_1007"))

    @test typeof(convert(UnipartiteNetwork, N)) <: UnipartiteNetwork

end


text:
## Mimi UI
using VegaLite
using FilePaths
import FileIO:save
export save

global app = nothing

#include functions and modules
include("getdataparts.jl")
include("buildspecs.jl")
include("buildmenu.jl")
include("results.jl")

"""
    explore(m::Model)

Produce a UI to explore the parameters and variables of Model `m` in an independent window.
"""
function explore(m::Model)

    if m.mi === nothing
        error("A model must be run before it can be plotted")
    end

    #start Electron app
    if app === nothing
        global app = Application()
    end

    #window options
    windowopts = Dict("width" => 1000, "height" => 1000)
    w = Window(app, joinpath(@__PATH__, p"mimiexplorer-app/build/index.html"), options = windowopts)

    #set async block to process messages
    @async for msg in msgchannel(w)
        if (msg["cmd"] == "display_spec")
            spec = _spec_for_item(m, Symbol(msg["comp_name"]), Symbol(msg["item_name"]))
            specJSON = JSON.json(spec)
            run(w, "display($specJSON)")
        end
        if (msg["cmd"] == "update_data")
            comp_name = msg["comp_name"];
            paths = _get_all_paths(m)
            comp_path = paths[Symbol(comp_name)];
            comp_def = find_comp(m, comp_path);
            menulist = menu_item_list(m, Symbol(comp_name), comp_def)
            menulistJSON = JSON.json(menulist);
            result = run(w, "setData($menulistJSON)");
        end
    end

    # Electron.toggle_devtools(w)

    #refresh tree view
    subcomplist = tree_view_values(m)
    subcomplistJSON = JSON.json(subcomplist)

    result = run(w, "setTreeChildren($subcomplistJSON)")

    #refresh data view
    datalist = menu_item_list(m)
    datalistJSON = JSON.json(datalist)
    result = run(w, "setData($datalistJSON)")

    return w

end

function explore(mi::ModelInstance)
    m = Model(mi)
    m.md.dirty = false # we need this to get explorer working, but it's a hack and should be temporary!
    explore(m)
end

"""
    explore(sim_inst::SimulationInstance; title="Electron", model_index::Int = 1, scen_name::Union{Nothing, String} = nothing, results_output_dir::Union{Nothing, String} = nothing)

Produce a UI to explore the output distributions of the saved variables in `SimulationInstance`
`sim` for results of model `model_index` and scenario with the name `scen_name` in
a Window with title `title`. The optional arguments default to a `model_index` of `1`, a `scen_name` of `nothing`
assuming there is no secenario dimension, and a window with title `Electron`.
The `results_output_dir` keyword argument refers to the main output directory as provided to `run`,
where all subdirectories are held. If provided, results are assumed to be stored there, otherwise it is
assumed that results are held in results.sim and not
in an output folder.

"""
function explore(sim_inst::SimulationInstance; title="Electron", model_index::Int = 1, scen_name::Union{Nothing, String} = nothing, results_output_dir::Union{Nothing, String} = nothing)

    # quick check
    if results_output_dir === nothing && isempty(sim_inst.results[model_index])
        error("Simulation instance results dictionaries are empty, if results were only saved to disk use the `results_output_dir` keyword argument.")
    end

    #start Electron app
    if app === nothing
        global app = Application()
    end

    #load main html file
    mainpath = replace(joinpath(@__DIR__, "assets", "main.html"), "\\" => "/")

    #window options
    windowopts = Dict("title" => title, "width" => 1100, "height" => 700)
    w = Window(app, URI(joinpath(@__PATH__, "assets", "main.html")), options = windowopts)

     #set async block to process messages
     @async for msg in msgchannel(w)
        comp_name = Symbol(msg["comp_name"]);
        item_name = Symbol(msg["item_name"]);

        if results_output_dir === nothing
            results = Mimi.get_sim_results(sim_inst, comp_name, item_name; model_index = model_index, scen_name = scen_name)
        else
            results = Mimi.get_sim_results(sim_inst, comp_name, item_name, results_output_dir; model_index = model_index, scen_name = scen_name)
        end

        spec = Mimi._spec_for_sim_item(sim_inst, comp_name, item_name, results; model_index = model_index)
        specJSON = JSON.json(spec)

        run(w, "display($specJSON)")

    end

     #refresh variable list
     menulist = menu_item_list(sim_inst)
     menulistJSON = JSON.json(menulist)

     result = run(w, "refresh($menulistJSON)")

     return w

end

# Helper function returns true if VegaLite is verison 3 or above, and false otherwise
function _is_VegaLite_v3()
    return isdefined(VegaLite, :vlplot) ? true : false
end

"""
    plot(m::Model, comp_name::Symbol, datum_name::Symbol; interactive::Bool = false)

Plot a specific `datum_name` (a `variable` or `parameter`) of Model `m`. If the
Bool `interactive` option is set to `true`, the plot will be interactive which will
limit saving options.
"""
function plot(m::Model, comp_name::Symbol, datum_name::Symbol; interactive::Bool = false)

    if m.mi === nothing
        error("A model must be run before it can be plotted")
    end

    spec = Mimi._spec_for_item(m, comp_name, datum_name, interactive=interactive)
    spec === nothing ? error("Spec cannot be built.") : VLspec = spec["VLspec"]

    return _is_VegaLite_v3() ? VegaLite.VLSpec(VLspec) : VegaLite.VLSpec{:plot}(VLspec)
end
"""
    plot(sim_inst::SimulationInstance, comp_name::Symbol, datum_name::Symbol; interactive::Bool = false, model_index::Int = 1, scen_name::Union{Nothing, String} = nothing, results_output_dir::Union{Nothing, String} = nothing)

Plot a specific `datum_name` that was one of the saved variables of `SimulationInstance` `sim_inst`
for results of model `model_index`, which defaults to 1, in scenario `scen_name`, which
defaults to `nothing` implying there is no scenario dimension. The `results_output_dir` keyword argument refers
to the main output directory as provided to `run`, where all subdirectories are held. If provided, results are
assumed to be stored there, otherwise it is assumed that results are held in results.sim and not
in an output folder.
"""
function plot(sim_inst::SimulationInstance, comp_name::Symbol, datum_name::Symbol; interactive::Bool = false, model_index::Int = 1, scen_name::Union{Nothing, String} = nothing, results_output_dir::Union{Nothing, String} = nothing)

    # quick check
    if results_output_dir === nothing && isempty(sim_inst.results[model_index])
        error("Simulation instance results dictionaries are empty, if results were only saved to disk use the `results_output_dir` keyword argument.")
    end

    if results_output_dir === nothing
        results = Mimi.get_sim_results(sim_inst, comp_name, datum_name; model_index = model_index, scen_name = scen_name)
    else
        results = Mimi.get_sim_results(sim_inst, comp_name, datum_name, results_output_dir; model_index = model_index, scen_name = scen_name)
    end

    spec = Mimi._spec_for_sim_item(sim_inst, comp_name, datum_name, results; interactive = interactive, model_index = model_index)
    spec === nothing ? error("Spec cannot be built.") : VLspec = spec["VLspec"]

    return _is_VegaLite_v3() ? VegaLite.VLSpec(VLspec) : VegaLite.VLSpec{:plot}(VLspec)

end


text:
__precompile__(true)

module ThermostattedDynamics

using ODE
using HDF5
using PyCall
using Cubature
using ForwardDiff
using StatsBase

import Distributions.Normal
import Distributions.Logistic
import ForwardDiff.derivative

export Thermostat, Potential, Integrator, Parameters, lyapunov_exponents, trajectory, marginaldistributions, hellingerdistance, sections,
save_lyapunov, save_trajectory, section

const sm = PyNULL()

function __init__()
    copy!(sm, pyimport_conda("statsmodels.api", "statsmodels"))
end

include("DensityDynamics/randominitialcondition.jl")
include("DensityDynamics/DDtypes.jl")
include("DensityDynamics/DDmethods.jl")
include("DensityDynamics/lyapunov.jl")
include("DensityDynamics/trajectory.jl")
include("DensityDynamics/runsimulations.jl")
include("analysis/marginaldistributions.jl")
include("analysis/hellingerdistance.jl")
include("analysis/sections.jl")

end

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
# verify transformation semigroup generated by wrappers
#
REAL_ID = [4 5 6; 3 4 5; 2 3 4]
COMPLEX_ID = REAL_ID * (1 + im)

#identity
#transpose
#adjoint
symmup(A) = Symmetric(A, :U)
symmlo(A) = Symmetric(A, :L)
hermup(A) = Hermitian(A, :U)
hermlo(A) = Hermitian(A, :L)
triaup(A) = UpperTriangular(A)
trialo(A) = LowerTriangular(A)
triaupu(A) = UnitUpperTriangular(A)
trialou(A) = UnitLowerTriangular(A)
diagonal(A) = Diagonal(A)
bidiup(A) = Bidiagonal(A, :U)
bidilo(A) = Bidiagonal(A, :L)
tridi(A) = Tridiagonal(A)
symmtridiup(A) = SymTridiagonal(Symmetric(A, :U))
symmtridilo(A) = SymTridiagonal(Symmetric(A, :L))
hermtridiup(A) = Hermitian(Tridiagonal(A), :U)
hermtridilo(A) = Hermitian(Tridiagonal(A), :L)

REAL_GEN = [identity, transpose, symmup, symmlo, triaup, trialo, triaupu, trialou,
            diagonal, bidiup, bidilo, tridi, symmtridiup, symmtridilo]

COMPLEX_GEN = [REAL_GEN..., adjoint, hermup, hermlo, hermtridiup, hermtridilo]

function closegen(ID::T, GEN::Vector) where T<:AbstractMatrix
    dic = Dict{T,Any}(ID => Set{Any}())
    modified = true
    len(set::Set) = isempty(set) ? 0 : length(first(set))
    setkg(k, gen) = isempty(dic[k]) ? (gen,) : (gen, first(dic[k])...)
    while modified
        modified = false
        keysprev = collect(keys(dic))
        for gen in GEN
            for k in keysprev
                gk = Matrix(gen(k))
                if !haskey(dic, gk)
                    dic[gk] = Set{Any}(Ref(setkg(k, gen)))
                    modified = true
                elseif len(dic[gk]) == len(dic[k]) + 1
                    push!(dic[gk], setkg(k, gen))
                end
            end
        end
    end
    dic
end

closereal() = closegen(REAL_ID, REAL_GEN)
closecomplex() = closegen(COMPLEX_ID, COMPLEX_GEN)


text:
@testset "cell" begin
	# 1D crystals
	cell = UnitCell(1, :lamellar, 10.0)
	@test cell.volume ≈ 10.0

	# 2D crystals
	cell = UnitCell(2, :square, 10.0)
	@test cell.volume ≈ 100.0

	cell = UnitCell(2, :rectangular, 10.0, 5.0)
	@test cell.volume ≈ 50.0

	cell = UnitCell(2, :hexagonal, 10.0)
	@test cell.volume ≈ 86.60254037844385

	cell = UnitCell(2, :oblique, 10.0, 5.0, pi/4)
	@test cell.volume ≈ 35.35533905932737

	# 3D crystals
	cell = UnitCell(3, :cubic, 10.0)
	@test cell.volume ≈ 1000.0

	cell = cell = UnitCell(3, :tetragonal, 10.0, 5.0)
	@test cell.volume ≈ 500.0

	cell = UnitCell(3, :orthorhombic, 10.0, 5.0, 2.0)
	@test cell.volume ≈ 100.0

	cell = UnitCell(3, :hexagonal, 10.0, 5.0)
	@test cell.volume ≈ 433.01270189221924

	cell = UnitCell(3, :monoclinic, 10.0, 5.0, 2.0, pi/4)
	@test cell.volume ≈ 70.71067811865474

	cell = UnitCell(3, :triclinic, 10.0, 5.0, 2.0, pi/2, pi/2, pi/2)
	@test cell.volume ≈ 100.0
end


text:
export lsh

import ..LSHSearch

using mlpack._Internal.io

import mlpack_jll
const lshLibrary = mlpack_jll.libmlpack_julia_lsh

# Call the C binding of the mlpack lsh binding.
function lsh_mlpackMain()
  success = ccall((:lsh, lshLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module lsh_internal
  import ..lshLibrary

import ...LSHSearch

# Get the value of a model pointer parameter of type LSHSearch.
function IOGetParamLSHSearch(paramName::String)::LSHSearch
  LSHSearch(ccall((:IO_GetParamLSHSearchPtr, lshLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type LSHSearch.
function IOSetParamLSHSearch(paramName::String, model::LSHSearch)
  ccall((:IO_SetParamLSHSearchPtr, lshLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeLSHSearch(stream::IO, model::LSHSearch)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeLSHSearchPtr, lshLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeLSHSearch(stream::IO)::LSHSearch
  buffer = read(stream)
  LSHSearch(ccall((:DeserializeLSHSearchPtr, lshLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    lsh(; [bucket_size, hash_width, input_model, k, num_probes, projections, query, reference, second_hash_size, seed, tables, true_neighbors, verbose])

This program will calculate the k approximate-nearest-neighbors of a set of
points using locality-sensitive hashing. You may specify a separate set of
reference points and query points, or just a reference set which will be used as
both the reference and query set.

For example, the following will return 5 neighbors from the data for each point
in `input` and store the distances in `distances` and the neighbors in
`neighbors`:

```julia
julia> using CSV
julia> input = CSV.read("input.csv")
julia> distances, neighbors, _ = lsh(k=5, reference=input)
```

The output is organized such that row i and column j in the neighbors output
corresponds to the index of the point in the reference set which is the j'th
nearest neighbor from the point in the query set with index i.  Row j and column
i in the distances output file corresponds to the distance between those two
points.

Because this is approximate-nearest-neighbors search, results may be different
from run to run.  Thus, the `seed` parameter can be specified to set the random
seed.

This program also has many other parameters to control its functionality; see
the parameter-specific documentation for more information.

# Arguments

 - `bucket_size::Int`: The size of a bucket in the second level hash.
      Default value `500`.

 - `hash_width::Float64`: The hash width for the first-level hashing in
      the LSH preprocessing. By default, the LSH class automatically estimates a
      hash width for its use.  Default value `0`.

 - `input_model::LSHSearch`: Input LSH model.
 - `k::Int`: Number of nearest neighbors to find.  Default value `0`.

 - `num_probes::Int`: Number of additional probes for multiprobe LSH; if
      0, traditional LSH is used.  Default value `0`.

 - `projections::Int`: The number of hash functions for each table
      Default value `10`.

 - `query::Array{Float64, 2}`: Matrix containing query points (optional).
 - `reference::Array{Float64, 2}`: Matrix containing the reference
      dataset.
 - `second_hash_size::Int`: The size of the second level hash table.
      Default value `99901`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `tables::Int`: The number of hash tables to be used.  Default value
      `30`.

 - `true_neighbors::Array{Int, 2}`: Matrix of true neighbors to compute
      recall with (the recall is printed when -v is specified).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `distances::Array{Float64, 2}`: Matrix to output distances into.
 - `neighbors::Array{Int, 2}`: Matrix to output neighbors into.
 - `output_model::LSHSearch`: Output for trained LSH model.

"""
function lsh(;
             bucket_size::Union{Int, Missing} = missing,
             hash_width::Union{Float64, Missing} = missing,
             input_model::Union{LSHSearch, Missing} = missing,
             k::Union{Int, Missing} = missing,
             num_probes::Union{Int, Missing} = missing,
             projections::Union{Int, Missing} = missing,
             query = missing,
             reference = missing,
             second_hash_size::Union{Int, Missing} = missing,
             seed::Union{Int, Missing} = missing,
             tables::Union{Int, Missing} = missing,
             true_neighbors = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, lshLibrary), Nothing, ());

  IORestoreSettings("K-Approximate-Nearest-Neighbor Search with LSH")

  # Process each input argument before calling mlpackMain().
  if !ismissing(bucket_size)
    IOSetParam("bucket_size", convert(Int, bucket_size))
  end
  if !ismissing(hash_width)
    IOSetParam("hash_width", convert(Float64, hash_width))
  end
  if !ismissing(input_model)
    lsh_internal.IOSetParamLSHSearch("input_model", convert(LSHSearch, input_model))
  end
  if !ismissing(k)
    IOSetParam("k", convert(Int, k))
  end
  if !ismissing(num_probes)
    IOSetParam("num_probes", convert(Int, num_probes))
  end
  if !ismissing(projections)
    IOSetParam("projections", convert(Int, projections))
  end
  if !ismissing(query)
    IOSetParamMat("query", query, points_are_rows)
  end
  if !ismissing(reference)
    IOSetParamMat("reference", reference, points_are_rows)
  end
  if !ismissing(second_hash_size)
    IOSetParam("second_hash_size", convert(Int, second_hash_size))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(tables)
    IOSetParam("tables", convert(Int, tables))
  end
  if !ismissing(true_neighbors)
    IOSetParamUMat("true_neighbors", true_neighbors, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("distances")
  IOSetPassed("neighbors")
  IOSetPassed("output_model")
  # Call the program.
  lsh_mlpackMain()

  return IOGetParamMat("distances", points_are_rows),
         IOGetParamUMat("neighbors", points_are_rows),
         lsh_internal.IOGetParamLSHSearch("output_model")
end


text:
using Colors, TestImages, ImageIO, ImageShow, FileIO, ImageCore

# see https://juliaimages.org

"""
size(nums) = (3, 256, 256)
typeof(nums) = Array{UInt8,3}
"""

peppers = FileIO.load("src/Julia/this-and-that/peppers_color_256.png")

#img = testimage("lena_color_256")
#img = testimage("lena_color_512")
img = peppers

# display(img)

p55=img[5,5]
@show p55
@show typeof(p55)
@show p55.r p55.g p55.b

p55_g = reinterpret(UInt8, p55.g) # 0x85 = 133

cn = channelview(img)
r = reinterpret(UInt8, cn) #
r2 = rawview(cn)
@show r == r2
r3 = normedview(r)
i3 = colorview(RGB, r3) # image again

ihsl = convert.(HSL, img) # hsl image
i_same = convert.(RGB{N0f8}, img)

@show r[:,5,5]

rc = collect(r)./255.0;
@show typeof(rc)
@show size(rc)

# normedview
# rawview

i2 = colorview(RGB, rc)

function isim(img::I) where I <: AbstractArray{C,2} where C <: Color # RGB
    isN0f8 = (eltype(img) <: RGB{N0f8})
    return (;isN0f8)
end
function isim(img)
    false
end
r1 = [1 2 3 4; 3 4 5 6]
g=r1.*2
b=r1.*3
rgb = zeros(3,2,4)
rgb[1,:,:] .= r1
rgb[2,:,:] .= g
rgb[3,:,:] .= b
rgb = UInt8.(rgb)
nrgb=normedview(rgb)
cl=colorview(RGB, nrgb)
@show isim(cl)


text:
"""
    shakhov(u, M, q, prim, Pr)
    shakhov(u, H, B, q, prim, Pr, K)
    shakhov(u, v, M, q, prim, Pr)
    shakhov(u, v, H, B, q, prim, Pr, K)
    shakhov(u, v, w, M, q, prim, Pr)

Shakhov non-equilibrium part

* @arg: particle velocity quadrature points
* @arg: discrete Maxwellian
* @arg: primitive variables, Prandtl number, heat flux, inner degree of freedom

"""
function shakhov(
    u::AbstractVector{X},
    M::AbstractVector{Y},
    q,
    prim::AbstractVector{Z},
    Pr,
) where {X<:AbstractFloat,Y<:AbstractFloat,Z<:Real} # 1F1V

    M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 - 5.0) *
       M

    return M_plus

end

#--- 2F1V ---#
function shakhov(
    u::AbstractVector{X},
    H::Y,
    B::Y,
    q,
    prim::AbstractVector{Z},
    Pr,
    K,
) where {
    X<:AbstractFloat,
    Y<:AbstractVector{<:AbstractFloat},
    Z<:Real,
}

    H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 + K - 5.0) *
       H
    B_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 + K - 3.0) *
       B

    return H_plus, B_plus

end

#--- 1F2V ---#
function shakhov(
    u::T,
    v::T,
    M::AbstractMatrix{X},
    q::AbstractVector{Y},
    prim::AbstractVector{Z},
    Pr,
) where {
    T<:AbstractArray{<:AbstractFloat,2},
    X<:AbstractFloat,
    Y<:AbstractFloat,
    Z<:Real,
}

    M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) - 5.0) *
       M

    return M_plus

end

#--- 2F2V ---#
function shakhov(
    u::T,
    v::T,
    H::X,
    B::X,
    q::AbstractVector{Y},
    prim::AbstractVector{Z},
    Pr,
    K,
) where {
    T<:AbstractArray{<:AbstractFloat,2},
    X<:AbstractArray{<:AbstractFloat,2},
    Y<:Real,
    Z<:Real,
}

    H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 5.0) *
       H
    B_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 3.0) *
       B

    return H_plus, B_plus

end

#--- 1F3V ---#
function shakhov(
    u::T,
    v::T,
    w::T,
    M::AbstractArray{X,3},
    q::AbstractVector{Y},
    prim::AbstractVector{Z},
    Pr,
) where {
    T<:AbstractArray{<:AbstractFloat,3},
    X<:AbstractFloat,
    Y<:Real,
    Z<:Real,
}

    M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       ((u - prim[2]) * q[1] + (v - prim[3]) * q[2] + (w - prim[4]) * q[3]) *
       (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) - 5.0) *
       M

    return M_plus

end

"""
    shakhov!(S, u, M, q, prim, Pr)
    shakhov!(SH, SB, u, H, B, q, prim, Pr, K)
    shakhov!(S, u, v, M, q, prim, Pr)
    shakhov!(SH, SB, u, v, H, B, q, prim, Pr, K)
    shakhov!(S, u, v, w, M, q, prim, Pr)

In-place Shakhov non-equilibrium part

* @arg: particle velocity quadrature points
* @arg: discrete Maxwellian
* @arg: primitive variables, Prandtl number, heat flux, inner degree of freedom

"""
function shakhov!(
    S::T1,
    u::AbstractVector{T2},
    M::T1,
    q,
    prim,
    Pr,
) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractFloat} # 1F1V

    @. S = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
       (u - prim[2]) *
       q *
       (2.0 * prim[end] * (u - prim[2])^2 - 5.0) *
       M

    return nothing

end

#--- 2F1V ---#
function shakhov!(
    SH::T1,
    SB::T1,
    u::AbstractVector{T2},
    H::T1,
    B::T1,
    q,
    prim,
    Pr,
    K,
) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractFloat}

    @. SH =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        (u - prim[2]) *
        q *
        (2.0 * prim[end] * (u - prim[2])^2 + K - 5.0) *
        H
    @. SB =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        (u - prim[2]) *
        q *
        (2.0 * prim[end] * (u - prim[2])^2 + K - 3.0) *
        B

    return nothing

end

#--- 1F2V ---#
function shakhov!(
    S::T1,
    u::T2,
    v::T2,
    M::T1,
    q::AbstractVector{T3},
    prim,
    Pr,
) where {
    T1<:AbstractArray{<:AbstractFloat,2},
    T2<:AbstractArray{<:AbstractFloat,2},
    T3<:Real,
}

    @. S =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) - 5.0) *
        M

    return nothing

end

#--- 2F2V ---#
function shakhov!(
    SH::T1,
    SB::T1,
    u::T2,
    v::T2,
    H::T1,
    B::T1,
    q::AbstractVector{T3},
    prim,
    Pr,
    K,
) where {
    T1<:AbstractArray{<:AbstractFloat,2},
    T2<:AbstractArray{<:AbstractFloat,2},
    T3<:Real,
}

    @. SH =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 5.0) *
        H
    @. SB =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2) + K - 3.0) *
        B

    return nothing

end

#--- 1F3V ---#
function shakhov!(
    S::T1,
    u::T2,
    v::T2,
    w::T2,
    M::T1,
    q::AbstractVector{T3},
    prim,
    Pr,
) where {
    T1<:AbstractArray{<:AbstractFloat,3},
    T2<:AbstractArray{<:AbstractFloat,3},
    T3<:Real,
}

    @. S =
        0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *
        ((u - prim[2]) * q[1] + (v - prim[3]) * q[2] + (w - prim[4]) * q[3]) *
        (2.0 * prim[end] * ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) - 5.0) *
        M

    return nothing

end


text:
struct StanModel{R,C,N}
  return_code::R
  chains::C
  cnames::N
end

struct StanODEData
end

function generate_differential_equation(f)
  theta_ex = MacroTools.postwalk(f.fex) do x
    if typeof(x) <: Expr && x.args[1] == :internal_var___p
      return Symbol("theta[$(x.args[2])]")
    else
      return x
    end
  end
  differential_equation = ""
  for i in 1:length(theta_ex.args)-1
    differential_equation = string(differential_equation,theta_ex.args[i], ";\n")
  end
  return differential_equation
end

function generate_priors(f,priors)
  priors_string = ""
  params = f.params
  if priors==nothing
    for i in 1:length(params)
      priors_string = string(priors_string,"theta[$i] ~ normal(0, 1)", " ; ")
    end
  else
    for i in 1:length(params)
      priors_string = string(priors_string,"theta[$i] ~",stan_string(priors[i]),";")
    end
  end
  priors_string
end

function generate_theta(n,priors)
  theta = ""
  for i in 1:n
    upper_bound = ""
    lower_bound = ""
    if maximum(priors[i]) != Inf
      upper_bound = string("upper=",maximum(priors[i]))
    end
    if minimum(priors[i]) != -Inf
      lower_bound = string("lower=",minimum(priors[i]))
    end
    if lower_bound != "" && upper_bound != ""
      theta = string(theta,"real","<$lower_bound",",","$upper_bound>"," theta$i",";")
    elseif lower_bound != ""
      theta = string(theta,"real","<$lower_bound",">"," theta$i",";")
    elseif upper_bound != ""
      theta = string(theta,"real","<","$upper_bound>"," theta$i",";")
    else
      theta = string(theta,"real"," theta$i",";")
    end
  end
  return theta
end

function stan_inference(prob::DiffEqBase.DEProblem,t,data,priors = nothing;alg=:rk45,
                            num_samples=1000, num_warmup=1000, reltol=1e-3,
                            abstol=1e-6, maxiter=Int(1e5),likelihood=Normal,
                            vars=(StanODEData(),InverseGamma(3,3)),nchains=1)
  length_of_y = length(prob.u0)
  length_of_params = length(vars)
  f = prob.f
  length_of_parameter = string(length(f.params))
  if alg ==:rk45
    algorithm = "integrate_ode_rk45"
  elseif alg == :bdf
    algorithm = "integrate_ode_bdf"
  else
    error("The choices for alg are :rk45 or :bdf")
  end
  hyper_params = ""
  tuple_hyper_params = ""
  setup_params = ""
  thetas = ""
  theta_string = generate_theta(length(priors),priors)
  for i in 1:length(priors)
    thetas = string(thetas,"theta[$i] <- theta$i",";")
  end
  for i in 1:length_of_params
    if isa(vars[i],StanODEData)
      tuple_hyper_params = string(tuple_hyper_params,"u_hat[t,:]",",")
    else
      dist = stan_string(vars[i])
      hyper_params = string(hyper_params,"sigma$(i-1) ~ $dist;")
      tuple_hyper_params = string(tuple_hyper_params,"sigma$(i-1)",",")
      setup_params = string(setup_params,"row_vector<lower=0>[$length_of_y] sigma$(i-1);")
    end
  end
  tuple_hyper_params = tuple_hyper_params[1:length(tuple_hyper_params)-1]
  differential_equation = generate_differential_equation(f)
  priors_string = string(generate_priors(f,priors))
  stan_likelihood = stan_string(likelihood)
  parameter_estimation_model = "
  functions {
    real[] sho(real t,real[] internal_var___u,real[] theta,real[] x_r,int[] x_i) {
      real internal_var___du[$length_of_y];
      $differential_equation
      return internal_var___du;
      }
    }
  data {
    real u0[$length_of_y];
    int<lower=1> T;
    real internal_var___u[T,$length_of_y];
    real t0;
    real ts[T];
  }
  transformed data {
    real x_r[0];
    int x_i[0];
  }
  parameters {
    $setup_params
    $theta_string
  }
  transformed parameters{
    real theta[$length_of_parameter];
    $thetas
  }
  model{
    real u_hat[T,$length_of_y];
    $hyper_params
    $priors_string
    u_hat = $algorithm(sho, u0, t0, ts, theta, x_r, x_i, $reltol, $abstol, $maxiter);
    for (t in 1:T){
      internal_var___u[t,:] ~ $stan_likelihood($tuple_hyper_params);
      }
  }
  "
  stanmodel = CmdStan.Stanmodel(num_samples=num_samples, num_warmup=num_warmup, name="parameter_estimation_model", model=parameter_estimation_model, nchains=nchains);
  parameter_estimation_data = Dict("u0"=>prob.u0, "T" => length(t), "internal_var___u" => data', "t0" => prob.tspan[1], "ts" => t)
  return_code, chains, cnames = CmdStan.stan(stanmodel, [parameter_estimation_data]; CmdStanDir=CMDSTAN_HOME)
  return StanModel(return_code, chains, cnames)
end

get_plot_object(bayesian_result::StanModel) = Mamba.plot(bayesian_result.chain_results)

function plot_chain(bayesian_result::StanModel,filename=nothing)
	p = get_plot_object(bayesian_result)
	if filename === nothing
		return Mamba.draw(p)
	else
		return Mamba.draw(p, filename=filename)
	end
end


text:
## regular expressions to use
whiteRe = r"\s*"
spaceRe = r"\s+"
nonSpaceRe = r"\S"
eqRe = r"\s*="
curlyRe = r"\s*\}"
#tagRe = r"#|\^|\/|>|\{|&|=|!"
tagRe = r"^[#^/<>{&=!]"

isWhitespace(x) = occursin(whiteRe, x)
function stripWhitepace(x)
    y = replace(x, r"^\s+" => "")
    replace(y, r"\s+$" => "")
end

## this is for falsy value
## Falsy is true if x is false, 0 length, "", ...
falsy(x::Bool) = !x
falsy(x::Array) = isempty(x)
falsy(x::AbstractString) = x == ""
falsy(x::Nothing) = true
falsy(x::Real) = x == 0
falsy(x) = (x == nothing) || false                #  default

## escape_html with entities

entityMap = [("&", "&"),
             ("<", "<"),
             (">", ">"),
             ("'", "'"),
             ("\"", """),
             ("/", "/")]

function escape_html(x)
    y = string(x)
    for (k,v) in entityMap
        y = replace(y, k => v)
    end
    y
end

## Make these work
function escapeRe(string)
    replace(string, r"[\-\[\]{}()*+?.,\\\^$|#\s]" => "\\\$&");
end

function escapeTags(tags)
   [Regex(escapeRe(tags[1]) * "\\s*"),
    Regex("\\s*" * escapeRe(tags[2]))]
end

## hueristic to avoid loading DataFrames
is_dataframe(x) = !isa(x, Dict) && !isa(x, Module) &&!isa(x, Array) && occursin(r"DataFrame", string(typeof(x)))


text:
# Spatial calculations

export
    vector1D,
    vector,
    wrap_coords,
    wrap_coords_vec,
    maxwell_boltzmann,
    random_velocities!

"""
    vector1D(c1, c2, side_length)

Displacement between two 1D coordinate values from c1 to c2, accounting for
the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinate accounting for the periodic boundaries.
"""
function vector1D(c1, c2, side_length)
    if c1 < c2
        return (c2 - c1) < (c1 - c2 + side_length) ? (c2 - c1) : (c2 - c1 - side_length)
    else
        return (c1 - c2) < (c2 - c1 + side_length) ? (c2 - c1) : (c2 - c1 + side_length)
    end
end

"""
    vector(c1, c2, box_size)

Displacement between two coordinate values from c1 to c2, accounting for
the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinates accounting for the periodic boundaries.
"""
vector(c1, c2, box_size) = vector1D.(c1, c2, box_size)

@generated function vector(c1::SVector{N}, c2::SVector{N}, box_size) where N
    quote
        Base.Cartesian.@ncall $N SVector{$N} i->vector1D(c1[i], c2[i], box_size[i])
    end
end

square_distance(i, j, coords, box_size) = sum(abs2, vector(coords[i], coords[j], box_size))

# Pad a vector to 3D to allow operations such as the cross product
function vector_pad3D(c1::SVector{2, T}, c2::SVector{2, T}, box_size::SVector{2, T}) where T
    SVector{3, T}(
        vector1D(c1[1], c2[1], box_size[1]),
        vector1D(c1[2], c2[2], box_size[2]),
        zero(T),
    )
end

vector_pad3D(c1::SVector{3}, c2::SVector{3}, box_size::SVector{3}) = vector(c1, c2, box_size)

# Trim a vector back to 2D if required
trim3D(v::SVector{3, T}, box_size::SVector{2}) where T = SVector{2, T}(v[1], v[2])
trim3D(v::SVector{3}, box_size::SVector{3}) = v

"""
    wrap_coords(c, side_length)

Ensure a 1D coordinate is within the simulation box and return the coordinate.
"""
wrap_coords(c, side_length) = c - floor(c / side_length) * side_length

"""
    wrap_coords_vec(c, box_size)

Ensure a coordinate is within the simulation box and return the coordinate.
"""
wrap_coords_vec(v, box_size) = wrap_coords.(v, box_size)

"""
    velocity(mass, temperature; dims=3)

Generate a random velocity from the Maxwell-Boltzmann distribution.
"""
function AtomsBase.velocity(mass, temp; dims::Integer=3)
    return SVector([maxwell_boltzmann(mass, temp) for i in 1:dims]...)
end

"""
    maxwell_boltzmann(mass, temperature)

Generate a random speed along one dimension from the Maxwell-Boltzmann distribution.
"""
function maxwell_boltzmann(mass, temp)
    T = typeof(convert(AbstractFloat, ustrip(temp)))
    k = unit(temp) == NoUnits ? one(T) : uconvert(u"u * nm^2 * ps^-2 * K^-1", T(Unitful.k))
    σ = sqrt(k * temp / mass)
    return rand(Normal(zero(T), T(ustrip(σ)))) * unit(σ)
end

"""
    random_velocities!(sys, temp)

Set the velocities of a `System` to random velocities generated from the
Maxwell-Boltzmann distribution.
"""
function random_velocities!(sys::System, temp)
    sys.velocities = [velocity(a.mass, temp) for a in sys.atoms]
    return sys
end


text:
using Color

# Test the colordiff function against example input and output from:
#
#   Sharma, G., Wu, W., & Dalal, E. N. (2005). The CIEDE2000 color‐difference
#   formula: Implementation notes, supplementary test data, and mathematical
#   observations. Color Research & Application, 30(1), 21–30. doi:10.1002/col
#

const abds = [
    ((50.0000,   2.6772, -79.7751), (50.0000,   0.0000, -82.7485),  2.0425),
    ((50.0000,   3.1571, -77.2803), (50.0000,   0.0000, -82.7485),  2.8615),
    ((50.0000,   2.8361, -74.0200), (50.0000,   0.0000, -82.7485),  3.4412),
    ((50.0000,  -1.3802, -84.2814), (50.0000,   0.0000, -82.7485),  1.0000),
    ((50.0000,  -1.1848, -84.8006), (50.0000,   0.0000, -82.7485),  1.0000),
    ((50.0000,  -0.9009, -85.5211), (50.0000,   0.0000, -82.7485),  1.0000),
    ((50.0000,   0.0000,   0.0000), (50.0000,  -1.0000,   2.0000),  2.3669),
    ((50.0000,  -1.0000,   2.0000), (50.0000,   0.0000,   0.0000),  2.3669),
    ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0009),  7.1792),
    ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0010),  7.1792),
    ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0011),  7.2195),
    ((50.0000,   2.4900,  -0.0010), (50.0000,  -2.4900,   0.0012),  7.2195),
    ((50.0000,  -0.0010,   2.4900), (50.0000,   0.0009,  -2.4900),  4.8045),
    ((50.0000,  -0.0010,   2.4900), (50.0000,   0.0010,  -2.4900),  4.8045),
    ((50.0000,  -0.0010,   2.4900), (50.0000,   0.0011,  -2.4900),  4.7461),
    ((50.0000,   2.5000,   0.0000), (50.0000,   0.0000,  -2.5000),  4.3065),
    ((50.0000,   2.5000,   0.0000), (73.0000,  25.0000, -18.0000), 27.1492),
    ((50.0000,   2.5000,   0.0000), (61.0000,  -5.0000,  29.0000), 22.8977),
    ((50.0000,   2.5000,   0.0000), (56.0000, -27.0000,  -3.0000), 31.9030),
    ((50.0000,   2.5000,   0.0000), (58.0000,  24.0000,  15.0000), 19.4535),
    ((50.0000,   2.5000,   0.0000), (50.0000,   3.1736,   0.5854),  1.0000),
    ((50.0000,   2.5000,   0.0000), (50.0000,   3.2972,   0.0000),  1.0000),
    ((50.0000,   2.5000,   0.0000), (50.0000,   1.8634,   0.5757),  1.0000),
    ((50.0000,   2.5000,   0.0000), (50.0000,   3.2592,   0.3350),  1.0000),
    ((60.2574, -34.0099,  36.2577), (60.4626, -34.1751,  39.4387),  1.2644),
    ((63.0109, -31.0961,  -5.8663), (62.8187, -29.7946,  -4.0864),  1.2530),
    ((61.2901,   3.7196,  -5.3901), (61.4292,   2.2480,  -4.9620),  1.8731),
    ((35.0831, -44.1164,   3.7933), (35.0232, -40.0716,   1.5901),  1.8645),
    ((22.7233,  20.0904, -46.6940), (23.0331,  14.9730, -42.5619),  2.0373),
    ((36.4612,  47.8580,  18.3852), (36.2715,  50.5065,  21.2231),  1.4146),
    ((90.8027,  -2.0831,   1.4410), (91.1528,  -1.6435,   0.0447),  1.4441),
    ((90.9257,  -0.5406,  -0.9208), (88.6381,  -0.8985,  -0.7239),  1.5381),
    (( 6.7747,  -0.2908,  -2.4247), ( 5.8714,  -0.0985,  -2.2286),  0.6377),
    ((2.0776,    0.0795,  -1.1350), ( 0.9033,  -0.0636,  -0.5514),  0.9082)]

const eps_cdiff = 0.01

metric = DE_2000()

for (i, (a, b, dexpect)) in enumerate(abds)
    @assert abs(dexpect - colordiff(Lab(a...), Lab(b...), metric)) < eps_cdiff
    @assert abs(dexpect - colordiff(Lab(b...), Lab(a...), metric)) < eps_cdiff
end


text:
using Printf
using Test

# Add tests below if you wish that they are not run together with all tests
const tests_not_to_run = Set{String}(map(uppercase, [
    "runtests.jl", # default test file
    "run_all_tests.jl", # this file
    "slow_test.jl", # too slow to run regularly
]))

function run_all_tests(test_name_regex = "")
    @testset "All tests" begin
        root = string(@__DIR__)
        tests_to_run = [joinpath(dir, file) for (dir, _, files) in walkdir(root) for file in files
            if occursin(Regex(test_name_regex), file) && is_test_script(dir, file) && !in(uppercase(file), tests_not_to_run)]

        for i = 1:length(tests_to_run)
            file = tests_to_run[i]
            test_name = replace(file, Regex("$root/?(.+).jl\$", "i") => s"\1")
            @printf("Running test %s (%d / %d)\n", test_name, i, length(tests_to_run))
            @time include(file)
        end
    end
end

function is_test_script(dir::AbstractString, file::AbstractString)
    # match files ending in ".jl" and not containing # or ~ (emacs temp files)
    if occursin(r"(?i)^[^#~]+\.jl$", file)
        src = read(joinpath(dir, file), String)

        pos = 1
        while pos <= length(src)
            expr, pos = Meta.parse(src, pos)
            expr === nothing && continue
            if expr.head == :incomplete
                msg = join(map(string, expr.args), "; ")
                throw(Meta.ParseError("While parsing file $file got invalid expression: $msg"))
            end
            contains_test_macro(expr) && return true
        end
    end
    return false
end

function contains_test_macro(expr::Expr)
    if expr.head == :macrocall && in(expr.args[1], [Symbol("@test") Symbol("@testset")])
        return true
    end
    return any(e -> contains_test_macro(e), filter(a -> isa(a, Expr), expr.args))
end

max_stars_repo_name: maxbennedich/julia-regression-analysis

max_issues_repo_name: maxbennedich/julia-regression-analysis

max_forks_repo_name: maxbennedich/julia-regression-analysis


text:
using ITensors

function insert(ψ::MPS, b::Integer, ψb::ITensor)
  dataψ′ = insert!(copy(ITensors.data(ψ)), b, ψb)
  ψ′ = MPS(dataψ′)
  lb = linkind(ψ, b - 1)
  l′b = sim(lb)
  δlb = δ(dag(lb), l′b)
  ψ′[b - 1] *= δlb
  ψ′[b] *= dag(δlb)
  return ψ′
end

N = 4
s = siteinds("S=1/2", N; conserve_qns=true)
ψ = randomMPS(s, n -> isodd(n) ? "↑" : "↓")

b = 2
sb = siteind("S=1/2"; conserve_qns=true)
ψb = onehot(sb => 1)
ψ′ = insert(ψ, b, ψb)

@show norm(prod(ψ) * ψb - prod(ψ′))


text:
struct SamplingLimit
    min::UInt16
    center::UInt16
    max::UInt16
end

struct SamplingDescription
    encoded_fov::NTuple{3,Float32}
    recon_fov::NTuple{3,Float32}
    encoded_matrix::NTuple{3,UInt16}
    recon_matrix::NTuple{3,UInt16}
    sampling_limits::NTuple{3,SamplingLimit}
end

struct ReconBuffer
    data::Array{ComplexF32}
    trajectory::Optional{Array{Float32}}
    density::Optional{Array{Float32}}
    header:: Array{AcquisitionHeader}
    sampling_description::SamplingDescription

end

struct ReconBit
    data::ReconBuffer
    ref::Optional{ReconBuffer}
end

struct ReconData
    bits::Vector{ReconBit}
end

Base.getindex(rd::ReconData,i::Int) = rd.bits[i]
Base.setindex(rd::ReconData,val::ReconBit,i::Int) = rd.bits[i] = val
Base.firstindex(rd::ReconData) = firstindex(rd.bits)
Base.lastindex(rd::ReconData) = lastindex(rd.bits)

MRD.read(io::IO, ::Type{ReconData}) = MRD.read_struct(io, ReconData)
MRD.read(io::IO, ::Type{ReconBit}) = MRD.read_struct(io, ReconBit)
MRD.read(io::IO, ::Type{ReconBuffer}) = MRD.read_struct(io, ReconBuffer)

function MRD.read(io::IO, ::Type{SamplingDescription})
    result = MRD.read_struct(io,SamplingDescription)
    MRD.read(io,UInt16) #Padding bytes
    return result
end

MRD.read(io::IO, ::Type{SamplingLimit}) = MRD.read_struct(io,SamplingLimit)


text:
function binary_op(op, W, @nospecialize(T))
  ty = LLVM_TYPES[T]
  if isone(W)
    V = T
  else
    ty = "<$W x $ty>"
    V = NTuple{W,VecElement{T}}
  end
  instrs = "%res = $op $ty %0, %1\nret $ty %res"
  call = :($LLVMCALL($instrs, $V, Tuple{$V,$V}, data(v1), data(v2)))
  W > 1 && (call = Expr(:call, :Vec, call))
  Expr(:block, Expr(:meta, :inline), call)
end

# Integer
for (op,f) ∈ [("add",:+),("sub",:-),("mul",:*),("shl",:<<)]
  ff = Symbol('v', op)
  fnsw = Symbol(ff,"_nsw")
  fnuw = Symbol(ff,"_nuw")
  fnw = Symbol(ff,"_nw")
  ff_fast = if op == "shl"
    @eval @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W, T)
  else
    ff_fast = :(FastMath.$(Symbol(op * "_fast")))
    @eval begin
      @generated $ff_fast(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W, T)
      @inline Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = $ff_fast(v1, v2)
    end
  end
  @eval begin
    @generated $fnsw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($(op * " nsw"), W, T)
    @generated $fnuw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($(op * " nuw"), W, T)
    @generated $fnw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($(op * " nsw nuw"), W, T)
    @generated $fnsw(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($(op * " nsw"), 1, T)
    @generated $fnuw(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($(op * " nuw"), 1, T)
    @generated $fnw(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($(op * " nsw nuw"), 1, T)
  end
end
for (op,f) ∈ [("div",:÷),("rem",:%)]
  @eval begin
    @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, W, T)
  end
end
@inline div_fast(x::T, y::T) where {T <: SignedHW} = Base.sdiv_int(x,y)
@inline div_fast(x::T, y::T) where {T <: UnsignedHW} = Base.udiv_int(x,y)
@inline div_fast(x::AbstractSIMD, y::AbstractSIMD) = x ÷ y
@inline vcld(x, y) = (div_fast((x - one(x)), y) + one(x))
@inline function divrem_fast(x, y)
  d = div_fast(x, y)
  r = x - (d * y)
  d, r
end
for (op,sub) ∈ [
  ("ashr",:SignedHW),
  ("lshr",:UnsignedHW),
  ("lshr",:IntegerTypesHW),
  ("and",:IntegerTypesHW),
  ("or",:IntegerTypesHW),
  ("xor",:IntegerTypesHW)
  ]
  ff = sub === :UnsignedHW ? :vashr : Symbol('v', op)
  @eval begin
    @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:$sub}  = binary_op($op, W, T)
    @generated $ff(v1::T, v2::T) where {T<:$sub}  = binary_op($op, 1, T)
  end
end

for (op,f,s) ∈ [("ashr",:>>,0x01),("lshr",:>>,0x02),("lshr",:>>>,0x03),("and",:&,0x03),("or",:|,0x03),("xor",:⊻,0x03)]
  fdef = Expr(:where, :(Base.$f(v1::Vec{W,T}, v2::Vec{W,T})), :W)
  if iszero(s & 0x01)
    @eval @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:UnsignedHW} = binary_op($op, W1, T, W2)
  elseif iszero(s & 0x02)
    @eval @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:SignedHW} = binary_op($op, W1, T, W2)
  else
    @eval @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W1, T, W2)
  end
  if op === "ashr" # replace first iteration with <<
    f = :(<<)
  end
  @eval begin
    @inline Base.$f(v::Vec, i::Integer) = ((x, y) = promote(v, i); $f(x, y))
    @inline Base.$f(i::Integer, v::Vec) = ((x, y) = promote(i, v); $f(x, y))
  end
end

for (op,f,ff) ∈ [("fadd",:(+),:add_fast),("fsub",:(-),:sub_fast),("fmul",:(*),:mul_fast),("fdiv",:(/),:div_fast),("frem",:(%),:rem_fast)]

  fop_fast = f === :(/) ? "fdiv fast" : op * ' ' * fast_flags(true)
  fop_contract = op * ' ' * fast_flags(false)
  @eval begin
    @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} = binary_op($fop_contract, W, T)
    @generated FastMath.$ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} = binary_op($fop_fast, W, T)
  end
end

@inline Base.div(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:FloatingTypes} = floor(signed(Base.uinttype(T)), FastMath.div_fast(v1, v2))

@inline Base.:(/)(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = float(a) / float(b)
@inline FastMath.div_fast(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = FastMath.div_fast(float(a), float(b))


text:
module FileFixturesTests

using Fixtures
using FactCheck

const TEMPDIR = tempdir()
const TEMPDIR_LEN = length(TEMPDIR)

facts("FileFixtures tests") do
  context("temporary filename") do
    temp_filename() do filename
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact ispath(filename) --> false
      touch(filename)
      @fact isfile(filename) --> true
      global remembered_filename = filename
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary filename turned into directory") do
    temp_filename() do filename
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact ispath(filename) --> false
      mkdir(filename)
      @fact isdir(filename) --> true
      global remembered_filename = filename
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary filename never created") do
    temp_filename() do filename
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact ispath(filename) --> false
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary filename with extension") do
    temp_filename(extension="csv") do filename
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact filename[end-2:end] --> "csv"
      @fact ispath(filename) --> false
      touch(filename)
      @fact isfile(filename) --> true
      global remembered_filename = filename
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary filename with create") do
    temp_filename(create=true) do filename
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact isfile(filename) --> true
      open(filename) do f
        @fact readlines(f) --> []
      end
      global remembered_filename = filename
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary filename with create from content") do
    const mycontent = "Hello World"
    temp_filename(create=mycontent) do filename
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact isfile(filename) --> true
      open(filename) do f
        @fact readlines(f) --> [mycontent]
      end
      global remembered_filename = filename
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary file for output") do
    temp_file() do file
      write(file, "hello world")
      const filename = file.name[7:end-1]
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact isfile(filename) --> true
      global remembered_filename = filename
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary file for input") do
    const mycontent = "Hello World"
    temp_file(content=mycontent, mode="r") do file
      @fact readlines(file) --> [mycontent]
      const filename = file.name[7:end-1]
      @fact filename[1:TEMPDIR_LEN] --> TEMPDIR
      @fact isfile(filename) --> true
      global remembered_filename = filename
    end
    @fact ispath(remembered_filename) --> false
  end

  context("temporary directory") do
    temp_dir() do dirname
      @fact dirname[1:TEMPDIR_LEN] --> TEMPDIR
      @fact isdir(dirname) --> true
      global remembered_dirname = dirname
    end
    @fact ispath(remembered_dirname) --> false
  end

  context("cleanup_path file") do
    const filename = tempname()
    cleanup_path(filename) do
      @fact ispath(filename) --> false
      touch(filename)
      @fact isfile(filename) --> true
    end
    @fact ispath(filename) --> false
  end

  context("cleanup_path directory") do
    const filename = tempname()
    cleanup_path(filename) do
      @fact ispath(filename) --> false
      mkdir(filename)
      @fact isdir(filename) --> true
    end
    @fact ispath(filename) --> false
  end

  context("cleanup_path missing (ignored)") do
    const filename = tempname()
    cleanup_path(filename, ignore_missing=true) do
      @fact ispath(filename) --> false
    end
    @fact ispath(filename) --> false
  end

  context("cleanup_path missing (not ignored)") do
    const filename = tempname()
    @fact_throws cleanup_path(filename) do
      @fact ispath(filename) --> false
    end
  end
end

end


text:
@testset "Genus representatives" begin
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x-1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [-717210130, 0, 0, 0, -55, 0, 0, 0, -1173298395869600]);
  gens = [[1, 0, 0], [1, 0, 0], [4//5, 2//5, 0], [4//5, 2//5, 0], [1493459//5517001, 0, 1//15778622860], [1493459//5517001, 0, 1//15778622860]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 6

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x-1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [-34, 0, 0, 0, -17, 0, 0, 0, -17192032]);
  gens = [[26, 0, 0], [234, 0, 0], [177//17, 2//17, 0], [177//17, 2//17, 0], [27//2, 0, 1//1768], [27//2, 0, 1//1768]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 93

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x-1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [-98, 0, 0, 0, -2, 0, 0, 0, -5829824]);
  gens = [[8//7, 0, 0], [8//7, 0, 0], [5//14, 13//2, 0], [45//14, 117//2, 0], [29//28, 3//4, 1//1456], [29//28, 3//4, 1//1456]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 114

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x^2-2;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [2, 0, 0, 0, -1, 0, 0, 0, -7436]);
  gens = [[13, 0, 0], [156*a+143, 0, 0], [3//2*a+5, 1, 0], [3//2*a+5, 1, 0], [21//2*a, 0, 1//26], [21//2*a, 0, 1//26]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 1

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1
  K, a = number_field(f)
  M = matrix(QQ, 2, 2, [47, 80, 80, 560])
  w = Int[]
  for i in 1:10
    M[2, 2] = 560 + i
    V = quadratic_space(K, M)
    L = lattice(V, identity_matrix(K, 2))
    g = length(Hecke._genus_representatives_binary_quadratic(L))
    push!(w, g)
  end
  @test w == Int[37, 16, 19, 11, 6, 3, 52, 7, 13, 21]

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x^2 - 2;
  K, a = number_field(f)
  D = matrix(K, 2, 2, [15, 2, 2, 32]);
  gens = [[1, 0], [1, 0], [0, 1], [0, 1]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  @test length(genus_representatives(L)) == 6
end


text:
function pcd!(
    rbm::DenseConvRBM,
    data_dense::AbstractArray,
    data_conv::AbstractArray;
    batchsize::Int = 1,
    epochs::Int = 1,
    optim = Flux.ADAM(),
    history::MVHistory = MVHistory(),
    wts = nothing,
    steps::Int = 1,
    vm_dense::AbstractArray = initial_fantasy_v(parts(rbm).dense, data_dense, batchsize),
    vm_conv::AbstractArray = initial_fantasy_v(parts(rbm).conv, data_conv, batchsize)
)
    @assert size(data_dense) == (size(rbm.visible_dense)..., _nobs(data_dense))
    @assert only(vsizes(parts(rbm).conv, data_conv).batch_size) == _nobs(data_dense)
    @assert size(vm_dense) == (size(rbm.visible_dense)..., batchsize)
    @assert only(vsizes(parts(rbm).conv, vm_conv).batch_size) == batchsize
    @assert isnothing(wts) || _nobs(data_dense) == _nobs(wts)

    for epoch in 1:epochs
        batches = minibatches(data_dense, data_conv, wts; batchsize)
        Δt = @elapsed for (vd_dense, vd_conv, wd) in batches
            _vm_dense, _vm_conv = RBMs.sample_v_from_v(rbm, vm_dense, vm_conv; steps)
            vm_dense .= _vm_dense
            vm_conv  .= _vm_conv
            gs = Zygote.gradient(rbm) do rbm
                contrastive_divergence(rbm, vd_dense, vd_conv, vm_dense, vm_conv; wd)
            end
            # removes keys with no gradient
            ∂ = Base.structdiff(only(gs), NamedTuple{(:stride, :pad, :dilation, :groups, :pool)})
            push!(history, :∂, RBMs.gradnorms(∂))
            RBMs.update!(rbm, RBMs.update!(∂, rbm, optim))
            push!(history, :Δ, RBMs.gradnorms(∂))
        end
        push!(history, :epoch, epoch)
        push!(history, :Δt, Δt)
        push!(history, :vm, copy(vm))
        @debug "epoch $epoch/$epochs ($(round(Δt, digits=2))s)"
    end
    return history
end

function contrastive_divergence(
    rbm::DenseConvRBM,
    vd_dense::AbstractArray, vd_conv::AbstractArray,
    vm_dense::AbstractArray, vm_conv::AbstractArray;
    wd = nothing, wm = nothing
)
    Fd = mean_free_energy(rbm, vd_dense, vd_conv; wts = wd)
    Fm = mean_free_energy(rbm, vm_dense, vm_conv; wts = wm)
    return Fd - Fm
end

function mean_free_energy(
    rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray; wts = nothing
)::Number
    F = RBMs.free_energy(rbm, v_dense, v_conv)
    if isempty(batch_size(rbm, v_dense, v_conv))
        wts::Nothing
        return F
    else
        return RBMs.wmean(F; wts)
    end
end

function RBMs.update!(∂::NamedTuple, x::DenseConvRBM, optim)
    for (k, g) in pairs(∂)
        if hasproperty(x, k)
            RBMs.update!(g, getproperty(x, k), optim)
        else
            g .= 0
        end
    end
    return ∂
end

function RBMs.update!(x::DenseConvRBM, ∂::NamedTuple)
    for (k, Δ) in pairs(∂)
        hasproperty(x, k) && RBMs.update!(getproperty(x, k), Δ)
    end
    return x
end

function initial_fantasy_v(rbm::RBM, data::AbstractArray, batchsize::Int)
    inputs = falses(size(data)[1:(end - 1)]..., batchsize)
    vm = RBMs.transfer_sample(visible(rbm), inputs)
    return oftype(data, vm)
end


text:
# Line geometry connects (x, y) coordinates with lines.

struct LineGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    # Do not reorder points along the x-axis.
    preserve_order::Bool

    order::Int

    tag::Symbol
end

function LineGeometry(default_statistic=Gadfly.Stat.identity();
                      preserve_order=false, order=2, tag=empty_tag)
    LineGeometry(default_statistic, preserve_order, order, tag)
end

"""
    Geom.line[(; preserve_order=false, order=2)]

Draw a line connecting the `x` and `y` coordinates.  Optionally plot multiple
lines according to the `group` or `color` aesthetics.  `order` controls whether
the lines(s) are underneath or on top of other forms.

Set `preserve_order` to `:true` to *not* sort the points according to their
position along the x axis, or use the equivalent [`Geom.path`](@ref) alias.
"""
const line = LineGeometry

### why would one ever want to set preserve_order to false here
"""
    Geom.contours[(; levels=15, samples=150, preserve_order=true)]

Draw contour lines of the 2D function, matrix, or DataFrame in the `z`
aesthetic.  This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.contour`](@ref); see the latter for more information.
"""
function contour(; levels=15, samples=150, preserve_order=true)
    return LineGeometry(Gadfly.Stat.contour(levels=levels, samples=samples),
                                            preserve_order=preserve_order)
end

# Only allowing identity statistic in paths b/c I don't think any
# any of the others will work with `preserve_order=true` right now
"""
    Geom.path

Draw lines between `x` and `y` points in the order they are listed.  This
geometry is equivalent to [`Geom.line`](@ref) with `preserve_order=true`.
"""
path() = LineGeometry(preserve_order=true)

"""
    Geom.density[(; bandwidth=-Inf)]

Draw a line showing the density estimate of the `x` aesthetic.
This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.density`](@ref); see the latter for more information.
"""
density(; bandwidth::Real=-Inf) =
    LineGeometry(Gadfly.Stat.density(bandwidth=bandwidth))

"""
    Geom.density2d[(; bandwidth=(-Inf,-Inf), levels=15)]

Draw a set of contours showing the density estimate of the `x` and `y`
aesthetics.  This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.density2d`](@ref); see the latter for more information.
"""
density2d(; bandwidth::Tuple{Real,Real}=(-Inf,-Inf), levels=15) =
    LineGeometry(Gadfly.Stat.density2d(bandwidth=bandwidth, levels=levels); preserve_order=true)

"""
    Geom.smooth[(; method:loess, smoothing=0.75)]

Plot a smooth function estimated from the line described by `x` and `y`
aesthetics.  Optionally group by `color` and plot multiple independent smooth
lines.  This geometry is equivalent to [`Geom.line`](@ref) with
[`Stat.smooth`](@ref); see the latter for more information.
"""
smooth(; method::Symbol=:loess, smoothing::Float64=0.75) =
    LineGeometry(Gadfly.Stat.smooth(method=method, smoothing=smoothing),
    order=5)

"""
    Geom.step[(; direction=:hv)]

Connect points described by the `x` and `y` aesthetics using a stepwise
function.  Optionally group by `color` or `group`.  This geometry is equivalent
to [`Geom.line`](@ref) with [`Stat.step`](@ref); see the latter for more
information.
"""
step(; direction::Symbol=:hv) = LineGeometry(Gadfly.Stat.step(direction=direction))

default_statistic(geom::LineGeometry) = geom.default_statistic

element_aesthetics(::LineGeometry) = [:x, :y, :color, :group, :linestyle]

function Gadfly.Geom.render(geom::LineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.line", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.line", aes, Geom.element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.group = IndirectArray(fill(1,length(aes.x)))
    default_aes.color = fill(theme.default_color, length(aes.x))
    default_aes.linestyle = fill(1, length(aes.x))
    aes = Gadfly.inherit(aes, default_aes)

    # Point order:
    p = 1:length(aes.x)
    !geom.preserve_order && (p = sortperm(aes.x))
    aes_x, aes_y, aes_color, aes_g = aes.x[p], aes.y[p], aes.color[p], aes.group[p]
    aes_linestyle = aes.linestyle[p]

    # Find the aesthetic with the most levels:
    aesv = [aes_g, aes_color, aes_linestyle]
    i1 = argmax([length(unique(a)) for a in aesv])
    aes_maxlvls = aesv[i1]

    # Concrete values?:
    cf = Gadfly.isconcrete.(aes_x) .& Gadfly.isconcrete.(aes_y)
    fcf = .!cf
    ulvls =  unique(aes_maxlvls[fcf])
    aes_concrete = zeros(Int, length(cf))
    for g in ulvls
        i = aes_maxlvls.==g
        aes_concrete[i] = cumsum(fcf[i])
    end

    aes_x, aes_y, aes_color, aes_g = aes_x[cf], aes_y[cf], aes_color[cf], aes_g[cf]
    aes_concrete, aes_linestyle = aes_concrete[cf], aes_linestyle[cf]

    # Render lines, using multivariate groupings:
    XT, YT, CT, GT, CNT = eltype(aes_x), eltype(aes_y), eltype(aes_color), eltype(aes_g), eltype(aes_concrete)
    LST = eltype(aes_linestyle)
    groups = collect((Tuple{GT, CT, LST, CNT}), zip(aes_g, aes_color, aes_linestyle, aes_concrete))
    ug = unique(groups)

    n = length(ug)
    lines = Vector{Vector{Tuple{XT,YT}}}(undef, n)
    line_colors = Vector{CT}(undef, n)
    line_styles = Vector{LST}(undef, n)
    linestyle_palette_length = length(theme.line_style)
    for (k,g) in enumerate(ug)
        i = groups.==[g]
        lines[k] = collect(Tuple{XT,YT}, zip(aes_x[i], aes_y[i]))
        line_colors[k] = first(aes_color[i])
        line_styles[k] = mod1(first(aes_linestyle[i]), linestyle_palette_length)
    end

    linestyles =  Gadfly.get_stroke_vector.(theme.line_style[line_styles])
    classes = svg_color_class_from_label.(aes.color_label(line_colors))
    ctx = context(order=geom.order)
    ctx = compose!(ctx, (context(), Compose.line(lines, geom.tag),
                        stroke(line_colors),
                        strokedash(linestyles),
                        svgclass(classes)), svgclass("geometry"))

    return compose!(ctx, fill(nothing), linewidth(theme.line_width))

end


text:
using CSV
using Plots
using Unitful
using UnitfulRecipes
using Measurements

# We're going to overengineer a bit.

struct Calibration
    wavelength::Vector{typeof(1.0u"nm")}
    counterAvg::Vector{Measurement}
    gradient
    yintercept
end

"""Load calibration data from path string and construct Calibration object."""
function Calibration(path::String)
    headers = ["order", "wavelength", "energy", "min", "max", "avg", "err"] # Custom column headers
    calData = CSV.File(path, header=headers, skipto=2) # skip to first row of data
    wavelength = Quantity.(calData.:wavelength, u"nm")
    counterAvg = calData.:avg .± calData.:err

    # Let's use linear least squares to find the fit
    # Can fit to wavelength = m1 * counts + m2
    # This makes our data kernel:
    kernel = [counterAvg ones(length(counterAvg))]
    # Generalised inverse
    genInv = inv(transpose(kernel) * kernel) * transpose(kernel)
    parameters = genInv * wavelength
    #display(parameters)

    Calibration(wavelength, counterAvg, parameters...)
end

# Expects to be passed an object of type Calibration
@userplot CalibrationPlot
@recipe function f(i::CalibrationPlot)
    cal = i.args[1]
    rawX = cal.counterAvg
    rawY = cal.wavelength

    legend := :bottomright
    xguide := "Counts"
    yguide := "Wavelength"
    grid := true
    minorgrid := true
    seriescolor := :auto

    @series begin
        seriestype := :scatter
        label := "Measurements"
        rawX, rawY
    end

    @series begin
        seriestype := :line
        label := "Fit"
        linestyle := :dash
        y = cal.gradient * rawX .+ cal.yintercept
        rawX, y
    end

end

cal = Calibration("data_calibration_v2.csv")
display(calibrationplot(cal))

λ(count, cal::Calibration=cal) = cal.gradient * count + cal.yintercept
E(λ) = uconvert(u"eV", u"h"* u"c0" / λ)

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

data = CSV.File("data_sample_v2.csv")

mutable struct OpticalData
    counter::Vector{Int}
    wavelength::Vector{typeof(1.0u"nm")}
    energy::Vector{typeof((1.0 ± 0.1)u"eV")}
    RawSi::Vector{typeof(1.0u"mV")}
    RawInP::Vector{typeof(1.0u"mV")}
    ScaledSi::Vector{<:Real}
    ScaledInP::Vector{<:Real}
end

"""Load data from path string and construct OpticalData object."""
function OpticalData(path::String)
    headers = ["counter", "wavelength", "energy", "spacer", "withoutInP", "withoutInPScaled", "withInP", "spacer2", "withInPScaled"] # Custom column headers
    data = CSV.File(path, header=headers, skipto=2) # skip to first row of data

    counter = data.:counter
    wavelength = Quantity.(data.:wavelength, u"nm")
    energy = E.(λ.(counter))
    RawSi = Quantity.(data.:withoutInP, u"mV")
    RawInP = Quantity.(data.:withInP, u"mV")
    ScaledSi = RawSi ./ maximum(RawSi)
    ScaledInP = RawInP ./ maximum(RawInP)

    OpticalData(counter, wavelength, energy, RawSi, RawInP, ScaledSi, ScaledInP)
end

data = OpticalData("data_sample_v2.csv")

# Expects to be passed x=energy, y=relative intensity and an optional label or any other plot options
@userplot OpticalDataPlot
@recipe function f(i::OpticalDataPlot)
    x, y = i.args

    legend := :topright
    xguide := "Incident Energy"
    yguide := "Relative Intensity"
    grid := true
    minorgrid := true

    x, y
end

SiPlot = opticaldataplot(data.energy, data.ScaledSi, label="Si")
InPPlot = opticaldataplot(data.energy, data.ScaledInP, label="InP", xlims=(-Inf, 1.5))
display(SiPlot)
display(InPPlot)

@userplot DirectnessCheck
@recipe function f(i::DirectnessCheck)
    energy, intensity = i.args
    invLog = log.(data.RawSi ./ intensity)

    legend := :topleft
    xguide := "Incident Energy"
    yguide := "f(αx)"
    grid := true
    minorgrid := true

    @series begin
        label := "Direct"
        energy, invLog .^ 2
    end

    @series begin
        label := "Indirect"
        energy, sqrt.(invLog)
    end

    @series begin
        # Filter out infinity and NaN rows
        arr = [energy invLog]
        replace!(arr, Inf=>NaN)
        nanrows = any(isnan.(arr), dims=2)
        arr = arr[.!vec(nanrows), :]
        arr[:,2] = arr[:,2] .^ 2

        line = arr[arr[:,1] .> 1.34u"eV", :]
        kernel = [ustrip.(line[:,1]) ones(length(line[:,1]))]
        genInv = inv(transpose(kernel) * kernel) * transpose(kernel)
        parameters = genInv * line[:,2]
        display(parameters)

        fline(x) = parameters[1].val * x + parameters[2].val
        x = LinRange(1.35, 1.385, 2)
        y = fline.(x)

        linestyle := :dash
        label := "Fit"

        x, y
    end

end

LinearPlot = directnesscheck(data.energy, data.RawInP, xlims=(1.1, 1.4))
display(LinearPlot)

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


text:
"""
    module GCSPET

Package providing IO, visualization and basic modelling components for the "Gantry Crane Scheduling Problem with External Trucks", as presented in:

Guo, Peng, et al. "Gantry crane scheduling in intermodal rail-road container terminals." **International Journal of Production Research** 56.16 (2018): 5419-5436.
"""
module GCSPET

    export  Job,                # Job
            id,
            loc,
            t_arrival,
            t_earliest_start,
            t_processing,
            jobtype,
            movetype,
            istruck,
            istrain,
            isload,
            isunload,
            t_waiting,
            t_truckwaiting,
            Crane,              # Crane
            speed,
            starting_position,
            zone,
            safety,
            Instance,           # Instance
            name,
            jobs,
            cranes,
            crane_starting_pos,
            parameters,
            njobs,
            ncranes,
            loadfactor,
            instancepath,
            instancelist,
            instancedir,
            generate_instance,
            precedence_graph,
            precedence_matrix,
            assignment_matrix,
            assignment_ranges

    using Reexport
    using SplitApplyCombine: group # model
    @reexport using LightGraphs       # model
    using DataDeps          # io
    using Luxor             # visualization
    using Colors            # visualization
    using Test              # validation

    include("model/job.jl")
    include("model/crane.jl")
    include("model/instance.jl")
    include("model/solution.jl")
    include("model/bounds.jl")
    include("model/utilities.jl")

    include("io/utilities.jl")
    include("io/dependencies.jl")
    include("io/io.jl")
    include("io/generator.jl")

    include("visualization/base.jl")

    function __init__()
        register(guo_instances)
    end

end


text:
"""
    DockerRunner

A `DockerRunner` represents Elliot's bowing of his head to the inevitability
that the QemuRunner just isn't ready for primetime yet, and he needs something
to use while he's on the plane to JuliaCon to whip up said JuliaCon presentation.
"""
mutable struct DockerRunner <: Runner
    docker_cmd::Cmd
    prefix::AbstractString
    rootfs_version::VersionNumber
    platform::Platform
end

docker_image(version::VersionNumber) = "julia_binarybuilder_rootfs:v$(version)"
docker_image(rootfs::CompilerShard) = docker_image(rootfs.version)

"""
    import_docker_image(rootfs::CompilerShard; verbose::Bool = false)

Checks to see if the given rootfs has been imported into docker yet; if it
hasn't, then do so so that we can run things like:

    docker run -ti binarybuilder_rootfs:v2018.08.27 /bin/bash

Which, after all, is the foundation upon which this whole doodad is built.
"""
function import_docker_image(rootfs::CompilerShard; verbose::Bool = false)
    if rootfs.archive_type != :targz
        throw(ArgumentError("Unable to import squashfs into docker!  Use .tar.gz shards!"))
    end

    # Does this image already exist?  If so, we're done!
    if success(`docker inspect --type=image $(docker_image(rootfs))`)
        if verbose
            @info("Docker base image already exists, skipping import...")
        end
        return
    end

    # Otherwise, import it!
    dockerfile_cmds = "ENTRYPOINT [\"/docker_entrypoint.sh\"]"
    run(`docker import $(download_path(rootfs)) -c $(dockerfile_cmds) $(docker_image(rootfs))`)
    return
end

function DockerRunner(workspace_root::String;
                      cwd = nothing,
                      platform::Platform = platform_key_abi(),
                      workspaces::Vector = [],
                      extra_env=Dict{String, String}(),
                      verbose::Bool = false)
    global use_ccache

    # Check to make sure we're not going to try and bindmount within an
    # encrypted directory, as that can trigger kernel bugs
    check_encryption(workspace_root; verbose=verbose)

    # Choose and prepare our shards
    shards = choose_shards(platform)
    prepare_shard.(shards; mount_squashfs = false, verbose=verbose)

    # Import docker image
    import_docker_image(shards[1]; verbose=verbose)

    # Construct environment variables we'll use from here on out
    envs = merge(platform_envs(platform), extra_env)

    # the workspace_root is always a workspace, and we always mount it first
    insert!(workspaces, 1, workspace_root => "/workspace")

    # If we're enabling ccache, then map in a named docker volume for it
    if use_ccache
        push!(workspaces, "binarybuilder_ccache" => "/root/.ccache")
    end

    # Construct docker command
    docker_cmd = `docker run --privileged `#--cap-add SYS_ADMIN`

    if cwd != nothing
        docker_cmd = `$docker_cmd -w $(abspath(cwd))`
    end

    # Add in read-only mappings and read-write workspaces
    for shard in shards[2:end]
        outside = realpath(mount_path(shard))
        inside = map_target(shard)
        docker_cmd = `$docker_cmd -v $(outside):$(inside):ro`
    end
    for (outside, inside) in workspaces
        if isdir(outside) || isfile(outside)
            outside = realpath(outside)
        end
        docker_cmd = `$docker_cmd -v $(outside):$inside`
    end

    # Build up environment mappings
    for (k, v) in envs
        docker_cmd = `$docker_cmd -e $k=$v`
    end

    # Finally, return the DockerRunner in all its glory
    return DockerRunner(docker_cmd, workspace_root, shards[1].version, platform)
end

"""
    chown_cleanup(dr::DockerRunner)

On Linux, the user id inside of the docker container doesn't correspond to ours
on the outside, so permissions get all kinds of screwed up.  To fix this, we
have to `chown -R \$(id -u):\$(id -g) \$prefix`, which really sucks, but is
still better than nothing.  This is why we prefer the UserNSRunner on Linux.
"""
function chown_cleanup(dr::DockerRunner; verbose::Bool = false)
    if !Sys.islinux()
        return
    end

    if verbose
        @info("chown'ing prefix back to us...")
    end
    run(`$(sudo_cmd()) chown $(getuid()):$(getgid()) -R $(dr.prefix)`)
end

function Base.run(dr::DockerRunner, cmd, logpath::AbstractString; verbose::Bool = false, tee_stream=stdout)
    did_succeed = true
    oc = OutputCollector(`$(dr.docker_cmd) $(docker_image(dr.rootfs_version)) $(cmd)`; verbose=verbose, tee_stream=tee_stream)
    did_succeed = wait(oc)

    if !isempty(logpath)
        # Write out the logfile, regardless of whether it was successful
        mkpath(dirname(logpath))
        open(logpath, "w") do f
            # First write out the actual command, then the command output
            println(f, cmd)
            print(f, merge(oc))
        end
    end

    # Cleanup permissions, if we need to.
    chown_cleanup(dr; verbose=verbose)

    # Return whether we succeeded or not
    return did_succeed
end

function run_interactive(dr::DockerRunner, cmd::Cmd; stdin = nothing, stdout = nothing, stderr = nothing)
    run_flags = (stdin === nothing && stdout === nothing && stderr === nothing) ? "-ti" : "-i"
    cmd = `$(dr.docker_cmd) $(run_flags) -i $(docker_image(dr.rootfs_version)) $(cmd)`
    if stdin isa AnyRedirectable
        cmd = pipeline(cmd, stdin=stdin)
    end
    if stdout isa AnyRedirectable
        cmd = pipeline(cmd, stdout=stdout)
    end
    if stderr isa AnyRedirectable
        cmd = pipeline(cmd, stderr=stderr)
    end

    if stdout isa IOBuffer
        if !(stdin isa IOBuffer)
            stdin = devnull
        end
        process = open(cmd, "r", stdin)
        @async begin
            while !eof(process)
                write(stdout, read(process))
            end
        end
        wait(process)
    else
        run(cmd)
    end

    # Cleanup permissions, if we need to.
    chown_cleanup(dr)
end


text:
import CLBLAS
import OpenCL

    const clblas = CLBLAS
    clblas.setup()

    const cl = OpenCL

    device, ctx, queue = clblas.get_next_compute_context()

    N = unsigned(5)
    data = Array(Float32, 5)

    data[1] = 5.55555
    data[2] = 3.333
    data[3] = 2.22
    data[4] = 4.4444
    data[5] = 1.1

    bufX = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=data)
    scratchBuff = cl.Buffer(Float32, ctx, :rw, 2 * length(data))
    iMax = cl.Buffer(cl.CL_uint, ctx, :w, 1)

    offiMax = offx = unsigned(0)
    incx = cl.cl_int(1)
    ncq = cl.cl_uint(1)
    clq = queue.id
    newl = cl.cl_uint(0)

    event = cl.UserEvent(ctx, retain=true)
    ptrEvent = [event.id]

    clblas.clblasiSamax( N, iMax.id, offiMax, bufX.id, offx, incx, scratchBuff.id, ncq, [clq], newl, C_NULL, ptrEvent);

    cl.api.clWaitForEvents(cl.cl_uint(1), ptrEvent)

    result = Array(cl.CL_uint, 1)
    cl.enqueue_read_buffer(queue, iMax, result, unsigned(0), nothing, true)

    result = int(floor(result[1]))

    expected = float32(5.55555)
    println(data[result])
    println("------------")
    println(expected)

    try
        @assert(isapprox(norm(expected - data[result]), zero(Float32)))
        info("success!")
    finally
        clblas.teardown()
    end