text:
compiler_env, script, output, base_sysimage = ARGS

if VERSION < v"0.7-"
    error("Unsupported Julia version $VERSION")
end

const Pkg =
    Base.require(Base.PkgId(Base.UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f"), "Pkg"))

Pkg.activate(compiler_env)
@info "Loading PackageCompiler..."
using PackageCompiler

@info "Installing PyCall..."
Pkg.activate(".")
Pkg.add([
    Pkg.PackageSpec(
        name = "PyCall",
        uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0",
    ),
])

if VERSION >= v"1.5-"
    mktempdir() do dir
        tmpimg = joinpath(dir, basename(output))
        @info "Compiling a temporary system image without `PyCall`..."
        create_sysimage(
            Symbol[];
            sysimage_path = tmpimg,
            project = ".",
            base_sysimage = isempty(base_sysimage) ? nothing : base_sysimage,
        )
        @info "Compiling system image..."
        create_sysimage(
            [:PyCall];
            sysimage_path = output,
            project = ".",
            precompile_execution_file = script,
            base_sysimage = tmpimg,
        )
    end
else
    @info "Compiling system image..."
    create_sysimage(
        [:PyCall],
        sysimage_path = output,
        project = ".",
        precompile_execution_file = script,
        base_sysimage = isempty(base_sysimage) ? nothing : base_sysimage,
    )
end

@info "System image is created at $output"

# Notes on two-stage system image monkey-patching for Julia 1.5
#
# Naive `@eval Base`-based monkey-patching stopped working as of Julia 1.5
# presumably because @eval-ing to another module during precompilation now
# throws an error: https://github.com/JuliaLang/julia/pull/35410
#
# We workaround this by creating the system image in two stages:
#
#   1.   The first stage is monkey-patching as done before but without
#       PyCall. This way, we don't get the error because julia does not
#       try to precompile any packages.
#
#   2.   The second stage is the inclusion of PyCall. At this point, we are
#       in a monkey-patched system image. So, it's possible to precompile
#       PyCall now.
#
# Note that even importing PackageCompiler.jl itself inside julia-py is
# problematic because (normally) it has to be precompiled. We avoid this
# problem by running PackageCompiler.jl under --compiled-modules=no. This is
# safe to do thanks to PackageCompiler.do_ensurecompiled
# (https://github.com/JuliaLang/PackageCompiler.jl/blob/3f0f4d882c560c4e4ccc6ab9a8b51ced380bb0d5/src/PackageCompiler.jl#L181-L188)
# using a custom function PackageCompiler.get_julia_cmd
# (https://github.com/JuliaLang/PackageCompiler.jl/blob/3f0f4d882c560c4e4ccc6ab9a8b51ced380bb0d5/src/PackageCompiler.jl#L113-L116)
# (instead of Base.julia_cmd) and ignores --compiled-modules=no of the current
# process.

max_stars_repo_name: dpinol/pyjulia

max_issues_repo_name: dpinol/pyjulia

max_forks_repo_name: dpinol/pyjulia


text:
using Test
using BEAST
using CompScienceMeshes

Γ = meshrectangle(1.0, 1.0, 0.5, 3)
X = raviartthomas(Γ)
SL = derive(TDMaxwell3D.singlelayer(speedoflight=1.0))

Δt, Nt = 1.01, 10
δ = timebasisdelta(Δt,Nt)
h = timebasisc0d1(Δt,Nt)

SLh = (SL, X⊗δ, X⊗h)
A = assemble(SLh...)

@show size(A)

ext: jl


text:
@testset "destination" begin

    p1 = Position([-75, 38.10096062273525])
    res1 = destination(p1, 100, 0).coordinates

    p2 = Position([-75, 39])
    res2 = destination(p2, 100, 180).coordinates

    p3 = Position([12, -54])
    res3 = rhumb_destination(p3, -100, 45).coordinates

    p4 = Position([-75, 38.10096062273525])
    res4 = rhumb_destination(p4, 100, 0).coordinates

    @test round.(res1; digits=6) == [-75.0, 39.000281]
    @test round.(res2; digits=5) == [-75, 38.10068]
    @test res3 == [10.90974456038191, -54.63591552764877]
    @test res4 == [-75, 39.00028098645979]
end


text:
using Flux, Plots
using Base.Iterators: repeated
using Flux: onehotbatch, onecold, crossentropy, throttle, @epochs
using LinearAlgebra

gridsize = 150;
tstart=0.0
tend=20.0
myfunc(x) = sin(x)+0.5*sin(2*x)+0.24*cos(x);
x = collect(range(tstart,stop=tend,length=gridsize));
y= myfunc.(x);

xtest=collect(tend:0.1:2*tend)
xpred=collect(tstart:0.1:2*tend)

data = []
for i in 1:length(x)
    push!(data, ([x[i]], y[i])) #have to input a list
end

Q = 20;
ann = Chain(Dense(1,Q,sin),Dense(Q,1));

function loss(x, y)
    pred=ann(x)
    loss=Flux.mse(ann(x), y)
    return loss
end

opt = ADAM()

ps=params(ann)

for j=1:20
    ann[1].W[j]=j
end
delete!(ps, ann[1].W) #Make the first layers weights (frequencies not trainable)
@epochs 3000 Flux.train!(loss,ps, data, opt)

plot(xpred, ann(xpred')', color=:red, lw=2.0, label="")
scatter!(x,myfunc.(x), color=:blue, legend=false)
scatter!(xtest,myfunc.(xtest), color=:green, lw=2.0, label="", markershape=:cross)
ylims!((-1.5,1.5))
xlims!((tstart,2*tend))
title!("Curve Fitting with Neural Networks")
xlabel!("Input")
ylabel!("Output")

max_stars_repo_name: khannay/CurveFittingNeuralNets

max_issues_repo_name: khannay/CurveFittingNeuralNets

max_forks_repo_name: khannay/CurveFittingNeuralNets


text:
"""
    print_boxed(io, msglines[, prefix="", suffix=""; color=:light_black, chars="[┌│└"])

Print each line in `msglines` to `io` with an enclosing box drawn to
the left, using the specified `color`. The first and last line will
additionally print the optional parameter `prefix` and `suffix`,
respectively. The box drawing characters can be customized by setting
the string `chars` to any combination of four characters.

# Examples

```jldoctest
julia> print_boxed(stdout, ["Hello", "World"], ">>>", "<<<")
┌ >>>  Hello
└  World <<<

julia> print_boxed(stdout, ["Hello world"], chars="(╭│╰")
(  Hello world

julia> print_boxed(stdout, ["Hello", "world"], chars="(╭│╰")
╭  Hello
╰  world
```
"""
function print_boxed(io::IO, msglines, prefix="", suffix="";
                     color=:light_black, chars="[┌│└")
    # This is mostly taken from the print routine behind @info &c from Base
    a,b,c,d = chars
    for (i,msg) in enumerate(msglines)
        boxstr = length(msglines) == 1 ? "$a " :
                 i == 1                ? "$b " :
                 i < length(msglines)  ? "$c " :
                                         "$d "

        printstyled(io, boxstr, bold=true, color=color)
        if i == 1 && !isempty(prefix)
            printstyled(io, prefix, " ", bold=true, color=color)
        end
        print(io, " ", msg)
        if i == length(msglines) && !isempty(suffix)
            !isempty(msg) && print(io, " ")
            printstyled(io, suffix, color=:light_black)
        end
        println(io)
    end
end

"""
    print_boxed(fun::Function, io, args...; kwargs...)

Block-version of [`print_boxed`](@ref) that via
[`redirect_output`](@ref) captures the output of `fun`, and prints
them in a block.

# Examples

```jldoctest
julia> print_boxed(stdout) do io
           println(io, "Hello")
           println(io, "World")
       end
┌  Hello
└  World
```
"""
print_boxed(fun::Function, io::IO, args...; kwargs...) =
    print_boxed(io, split(rstrip(redirect_output(fun)), "\n"), args...; kwargs...)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -4.116064765345723 <= q <= -0.12537054293327043)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.6620590692147021 + -0.9580027962881332*1/(1 + exp(-(-0.38090481240177176 + -0.02939531016329644*1/(1 + exp(-(0.2067810723024852 + 0.7469716803367321*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.3209180228281623*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))) + -0.9724733578803808*1/(1 + exp(-(-0.4274287128450034 + -0.22403315015562253*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.6136110647248842*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3])))))))))) + -0.699289524220764*1/(1 + exp(-(0.4498619401955626 + 0.5463276179102481*1/(1 + exp(-(0.2067810723024852 + 0.7469716803367321*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.3209180228281623*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))) + -0.02415027804893688*1/(1 + exp(-(-0.4274287128450034 + -0.22403315015562253*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.6136110647248842*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))))))))) + 1/(1 + exp(-(-0.8107501452143855 + 0.008470632510253484*1/(1 + exp(-(-0.38090481240177176 + -0.02939531016329644*1/(1 + exp(-(0.2067810723024852 + 0.7469716803367321*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.3209180228281623*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))) + -0.9724733578803808*1/(1 + exp(-(-0.4274287128450034 + -0.22403315015562253*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.6136110647248842*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3])))))))))) + 0.39739497887500974*1/(1 + exp(-(0.4498619401955626 + 0.5463276179102481*1/(1 + exp(-(0.2067810723024852 + 0.7469716803367321*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.3209180228281623*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))) + -0.02415027804893688*1/(1 + exp(-(-0.4274287128450034 + -0.22403315015562253*1/(1 + exp(-(-0.2319263638689142 + 0.46965024077813977*$(x[1]) + 0.028636456496167195*$(x[2]) + 0.15542595955334626*$(x[3])))) + -0.6136110647248842*1/(1 + exp(-(-0.2352871682017641 + -0.6225824541901428*$(x[1]) + -0.2700435280904796*$(x[2]) + -0.7791025613764826*$(x[3]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
ENV["GKSwstype"] = 322 # workaround for gr segfault on GH actions
using Documenter, LowLevelParticleFilters

using Plots
gr()

makedocs(
      sitename = "LowLevelParticleFilters Documentation",
      doctest = false,
      modules = [LowLevelParticleFilters],
      pages = [
            "Home" => "index.md",
            "Benchmark" => "benchmark.md",
            "High-performance distributions" => "distributions.md",
            "API" => "api.md",
      ],
      format = Documenter.HTML(prettyurls = haskey(ENV, "CI")),
) # Due to lots of plots, this will just have to be run on my local machine

deploydocs(
      deps = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-cinder"),
      repo = "github.com/baggepinnen/LowLevelParticleFilters.jl.git",
)


text:
function test_inplace_layer(backend::Backend, T, eps)
  println("-- Testing InplaceLayer on $(typeof(backend)){$T}...")
  registry_reset(backend) # clear layer registry

  ratio = convert(T, 0.6)
  width, height, channels, batch_size = (3,4,5,6)
  input = rand(T, width, height, channels, batch_size)

  layer_data = MemoryDataLayer(tops=[:data], batch_size=batch_size,
      data = Array[input])
  layer_ip1 = InnerProductLayer(name="ip1", tops=[:ip1], bottoms=[:data], output_dim=8)
  layer_dropout1 = DropoutLayer(bottoms=[:ip1])
  layer_dropout2 = DropoutLayer(bottoms=[:ip1])
  layer_ip2 = InnerProductLayer(name="ip2", tops=[:ip2], bottoms=[:ip1], output_dim=1)

  # in random order but dropout1 is put before dropout2
  layers = [layer_ip2, layer_data, layer_ip1, layer_dropout1, layer_dropout2]

  println("    > Setup")
  net = Net("test-inplace", backend, layers)

  # make sure toplogical sort is working properly
  @test net.layers[1] == layer_data
  @test net.layers[2] == layer_ip1
  @test net.layers[3] == layer_dropout1
  @test net.layers[4] == layer_dropout2
  @test net.layers[5] == layer_ip2

  init(net)

  println("    > Forward")
  forward(net)

  ############################################################
  # manually do forward... =.=
  ############################################################
  #-- layer_ip1
  weight = zeros(T, size(net.states[2].W))
  copy!(weight, net.states[2].W)
  output = reshape(weight, size(weight,1),size(weight,2))' *
      reshape(input, width*height*channels, batch_size)

  #-- layer_dropout1
  rand_vals = zeros(T, size(output))
  copy!(rand_vals, net.states[3].rand_vals)
  output .*= (rand_vals .> layer_dropout1.ratio) / (1-layer_dropout1.ratio)

  #-- layer_dropout2
  copy!(rand_vals, net.states[4].rand_vals)
  output .*= (rand_vals .> layer_dropout2.ratio) / (1-layer_dropout2.ratio)
  output = convert(Array{T}, output)

  #-- middle-term check
  got_output = similar(output)
  copy!(got_output, net.states[2].blobs[1])
  @test all(abs.(output - got_output) .< eps)

  #-- layer_ip2
  weight = zeros(T, size(net.states[5].W))
  copy!(weight, net.states[5].W)
  output = reshape(weight, size(weight,1), size(weight,2))' * output
  got_output = similar(output)
  copy!(got_output, net.states[5].blobs[1])
  @test all(abs.(output - got_output) .< eps)

  println("    > Backward")
  top_diff = rand(T, size(output))
  copy!(net.states[5].blobs_diff[1], top_diff)
  backward(net)
  ############################################################
  # manually do backward... =.=
  ############################################################
  #-- layer_ip2
  weight = zeros(T, size(net.states[5].W))
  copy!(weight, net.states[5].W)
  grad = reshape(weight, size(weight,1), size(weight,2)) * top_diff

  #-- layer_dropout2
  copy!(rand_vals, net.states[4].rand_vals)
  grad .*= (rand_vals .> layer_dropout2.ratio) / (1-layer_dropout2.ratio)

  #-- layer_dropout1
  copy!(rand_vals, net.states[3].rand_vals)
  grad .*= (rand_vals .> layer_dropout1.ratio) / (1-layer_dropout1.ratio)
  grad = convert(Array{T}, grad)

  got_grad = similar(grad)
  copy!(got_grad, net.states[2].blobs_diff[1])
  @test all(abs.(grad - got_grad) .< eps)

  destroy(net)
end
function test_inplace_layer(backend::Backend)
  test_inplace_layer(backend, Float64, 1e-10)
  test_inplace_layer(backend, Float32, 1e-4)
end

if test_cpu
  test_inplace_layer(backend_cpu)
end
if test_gpu
  test_inplace_layer(backend_gpu)
end


text:
@testset "Moduels" begin

    Symbol_array =[:sym, :sym2, :sym3]
    String_array = ["a", "b", "c"]

    mt_a = Mod("a")
    mt_b = Mod("b")
    mt_c = Mod("c")
    mt_d = Mod("d")
    mt_e = Mod("e")

    mt_symb = Mod(:sym)

    moduel_test_Symbol_array_non_unique =[:sym, :sym2, :sym3, :sym2]
    moduel_test_String_array_non_unique= ["a", "b", "c" ,"a"]

    @testset "SetMod_single"  begin
        @test mt_symb.m == :sym
        @test mt_a.m == "a"
        @test mt_symb |> length  == 1
        @test mt_a |> length == 1
    end

    Group_Mod_test = Group_Mod([mt_a,mt_b,mt_c,mt_d,mt_e])
    Group_Mod_test_order = Group_Mod([mt_b,mt_a,mt_c,mt_d,mt_e,mt_c])

    @testset "Group_mod" begin
        @test Group_Mod_test_order.m == Group_Mod_test.m
    end

    @testset "group_mod" begin
        @test group_mod(["b", "a", "c","c","e","d"]).m == Group_Mod([mt_a,mt_b,mt_c,mt_d,mt_e]).m
        @test group_mod([:b, :a, :c,:c,:e,:d]).m == Group_Mod([Mod(:a),Mod(:b),Mod(:c),Mod(:d),Mod(:e)]).m
    end

    mt_String_array = group_mod(["b", "a", "c","c","e","d"])
    @testset "length" begin
        @test length(mt_String_array) == 5
        @test length(Mod("a")) == 1
    end

    @testset "isless" begin
        @test isless(Mod(:b),Mod(:a)) == false
        @test isless(Mod(:c),Mod(:e)) == true
        @test isless(Mod("a"),Mod("begin")) == true
    end

    @testset "isless" begin
        @test isless(Mod(:b),Mod(:a)) == false
        @test isless(Mod(:c),Mod(:e)) == true
        @test isless(Mod("a"),Mod("begin")) == true
    end

    test_collect = collect(Group_Mod_test_order)
    @testset "collect" begin
        @test length(test_collect) == 5
        @test eltype(test_collect) == Mod{String}
    end

end


text:
module Variogram
using Distances
using SmoothingSplines
using Plots

export variogram
export getdecorrelationlength

"""
    Utility type to collect values for a semivariogram bin
"""
type BuildBin
    sumofsquares::Float64
    entries::Int

    BuildBin() = new(0.0, 0)
end

"""
   Add a value to a mean
"""
function add(x::BuildBin, value::Float64)
    x.sumofsquares += value^2
    x.entries += 1
end

"""
   Calculate the final value of a mean
"""
function mean(x::BuildBin)
    result::Float64 = NaN

    if x.entries > 0
        result = x.sumofsquares / (2 * x.entries)

    end

    return result
end

"""
    Computes a variogram for a set of geographical points.

    The passed in `pointdata` must be an array in which is row
    contains `lon`, `lat`, `value`.

    `binsize` defines the size of the bins in the variogram,
    in kilometres.

"""
function variogram(pointdata::Array{Float64,2}, binsize::Int64)

    @assert length(pointdata[1, :]) == 3 "Point data must contain 3 columns"
    @assert binsize > 0 "Bin size must be greater than zero"

    # Create the array of bins.
    # The Earth is 40,075km in diameter, so the maximum
    # distance between points is 20,037.5km.
    # So the number of bins is ceil(20,037.5 / binsize).
    bins::Array{BuildBin} = [BuildBin() for i = 1:ceil(20037.5 / binsize)]
    bin_distances::Array{Int64} = collect(0:binsize:20037.5)

    # Go through every pair of points, and
    # add the difference in value between each pair
    # to the appropriate bin
    npoints::Int64 = size(pointdata)[1]

    largest_bin::Int64 = 0

    for i in 1:npoints - 1

        point1::Array{Float64} = pointdata[i, 1:2]
        value1::Float64 = pointdata[i, 3]

        for j in i:npoints

            point2::Array{Float64} = pointdata[j, 1:2]

            distance::Float64 = evaluate(Distances.Haversine(6371), point1, point2)
            bin::Int64 = Int64(floor(distance / binsize)) + 1
            if bin > largest_bin
                largest_bin = bin
            end
            add(bins[bin], abs(pointdata[j, 3] - value1))
        end
    end

    bin_means::Array{Float64} = [mean(i) for i in bins]

    return [bin_distances[1:largest_bin] bin_means[1:largest_bin]]
end

"""
    Get the decorrelation length of a variogram.

    Returns -1 if no decorrelation length can be Calculated
"""
function getdecorrelationlength(variogram::Array{Float64, 2}, plot::Bool)

    local decorrelationlength::Int64 = -1
    fillgaps!(variogram)

    if size(variogram)[1] > 2
        local spline = SmoothingSplines.fit(SmoothingSplines.SmoothingSpline, variogram[:,1], variogram[:,2], 1e6)
        local smoothed = SmoothingSplines.predict(spline)

        if plot
            Plots.plot!(variogram[:,1], smoothed, linewidth=2)
        end

        limitindex::Int64 = getdecorrelationlimit(smoothed)
        if limitindex != -1
            decorrelationlength = variogram[limitindex, 1]
        end
    end

    return decorrelationlength
end

"""
    Fill gaps in a variogram using linear interpolation
"""
function fillgaps!(variogram::Array{Float64,2})

    local lastlag::Int64 = 0
    local currentlag::Int64 = 1

    while currentlag < size(variogram)[1]

        if !isnan(variogram[currentlag, 2])
            lastlag = currentlag
            currentlag +=  1
        else
            local nextgoodlag::Int64 = currentlag + 1
            while isnan(variogram[nextgoodlag, 2])
                nextgoodlag = nextgoodlag + 1
            end

            local gaplength::Int64 = nextgoodlag - lastlag
            local gaprange::Float64 = variogram[nextgoodlag, 2] - variogram[lastlag, 2]
            local stepsize::Float64 = gaprange / gaplength

            local fillstep::Int64 = 0
            for filllag in currentlag:nextgoodlag - 1
                fillstep += 1
                variogram[filllag, 2] = variogram[lastlag, 2] + (stepsize * fillstep)
            end

            currentlag = nextgoodlag + 1
            lastlag = currentlag - 1
        end
    end
end

"""
    Calculate the decorrelation limit for a variogram
    as the first step whose delta is less than 5%
    of the first positive delta of the varioagram.
"""
function getdecorrelationlimit(spline::Array{Float64, 1})

    local deltas::Array{Float64} = [(spline[i + 1] - spline[i]) for i in 1:length(spline) - 1]

    local position::Int64 = 1
    local firstpositivedelta::Float64 = 0

    # Find the first positive delta
    while firstpositivedelta == 0 && position <= length(deltas)
        if deltas[position] > 0
            firstpositivedelta = deltas[position]
        end
        position += 1
    end

    # Find the first delta that's smaller than 5% of the first positive delta
    local decorrelationposition::Int64 = -1

    while decorrelationposition == -1 && position < length(deltas)
        if deltas[position] <= 0 || deltas[position] / firstpositivedelta <= 0.05
            decorrelationposition = position
        else
            position += 1
        end
    end

    return decorrelationposition
end

end #Module


text:
is_constant_cache(cache::OrdinaryDiffEqConstantCache) = true
is_constant_cache(cache::OrdinaryDiffEqCache) = false
is_constant_cache(cache::CompositeCache) = is_constant_cache(cache.caches[1])

function DiffEqBase.unwrap_cache(integrator::ODEIntegrator, is_stiff)
  alg   = integrator.alg
  cache = integrator.cache
  iscomp = alg isa CompositeAlgorithm
  if !iscomp
    return cache
  elseif alg.choice_function isa AutoSwitch
    num = is_stiff ? 2 : 1
    return cache.caches[num]
  else
    return cache.caches[integrator.cache.current]
  end
end

@deprecate alg_cache(alg::Union{OrdinaryDiffEqAlgorithm,DAEAlgorithm},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{iip}}) where iip alg_cache(alg,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,Val(iip))


text:
#MDDatasets show functions
#-------------------------------------------------------------------------------

#==DataF1 show functions
===============================================================================#
#show() a compact version of a vector:
function _showcompact(io::IO, x::Vector{T}) where T<:Number
	maxelem = 10 #WANTCONST
	if length(x)>maxelem
		print(io, "[")
		for i in 1:(maxelem-3)
			print(io, x[i], ",")
		end
		print(io, "...")
		for v in x[end-1:end]
			print(io, ",", v)
		end
		print(io, "]")
	else
		show(io, x)
	end
end

#Don't show module name/subtypes for DataF1:
function Base.show(io::IO, ds::DataF1{TX,TY}) where {TX<:Number, TY<:Number}
	print(io, "DataF1(x=")
		_showcompact(io, ds.x)
		print(io, ",y=")
		_showcompact(io, ds.y)
		print(io, ")")
end

#==Misc show functions
===============================================================================#
#==
function Base.show(io::IO, d::DataScalar{T}) where T
	print(io, "DataScalar{$T}($(d.v))")
end
==#

function Base.show(io::IO, p::Point2D{TX,TY}) where {TX,TY}
	print(io, "Point2D{$TX,$TY}($(p.x), $(p.y))")
end

function Base.show(ctx::IOContext, p::Point2D)
	if get(ctx.dict, :compact, false) #Show compact:
		print(io, "pt($(p.x), $(p.y))")
	else
		show(ctx.io, p)
	end
end

function Base.show(io::IO, x::CrossType)
	proplist = String[]
	x.v&XINGTYPE_RISE != 0 ? push!(proplist, "rise") : nothing
	x.v&XINGTYPE_FALL != 0 ? push!(proplist, "fall") : nothing
	x.v&XINGTYPE_SING != 0 ? push!(proplist, "sing") : nothing
	x.v&XINGTYPE_FLAT != 0 ? push!(proplist, "flat") : nothing
	x.v&XINGTYPE_THRU != 0 ? push!(proplist, "thru") : nothing
	x.v&XINGTYPE_REV != 0 ? push!(proplist, "rev") : nothing
	x.v&XINGTYPE_FIRSTLAST != 0 ? push!(proplist, "firstlast") : nothing

	print(io, "CrossType(", join(proplist, ", "), ")")
end

#Last line


text:
abstract type AbstractConnectivityModel end

@inline (mdl::AbstractConnectivityModel)(n_vertices, n_pins, n_nets, k) = mdl(n_vertices, n_pins, n_nets)

struct AffineConnectivityModel{Tv} <: AbstractConnectivityModel
    α::Tv
    β_vertex::Tv
    β_pin::Tv
    β_net::Tv
end

function AffineConnectivityModel(; α = false, β_vertex = false, β_pin = false, β_net = false)
    AffineConnectivityModel(promote(α, β_vertex, β_pin, β_net)...)
end

@inline cost_type(::Type{AffineConnectivityModel{Tv}}) where {Tv} = Tv

(mdl::AffineConnectivityModel)(n_vertices, n_pins, n_nets) = mdl.α + n_vertices * mdl.β_vertex + n_pins * mdl.β_pin + n_nets * mdl.β_net

function bound_stripe(A::SparseMatrixCSC, K, mdl::AffineConnectivityModel)
    return bound_stripe(A, K, oracle_stripe(StepHint(), mdl, A))
end
function bound_stripe(A::SparseMatrixCSC, K, ocl::AbstractOracleCost{<:AffineConnectivityModel})
    m, n = size(A)
    N = nnz(A)
    mdl = oracle_model(ocl)
    @assert mdl.β_vertex >= 0
    @assert mdl.β_pin >= 0
    @assert mdl.β_net >= 0
    c_hi = ocl(1, n + 1)
    c_lo = mdl.α + fld(c_hi - mdl.α, K)
    return (c_lo, c_hi)
end

struct ConnectivityOracle{Ti, Net, Mdl} <: AbstractOracleCost{Mdl}
    pos::Vector{Ti}
    net::Net
    mdl::Mdl
end

oracle_model(ocl::ConnectivityOracle) = ocl.mdl

function oracle_stripe(hint::AbstractHint, mdl::AbstractConnectivityModel, A::SparseMatrixCSC{Tv, Ti}; net=nothing, adj_A=nothing, kwargs...) where {Tv, Ti}
    @inbounds begin
        m, n = size(A)
        pos = A.colptr
        if net === nothing
            net = netcount(hint, A; kwargs...)
        end
        return ConnectivityOracle(pos, net, mdl)
    end
end

@inline function (cst::ConnectivityOracle{Ti, Mdl})(j, j′, k...) where {Ti, Mdl}
    @inbounds begin
        w = cst.pos[j′] - cst.pos[j]
        d = cst.net[j, j′]
        return cst.mdl(j′ - j, w, d, k...)
    end
end

@inline function (stp::Step{Ocl})(_j, _j′, _k...) where {Ti, Mdl, Ocl <: ConnectivityOracle{Ti, Mdl}}
    @inbounds begin
        cst = stp.ocl
        j = destep(_j)
        j′ = destep(_j′)
        k = maptuple(destep, _k...)
        w = cst.pos[j′] - cst.pos[j]
        d = Step(cst.net)(_j, _j′)
        return cst.mdl(j′ - j, w, d, k...)
    end
end

#=
function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::DynamicTotalChunker{F}, args...; kwargs...) where {F<:AbstractConnectivityModel, Tv, Ti}
    return pack_stripe(A, DynamicTotalChunker(ConstrainedCost(method.f, FeasibleCost(), Feasible())), args..., kwargs...)
end

function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::DynamicTotalChunker{<:ConstrainedCost{F}}, args...; n_nets = nothing, kwargs...) where {F<:AbstractConnectivityModel, Tv, Ti}
    @inbounds begin
        # matrix notation...
        # i = 1:m rows, j = 1:n columns
        m, n = size(A)

        A_pos = A.colptr
        A_idx = A.rowval

        f = method.f.f
        w = oracle_stripe(StepHint(), method.f.w, A, args...) #TODO reverse outer loop so step oracles can have stationary j′
        w_max = method.f.w_max

        Δ_net = zeros(Int, n + 1) # Δ_net is the number of additional distinct entries we see as our part size grows.
        hst = fill(n + 1, m) # hst is the last time we saw some nonzero
        cst = Vector{cost_type(f)}(undef, n + 1) # cst[j] is the best cost of a partition from j to n
        spl = Vector{Int}(undef, n + 1)
        Δ_net[n + 1] = 0
        cst[n + 1] = zero(cost_type(f))
        for j = n:-1:1
            d = A_pos[j + 1] - A_pos[j] # The number of distinct nonzero blocks in each candidate part
            Δ_net[j] = d
            for q = A_pos[j] : A_pos[j + 1] - 1
                i = A_idx[q]
                j′ = hst[i]
                Δ_net[j′] -= 1
                hst[i] = j
            end
            best_c = cst[j + 1] + f(1, d, d)
            best_d = d
            best_j′ = j + 1
            for j′ = j + 2 : n + 1
                if w(j, j′) > w_max
                    break
                end
                d += Δ_net[j′ - 1]
                c = cst[j′] + f(j′ - j, A_pos[j′] - A_pos[j], d)
                if c < best_c
                    best_c = c
                    best_d = d
                    best_j′ = j′
                end
            end
            cst[j] = best_c
            spl[j] = best_j′
        end

        K = 0
        j = 1
        while j != n + 1
            j′ = spl[j]
            K += 1
            spl[K] = j
            j = j′
        end
        spl[K + 1] = j
        resize!(spl, K + 1)
        return SplitPartition{Ti}(K, spl)
    end
end
=#


text:
println("Hello EPL")

@Base.kwdef mutable struct Executor
    num_sensors::Int = 72
    mc_values::Vector{Int64} = []
    gc_stored_odors::Vector{Vector{Int64}} = [[]]
    mc_to_gc_fanin::Int
    max_test_iters::Int
    minority_vote_ratio::Float16
end

exec = Executor(num_sensors=72, mc_to_gc_fanin=7, max_test_iters=3, minority_vote_ratio=0.03)
# print(typeof(exec.gc_stored_odors))

function train(exec::Executor, train_data)
    exec.gc_stored_odors = train_data
end

train_data = [[3, 4, 5], [7, 8, 9], [2, 4, 6]]
train(exec, train_data)
println("stored odors: ", exec.gc_stored_odors)
println("num_sensors: $(exec.num_sensors)")

max_stars_repo_name: epl-actors/epl-julia

max_issues_repo_name: epl-actors/epl-julia

max_forks_repo_name: epl-actors/epl-julia


text:
import Statistics
using Test, Random
using NiLang, NiLang.AD
using Distributions

@testset "statistics" begin
    x = randn(100)
    y = randn(100)
    @test i_mean_sum(0.0, 0.0, x)[1] ≈ Statistics.mean(x)
    @test all(i_var_mean_sum(0.0, 0.0, 0.0, 0.0, copy(x)) .≈ (Statistics.var(x), Statistics.var(x)*99, Statistics.mean(x), sum(x), x))
    @test all(isapprox.((~i_var_mean_sum)(i_var_mean_sum(0.0, 0.0, 0.0, 0.0, copy(x))...), (0.0, 0.0, 0.0, 0.0, x), atol=1e-8))
    @test all(i_cor_cov(0.0, 0.0, copy(x), copy(y)) .≈ (Statistics.cor(x,y), Statistics.cov(x,y), x, y))
    @test all(isapprox.((~i_cor_cov)(i_cor_cov(0.0, 0.0, copy(x), copy(y))...), (0.0, 0.0, x, y), atol=1e-8))
end

@testset "normal log pdf" begin
    out = 0.0
    x = 1.0
    μ = 0.3
    σ = 1.5
    l1 = i_normal_logpdf(out, x, μ, σ)
    distri = Normal(μ, σ)
    l2 = logpdf(distri, x)
    @test l1[1] ≈ l2
    @test check_inv(i_normal_logpdf, (out, x, μ, σ))
end


text:
"""
Based on
https://hyperopt.github.io/hyperopt/
"""
module TestQuantUniform

using Test
using TreeParzen

# Simple quantuniform
# Minimize the objective over the space
space = Dict(:x => HP.QuantUniform(:x, 0.0, 1.0, 2.0))
best = fmin(params -> params[:x]^2, space, 100)
@test 0 <= best[:x] <= 1

end # module TestQuantUniform
true


text:
using Test

include("nthprime.jl")
include("example.jl")

@testset "first prime" begin
    @test main(1) == 2
end

@testset "no zeroth" begin
    @test main(0) == "there is no zeroth prime"
end

@testset "second prime" begin
    @test main(2) == 3
end

@testset "sixth prime" begin
    @test main(6) == 13
end

@testset "big prime" begin
    @test main(10001) == 104743
end

max_stars_repo_name: KrymoreAkbota/julia

max_issues_repo_name: KrymoreAkbota/julia

max_forks_repo_name: KrymoreAkbota/julia


text:
using Test
using Sirion

function testBackAndForthPressureEnergy()
    d_list = [10.0^n for n in -5:5]
    p_list = [10.0^n for n in -5:5]
    eos = Sirion.IdealGas(5.0/3.0)
    for d in d_list
        for p in p_list
            e = Sirion.calcSpecificThermalEnergy(eos, d, p)
            pr = Sirion.calcPressure(eos, d, e)
            if !(pr ≈ p)
                return false
            end
        end
    end
    return true
end

@test testBackAndForthPressureEnergy()

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
"""
    get_value(obj::AbstractObject)

Returns the value saved by `obj`
"""
function get_value(obj::AbstractObject)
    return obj.result[1]
end

"""
    val(x)

`val` is just a short-hand for [`get_value`](@ref)
"""
val(x) = get_value(x)

get_position(p::Point) = p
get_position(t::Transformation) = t.point

"""
    get_position(obj::Object)

Get access to the position that got saved in a previous object.

# Returns
- `Point`: the point stored by a previous object.

# Throws
- If the function of Object didn't return a Point or Transformation
"""
function get_position(obj::Object)
    return get_position(obj.result[1])
end

"""
    pos(x)

`pos` is just a short-hand for [`get_position`](@ref)
"""
pos(x) = get_position(x)

# As it is just the number tuple -> return it
get_scale(x::Tuple{<:Number,<:Number}) = Scale(x...)

# If just the number -> return it as a tuple
get_scale(x::Number) = Scale(x, x)

get_scale(s::Scale) = s
get_scale(t::Transformation) = Scale(t.scale.x, t.scale.y)

function get_scale(s::Symbol)
    @assert s == :current_scale
    cs = get_current_setting()
    return cs.current_scale
end

"""
    get_scale(obj::AbstractObject)

Get access to the scaling that got saved in a previous object.

# Returns
- `Scaling`: the scale stored by a previous object.
"""
function get_scale(obj::AbstractObject)
    return get_scale(obj.result[1])
end

"""
    scl(x)

`scl` is just a short-hand for [`get_scale`](@ref)
"""
scl(x) = get_scale(x)

get_angle(x::Number) = x
get_angle(t::Transformation) = t.angle

"""
    get_angle(obj::AbstractObject)

Get access to the angle that got saved in a previous object.

# Returns
- `Float64`: the angle stored by a previous object i.e via `return Transformation(p, angle)`
"""
function get_angle(obj::AbstractObject)
    return get_angle(obj.result[1])
end

"""
    ang(x)

`ang` is just a short-hand for [`get_angle`](@ref)
"""
ang(x) = get_angle(x)


text:
"""
install a given target. Currently only supports "path",
future will support installing a specific Julia version
and Comonicon applications.

# Arguments

- `target`: install given target, can be "path".
"""
@cast function install(target::String)
    if target == "path"
        comonicon_install_path()
    else
        error("unkown target: $target")
    end
    return
end


text:
module ComputeGraphs

# command definitions
include("cmds.jl")
# compute graph data structure
include("tree.jl")
# tree optimisation functions
include("tree_opt.jl")
# tree statistics
include("tree_stats.jl")
# functions to parse dsl file
include("dsl.jl")

end


text:
using PyPlot
using Jagot.plotting
plot_style("ggplot")

using PyCall
Cycler = pyimport("cycler")
plt.rc("axes", prop_cycle=(plt.rcParams["axes.prop_cycle"] +
                           Cycler.cycler("linestyle", ["-","--",":","-.",".","-","--"])))

using Statistics
using Random

using Atoms
using AtomicLevels
using AtomicPotentials
using SCF

using LinearAlgebra
using FiniteDifferencesQuasi
using FEDVRQuasi
using IntervalSets

function mean_color(color::String)
    c = parse(Colorant, color)
    mean([c.r,c.g,c.b])
end

lerp(a,b,t) = (1-t)*a + t*b

mean_position(x, ϕ) = ϕ'*Diagonal(x)*ϕ/(ϕ'ϕ)

function cfigure(fun::Function, figname; clear=true, tight=true, kwargs...)
    figure(figname; kwargs...)
    clear && clf()
    fun()
    tight && tight_layout()
end

function csubplot(fun::Function, args...; nox=false, noy=false)
    ax = subplot(args...)
    fun()
    if nox
        no_tick_labels(:x)
        xlabel("")
    end
    if noy
        no_tick_labels(:y)
        ylabel("")
    end
    ax
end

function plot_orbitals(atom::Atom; nrplot=1000, plot_basis=false,
                       fig_name=string(atom.potential))
    R = radial_basis(atom)
    d = axes(R,1).domain
    r = range(leftendpoint(d), stop=rightendpoint(d), length=nrplot)
    χ = R[r, :]

    n = 0
    ℓ = 0
    i = 0
    orbitals = map(atom.orbitals) do o
        on,oℓ = Atoms.getn(o),Atoms.getℓ(o)
        if on == n && oℓ == ℓ
            i
        else
            n = on
            ℓ = oℓ
            i += 1
        end
    end
    m = orbitals[end] + plot_basis

    cfigure("$(fig_name)") do
        for i in 1:m
            csubplot(m,1,i, nox=i<m) do
                js = findall(isequal(i), orbitals)
                for j in js
                    plot(r, abs2.(χ*view(atom, j).args[2]), label="$(atom.orbitals[j])")
                end
                legend(framealpha=0.75,ncol=2)
                xscale("log")
                i == m && xlabel(L"$r$ [au]")
                spatial_orbs = unique([Atoms.getspatialorb(atom.orbitals[j])
                                       for j in js])
                ylabel(join(string.(spatial_orbs), ", "))
            end
        end
        if plot_basis
            csubplot(m,1,m) do
                plot(r, χ)
                xscale("log")
                xlabel(L"$r$ [au]")
            end
        end
    end
end

function get_atom_grid(grid_type, rₘₐₓ, ρ, nucleus; fedvr_order=10)
    Z = charge(nucleus)
    amend_order=nucleus isa PointCharge # For correct boundary
                                        # conditions at r=0.

    if grid_type == :fedvr
        # FEDVR is more accurate, but can be expensive to use
        N = max(ceil(Int, rₘₐₓ/(ρ*fedvr_order)),2)
        t = range(0.0, stop=rₘₐₓ, length=N)
        amended_order = vcat(fedvr_order+5, fill(fedvr_order,length(t)-2))
        FEDVR(t, amend_order ? amended_order : fedvr_order)[:,2:end-1]
    else
        # Finite-differences are much lighter, but may require very
        # fine grids to converge.
        N = ceil(Int, rₘₐₓ/ρ + 1/2)
        args = amend_order ? () : (zero(ρ),)
        RadialDifferences(N, ρ, Z, args...)
    end
end

function hydrogen()
    nucleus = pc"H"
    R = get_atom_grid(:fedvr, 10.0, 0.1, nucleus)
    gst = ground_state(nucleus)
    atom = Atom(R, [spin_configurations(gst)[1]], nucleus)
    plot_orbitals(atom)
    savefig("docs/src/figures/hydrogen.svg")
end

function helium()
    nucleus = pc"He"
    R = get_atom_grid(:fedvr, 10.0, 0.1, nucleus)
    gst = ground_state(nucleus)
    atom = Atom(R, [spin_configurations(gst)[1]], nucleus)
    fock = Fock(atom)
    optimize!(fock)
    plot_orbitals(atom)
    begin
        hydrogenics = [Atom(R, [spin_configurations(ground_state(n))[1]], n)
                       for n in [pc"H", pc"He"]]
        d = axes(R,1).domain
        r = range(leftendpoint(d), stop=rightendpoint(d), length=1000)
        χ = R[r, :]
        for (hyd,label) in zip(hydrogenics, ["H 1s", "He+ 1s"])
            plot(r, abs2.(χ*view(hyd, 1).args[2]), label=label)
        end
        legend(framealpha=0.75,ncol=2)
    end
    savefig("docs/src/figures/helium.svg")
end

function beryllium(grid_type)
    nucleus = pc"Be"
    R = get_atom_grid(grid_type, 15.0, 0.1, nucleus)
    gst = ground_state(nucleus)
    atom = Atom(R, [spin_configurations(gst)[1]], nucleus)
    fock = Fock(atom)
    optimize!(fock,ω=0.999,ωmax=1-1e-3,scf_method=:arnoldi)
    plot_orbitals(atom)
    savefig("docs/src/figures/beryllium-$(grid_type).svg")
end

macro echo(expr)
    println(expr)
    :(@time $expr)
end

@info "Documentation plots"
mkpath("docs/src/figures")
@echo hydrogen()
@echo helium()
@echo beryllium(:fedvr)
@echo beryllium(:fd)


text:
zero(::Type{Sonum8}) = Sonum8(0x00)
zero(::Type{Sonum16}) = Sonum16(0x0000)

one(::Type{Sonum8}) = Sonum8(1.0)
one(::Type{Sonum16}) = Sonum16(1.0)

minusone(::Type{Sonum8}) = Sonum8(-1.0)
minusone(::Type{Sonum16}) = Sonum16(-1.0)

notareal(::Type{Sonum8}) = Sonum8(0x80)
notareal(::Type{Sonum16}) = Sonum16(0x8000)
notareal(x::AbstractSonum) = notareal(typeof(x))

floatmin(::Type{Sonum8}) = Sonum8(0x01)
floatmax(::Type{Sonum8}) = Sonum8(0x7F)

floatmin(::Type{Sonum16}) = Sonum16(0x0001)
floatmax(::Type{Sonum16}) = Sonum16(0x7FFF)
floatmin(x::AbstractSonum) = floatmin(typeof(x))
floatmax(x::AbstractSonum) = floatmax(typeof(x))

isfinite(x::T) where {T<:AbstractSonum} = x != notareal(T)


text:
module JuliaSandbox

include("internal_module.jl")

end

max_stars_repo_name: erwanlecarpentier/JuliaSandbox

max_issues_repo_name: erwanlecarpentier/JuliaSandbox

max_forks_repo_name: erwanlecarpentier/JuliaSandbox


text:
@testset "minimal_overlap_test" begin

    @test_broken MISSING_overlap(emptyinterval(), emptyinterval()) === bothEmpty

    @test_broken MISSING_overlap(emptyinterval(), interval(1.0,2.0)) === firstEmpty

    @test_broken MISSING_overlap(interval(1.0,2.0), emptyinterval()) === secondEmpty

    @test_broken MISSING_overlap(interval(-Inf,2.0), interval(3.0,Inf)) === before

    @test_broken MISSING_overlap(interval(-Inf,2.0), interval(3.0,4.0)) === before

    @test_broken MISSING_overlap(interval(2.0,2.0), interval(3.0,4.0)) === before

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(3.0,4.0)) === before

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(3.0,3.0)) === before

    @test_broken MISSING_overlap(interval(2.0,2.0), interval(3.0,3.0)) === before

    @test_broken MISSING_overlap(interval(2.0,2.0), interval(3.0,Inf)) === before

    @test_broken MISSING_overlap(interval(-Inf,2.0), interval(2.0,3.0)) === meets

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(2.0,3.0)) === meets

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(2.0,Inf)) === meets

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(1.5,2.5)) === overlaps

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(1.0,Inf)) === starts

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(1.0,3.0)) === starts

    @test_broken MISSING_overlap(interval(1.0,1.0), interval(1.0,3.0)) === starts

    @test_broken MISSING_overlap(interval(1.0,2.0), entireinterval()) === containedBy

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(-Inf,3.0)) === containedBy

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(0.0,3.0)) === containedBy

    @test_broken MISSING_overlap(interval(2.0,2.0), interval(0.0,3.0)) === containedBy

    @test_broken MISSING_overlap(interval(2.0,2.0), interval(0.0,Inf)) === containedBy

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(-Inf,2.0)) === finishes

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(0.0,2.0)) === finishes

    @test_broken MISSING_overlap(interval(2.0,2.0), interval(0.0,2.0)) === finishes

    @test_broken MISSING_overlap(interval(1.0,2.0), interval(1.0,2.0)) === equals

    @test_broken MISSING_overlap(interval(1.0,1.0), interval(1.0,1.0)) === equals

    @test_broken MISSING_overlap(interval(-Inf,1.0), interval(-Inf,1.0)) === equals

    @test_broken MISSING_overlap(entireinterval(), entireinterval()) === equals

    @test_broken MISSING_overlap(interval(3.0,4.0), interval(2.0,2.0)) === after

    @test_broken MISSING_overlap(interval(3.0,4.0), interval(1.0,2.0)) === after

    @test_broken MISSING_overlap(interval(3.0,3.0), interval(1.0,2.0)) === after

    @test_broken MISSING_overlap(interval(3.0,3.0), interval(2.0,2.0)) === after

    @test_broken MISSING_overlap(interval(3.0,Inf), interval(2.0,2.0)) === after

    @test_broken MISSING_overlap(interval(2.0,3.0), interval(1.0,2.0)) === metBy

    @test_broken MISSING_overlap(interval(2.0,3.0), interval(-Inf,2.0)) === metBy

    @test_broken MISSING_overlap(interval(1.5,2.5), interval(1.0,2.0)) === overlappedBy

    @test_broken MISSING_overlap(interval(1.5,2.5), interval(-Inf,2.0)) === overlappedBy

    @test_broken MISSING_overlap(interval(1.0,Inf), interval(1.0,2.0)) === startedBy

    @test_broken MISSING_overlap(interval(1.0,3.0), interval(1.0,2.0)) === startedBy

    @test_broken MISSING_overlap(interval(1.0,3.0), interval(1.0,1.0)) === startedBy

    @test_broken MISSING_overlap(interval(-Inf,3.0), interval(1.0,2.0)) === contains

    @test_broken MISSING_overlap(entireinterval(), interval(1.0,2.0)) === contains

    @test_broken MISSING_overlap(interval(0.0,3.0), interval(1.0,2.0)) === contains

    @test_broken MISSING_overlap(interval(0.0,3.0), interval(2.0,2.0)) === contains

    @test_broken MISSING_overlap(interval(-Inf,2.0), interval(1.0,2.0)) === finishedBy

    @test_broken MISSING_overlap(interval(0.0,2.0), interval(1.0,2.0)) === finishedBy

    @test_broken MISSING_overlap(interval(0.0,2.0), interval(2.0,2.0)) === finishedBy

end

@testset "minimal_overlap_dec_test" begin

    @test_broken MISSING_overlap(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(emptyinterval(), trv)) === bothEmpty

    @test_broken MISSING_overlap(DecoratedInterval(emptyinterval(), trv), DecoratedInterval(interval(1.0,2.0), com)) === firstEmpty

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(emptyinterval(), trv)) === secondEmpty

    @test_broken MISSING_overlap(DecoratedInterval(interval(2.0,2.0), def), DecoratedInterval(interval(3.0,4.0), def)) === before

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), dac), DecoratedInterval(interval(3.0,4.0), com)) === before

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), com), DecoratedInterval(interval(3.0,3.0), trv)) === before

    @test_broken MISSING_overlap(DecoratedInterval(interval(2.0,2.0), trv), DecoratedInterval(interval(3.0,3.0), def)) === before

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(2.0,3.0), def)) === meets

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), dac), DecoratedInterval(interval(1.5,2.5), def)) === overlaps

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(1.0,3.0), com)) === starts

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,1.0), trv), DecoratedInterval(interval(1.0,3.0), def)) === starts

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(0.0,3.0), dac)) === containedBy

    @test_broken MISSING_overlap(DecoratedInterval(interval(2.0,2.0), trv), DecoratedInterval(interval(0.0,3.0), def)) === containedBy

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), trv), DecoratedInterval(interval(0.0,2.0), com)) === finishes

    @test_broken MISSING_overlap(DecoratedInterval(interval(2.0,2.0), def), DecoratedInterval(interval(0.0,2.0), dac)) === finishes

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,2.0), def), DecoratedInterval(interval(1.0,2.0), def)) === equals

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,1.0), dac), DecoratedInterval(interval(1.0,1.0), dac)) === equals

    @test_broken MISSING_overlap(DecoratedInterval(interval(3.0,4.0), trv), DecoratedInterval(interval(2.0,2.0), trv)) === after

    @test_broken MISSING_overlap(DecoratedInterval(interval(3.0,4.0), def), DecoratedInterval(interval(1.0,2.0), def)) === after

    @test_broken MISSING_overlap(DecoratedInterval(interval(3.0,3.0), com), DecoratedInterval(interval(1.0,2.0), dac)) === after

    @test_broken MISSING_overlap(DecoratedInterval(interval(3.0,3.0), def), DecoratedInterval(interval(2.0,2.0), trv)) === after

    @test_broken MISSING_overlap(DecoratedInterval(interval(2.0,3.0), def), DecoratedInterval(interval(1.0,2.0), trv)) === metBy

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.5,2.5), com), DecoratedInterval(interval(1.0,2.0), com)) === overlappedBy

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,3.0), dac), DecoratedInterval(interval(1.0,2.0), def)) === startedBy

    @test_broken MISSING_overlap(DecoratedInterval(interval(1.0,3.0), com), DecoratedInterval(interval(1.0,1.0), dac)) === startedBy

    @test_broken MISSING_overlap(DecoratedInterval(interval(0.0,3.0), com), DecoratedInterval(interval(1.0,2.0), dac)) === contains

    @test_broken MISSING_overlap(DecoratedInterval(interval(0.0,3.0), com), DecoratedInterval(interval(2.0,2.0), def)) === contains

    @test_broken MISSING_overlap(DecoratedInterval(interval(0.0,2.0), def), DecoratedInterval(interval(1.0,2.0), trv)) === finishedBy

    @test_broken MISSING_overlap(DecoratedInterval(interval(0.0,2.0), dac), DecoratedInterval(interval(2.0,2.0), def)) === finishedBy

end


text:
import LegibleLambdas: LegibleLambda
export MathGate, mathgate

struct MathGate{N, F <: Union{LegibleLambda, Function}} <: PrimitiveBlock{N}
    f::F
end

function MathGate{N}(f::Union{LegibleLambda, Function}) where N
    return MathGate{N, typeof(f)}(f)
end

"""
    mathgate(nbits, f)

Create a [`MathGate`](@ref) with a math function `f` and number of bits.

    mathgate(f) -> f(n)

Lazy curried version of `mathgate`.

# Example

We can make a classical toffoli gate on quantum register.

```jldoctest; setup=:(using YaoBlocks, YaoArrayRegister, BitBasis)
julia> r = ArrayReg(bit"110")
ArrayReg{1, Complex{Float64}, Array...}
    active qubits: 3/3

julia> function toffli(b::BitStr)
           t = @inbounds b[1] ⊻ (b[3] & b[2])
           return @inbounds bit_literal(t, b[2], b[3])
       end
toffli (generic function with 1 method)

julia> g = mathgate(3, toffli)
┌ Warning: `MathGate` will be moved to `YaoExtensions.jl` in the next release.
└ @ YaoBlocks ~/.julia/packages/YaoBlocks/XXC9Y/src/primitive/math_gate.jl:8
mathgate(toffli; nbits=3)

julia> apply!(r, g) == ArrayReg(bit"111")
true

```
"""
mathgate(nbits::Int, f) = MathGate{nbits}(f)
mathgate(f::Union{LegibleLambda, Function}) = @λ(nbits->matgate(nbits, f))

function apply!(r::ArrayReg, m::MathGate{N, F}) where {N, F}
    nstate = zero(r.state)
    for b in basis(BitStr64{N})
        b2 = m.f(b)
        nstate[Int64(b2)+1, :] = view(r.state, Int64(b)+1, :)
    end
    r.state .= nstate
    return r
end

# TODO: use trait to correct this
function mat(::Type{T}, m::MathGate{N}) where {T, N}
    L = 1<<N
    vals = zeros(T, L)
    perm = zeros(Int, L)
    for b in basis(BitStr64{N})
        b2 = m.f(b)
        vals[Int64(b2)+1] += 1
        perm[Int64(b2)+1] = Int64(b) + 1
    end
    any(==(0), vals) && throw(ArgumentError("This `MathGate` is not unitary! Failed converting to a `PermMatrix`! maybe use `applymatrix` to check your block?"))
    return PermMatrix(perm, vals)
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 426b5c90-7c3f-11eb-3f2f-9dde1b32247f
using Luxor, Colors

# ╔═╡ d1c962e0-7c44-11eb-3d44-458c17db2497
using PlutoUI

# ╔═╡ 549df61e-7c3f-11eb-084e-c569635b421c
function barnsley_fern(n, hueshift=false)
	🐢 = Turtle()
	Penup(🐢)
	Pencolor(🐢, "green")

	p = Point(0, 0)
	for i in 1:n
		Reposition(🐢, Point((65p.x - 10), -(37p.y - 180)))
		Circle(🐢, 0.5)
		if hueshift
			HueShift(🐢, 0.0005)
		end
		prob = rand()
		if prob < 0.01
			p = Point(0, 0.16p.y)
		elseif prob < 0.86
			p = Point(0.85p.x + 0.04p.y, -0.04p.x + 0.85p.y + 1.6)
		elseif prob < 0.93
			p = Point(0.2p.x - 0.26p.y, 0.23p.x + 0.22p.y + 1.6)
		else
			p = Point(-0.15p.x + 0.28p.y, 0.26p.x + 0.24p.y + 0.44)
		end
	end
end

# ╔═╡ e523e4f0-7c44-11eb-3ec5-5fb18d2e5b46
@bind n Slider(1:1000:200000)

# ╔═╡ 624affb4-7c3f-11eb-382c-ff02831d8a5f
@png barnsley_fern(n) 400 400 "barnsley_fern.png"

# ╔═╡ 03178496-7c48-11eb-17a2-0933ef395323
@png barnsley_fern(n, true) 400 400 "barnsley_fern_he.png"

# ╔═╡ b3f0fd0a-7c3f-11eb-1639-31ab166f1b37
rand()

# ╔═╡ 140610ae-7c40-11eb-1763-832cd66375f7
fieldnames(Point)

# ╔═╡ Cell order:
# ╠═426b5c90-7c3f-11eb-3f2f-9dde1b32247f
# ╠═d1c962e0-7c44-11eb-3d44-458c17db2497
# ╠═549df61e-7c3f-11eb-084e-c569635b421c
# ╠═624affb4-7c3f-11eb-382c-ff02831d8a5f
# ╠═03178496-7c48-11eb-17a2-0933ef395323
# ╠═e523e4f0-7c44-11eb-3ec5-5fb18d2e5b46
# ╠═b3f0fd0a-7c3f-11eb-1639-31ab166f1b37
# ╠═140610ae-7c40-11eb-1763-832cd66375f7


text:
"""
    Aqua.test_stale_deps(package; [ignore])

Test that `package` loads all dependencies listed in `Project.toml`.

!!! note "Known bug"

    Currently, `Aqua.test_stale_deps` does not detect stale
    dependencies when they are stdlib.  This is considered a bug and
    may be fixed in the future.  Such a release is considered
    non-breaking.

# Arguments
- `packages`: a top-level `Module`, a `Base.PkgId`, or a collection of
  them.

# Keyword Arguments
- `ignore::Vector{Symbol}`: names of dependent packages to be ignored.
"""
test_stale_deps
function test_stale_deps(packages; kwargs...)
    @testset "$(result.label)" for result in analyze_stale_deps(packages, kwargs)
        @debug result.label result
        @test result ⊜ true
    end
end

analyze_stale_deps(packages; kwargs...) = analyze_stale_deps(packages, kwargs)
analyze_stale_deps(packages, kwargs) =
    [_analyze_stale_deps_1(pkg; kwargs...) for pkg in aspkgids(packages)]

function _analyze_stale_deps_1(pkg::PkgId; ignore::AbstractArray{Symbol} = Symbol[])
    label = "$pkg"

    result = root_project_or_failed_lazytest(pkg)
    result isa LazyTestResult && return result
    root_project_path = result

    @debug "Parsing `$root_project_path`"
    prj = TOML.parsefile(root_project_path)
    raw_deps = get(prj, "deps", nothing)
    if raw_deps === nothing
        return LazyTestResult(label, "No `deps` table in `$root_project_path`", true)
    end
    deps = [PkgId(UUID(v), k) for (k, v) in raw_deps]

    code = """
    $(Base.load_path_setup_code())
    Base.require($(reprpkgid(pkg)))
    for pkg in keys(Base.loaded_modules)
        pkg.uuid === nothing || println(pkg.uuid)
    end
    """
    cmd = Base.julia_cmd()
    output = read(`$cmd --startup-file=no --color=no -e $code`, String)
    @debug("Checked modules loaded in a separate process.", cmd, Text(code), Text(output))
    loaded_uuids = map(UUID, eachline(IOBuffer(output)))

    return _analyze_stale_deps_2(;
        pkg = pkg,
        deps = deps,
        loaded_uuids = loaded_uuids,
        ignore = ignore,
    )
end

# Side-effect -free part of stale dependency analysis.
function _analyze_stale_deps_2(;
    pkg::PkgId,
    deps::AbstractArray{PkgId},
    loaded_uuids::AbstractArray{UUID},
    ignore::AbstractArray{Symbol},
)
    label = "$pkg"
    deps_uuids = [p.uuid for p in deps]
    pkgid_from_uuid = Dict(p.uuid => p for p in deps)

    stale_uuids = setdiff(deps_uuids, loaded_uuids)
    stale_pkgs = [pkgid_from_uuid[uuid] for uuid in stale_uuids]
    stale_pkgs = [p for p in stale_pkgs if !(Symbol(p.name) in ignore)]

    if isempty(stale_pkgs)
        return LazyTestResult(
            label,
            """
            All packages in `deps` are loaded via `using $(pkg.name)`.
            """,
            true,
        )
    end

    stale_msg = join(("* $p" for p in stale_pkgs), "\n")
    msglines = [
        "Some package(s) in `deps` of $pkg are not loaded during via" *
        " `using $(pkg.name)`.",
        stale_msg,
        "",
        "To ignore from stale dependency detection, pass the package name to" *
        " `ignore` keyword argument of `Aqua.test_stale_deps`",
    ]
    return LazyTestResult(label, join(msglines, "\n"), false)
end


text:
@time @testset "Systems.NonDeterministicInput" begin include("unit_NonDeterministicInput.jl") end
@time @testset "Systems.System" begin include("unit_System.jl") end


text:
using ApproxFun, Base.Test

import ApproxFun: expα, expβ, expγ

#
# This tests the numerical stability of the formulæ appearing in Eq. (2.5) of:
#
# A.-K. Kassam and L. N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput., 26:1214--1233, 2005.
#

x = logspace(-15,2,18)

@test norm(expα(x)./expα(big(x))-1,Inf) < eps()
@test norm(expβ(x)./expβ(big(x))-1,Inf) < eps()
@test norm(expγ(x)./expγ(big(x))-1,Inf) < eps()


text:
function PI{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_pi), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
PI(::Type{ArbFloat}) = PI(ArbFloat{precision(ArbFloat)})

function SQRTPI{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_sqrt_pi), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
SQRTPI(::Type{ArbFloat}) = SQRTPI(ArbFloat{precision(ArbFloat)})

function LOG2{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_log2), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
LOG2(::Type{ArbFloat}) = LOG2(ArbFloat{precision(ArbFloat)})

function LOG10{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_log10), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
LOG10(::Type{ArbFloat}) = LOG10(ArbFloat{precision(ArbFloat)})

function EXP1{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_e), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
EXP1(::Type{ArbFloat}) = EXP1(ArbFloat{precision(ArbFloat)})

function EULER{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_euler), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
EULER(::Type{ArbFloat}) = EULER(ArbFloat{precision(ArbFloat)})

function CATALAN{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_catalan), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
CATALAN(::Type{ArbFloat}) = CATALAN(ArbFloat{precision(ArbFloat)})

function KINCHIN{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_kinchin), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
KHINCHIN(::Type{ArbFloat}) = KINCHIN(ArbFloat{precision(ArbFloat)})

function GLAISHER{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_apery), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
GLAISHER(::Type{ArbFloat}) = GLAISHER(ArbFloat{precision(ArbFloat)})

function APERY{P}(::Type{ArbFloat{P}})
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_apery), Void, (Ptr{ArbFloat}, Int), &z, P)
    z
end
APERY(::Type{ArbFloat}) = APERY(ArbFloat{precision(ArbFloat)})


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type EffSampleSizeAlgorithm

Abstract type for integrated autocorrelation length estimation algorithms.
"""
abstract type EffSampleSizeAlgorithm end
export EffSampleSizeAlgorithm

"""
    bat_eff_sample_size(
        v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
        [algorithm::EffSampleSizeAlgorithm]
    )

    bat_eff_sample_size(
        smpls::DensitySampleVector,
        [algorithm::EffSampleSizeAlgorithm]
    )

Estimate effective sample size estimation for variate series `v`, resp.
density samples `smpls`, separately for each degree of freedom.

Returns a NamedTuple of the shape

```julia
(result = eff_sample_size, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_eff_sample_size`, add methods to
    `bat_eff_sample_size_impl` instead.
"""
function bat_eff_sample_size end
export bat_eff_sample_size

function bat_eff_sample_size_impl end

function bat_eff_sample_size(
    v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
    algorithm = bat_default_withdebug(bat_eff_sample_size, Val(:algorithm), v)
)
    r = bat_eff_sample_size_impl(v, algorithm)
    result_with_args(r, (algorithm = algorithm,))
end

function bat_eff_sample_size(
    smpls::DensitySampleVector,
    algorithm = bat_default_withdebug(bat_eff_sample_size, Val(:algorithm), smpls);
    kwargs...
)
    r = bat_eff_sample_size_impl(smpls, algorithm; kwargs...)
    result_with_args(r, (algorithm = algorithm,))
end

function argchoice_msg(::typeof(bat_eff_sample_size), ::Val{:algorithm}, x::EffSampleSizeAlgorithm)
    "Using integrated autocorrelation length estimator $x"
end


text:
export GF, normalmatrix, allzeros

category_trait(::Type{<:GaloisField}) = FieldTrait

"""
    GF(p, r; mod=nothing, nr=0)

Create a class of element type `GaloisField` of order `p^r`.
`p` must be a prime integer and `r` a positive integer.
If `nr != 0` is given, it triggers the use of an alternate modulus
for the `GFImpl` class.

Optionally either the modulus polynomial `mod` or a integer search modifier `nr` may be given
to control the selection of the modulus polynomial.
"""
function GF(p::Integer, r::Integer; nr::Integer=0, mod=nothing)
    Q = GFImpl(p, r, nr=nr, mod=mod)
    ord = order(Q)
    r = dimension(Q)
    r == 1 && mod === nothing && return Q

    T = mintype_for(p, r, true)
    gen = first(Iterators.filter(x -> order(x) == ord-1, Q))

    function gfclass(p, ord, gen)
        c = fill(gen, ord)
        c[1] = c[2] = 1
        cumprod!(c, c)
        c[1] = 0
        exptable = [T(tonumber(x, p)) for x in c]
        logtable = invperm(exptable .+ 1) .- 1
        zechtable = logtable[[T(tonumber(x + 1, p)) for x in c] .+ 1]
        exptable, zechtable
    end

    Id = (p, r, ord, tonumber(gen, p))

    new_class(gfclass, GaloisField{Id,T,Q}, p, ord, gen)
end
function GF(n::Integer; mod=nothing, nr=0)
    f = factor(n)
    length(f) == 1 || throw(ArgumentError("$n is not p^r with p prime and r >= 1"))
    p, r = f.pe[1]
    GF(p, r, mod=mod, nr=nr)
end

"""
    GaloisField{Id,T,Q}[num::Integer]

Ring element constructor. `num` is *not* the canonical homomorphism, but enumerates
all elements of `GF(p, r)` in `0:p^r-1`.
The numbers `G[0:p-1]` correspond to the base field, and `G[p]` to the polynomial `x` in
the representation of `Q`.
"""
function Base.getindex(::Type{G}, num::Integer) where G<:GaloisField
    tv = log_calc(num, G) + 1
    G(tv, NOCHECK)
end
function Base.getindex(::Type{G}, ind::AbstractRange{<:Integer}) where G<:GaloisField
    p, r, ord = characteristic(G), dimension(G), order(G)
    zt = gettypevar(G).zechtable
    logt = logg0(G)
    v = Vector{G}(undef, length(ind))
    for (i, k) in enumerate(ind)
        tv = log_calc(tocoeffs(k, p, r, ord, Int), ord, zt, logt) + 1
        v[i] = G(tv, NOCHECK)
    end
    v
end

Base.collect(::Type{G}) where G<:GaloisField = G[0:order(G)-1]

function GaloisField{Id,T,Q}(a::GaloisField{Id,T,Q}) where {Id,T,Q}
    GaloisField{Id,T,Q}(a.val, NOCHECK)
end
function GaloisField{Id,T,Q}(a::GaloisField{Id2,T,Q2}) where {Id,T,Q,Id2,Q2}
    GaloisField{Id,T,Q}(a.val, NOCHECK)
end
(::Type{G})(q::Q) where {Id,T,Q,G<:GaloisField{Id,T,Q}} = G[tonumber(q, characteristic(Q))]
(::Type{G})(q::Q) where {Id,T,Q,G<:GaloisField{Id,T,<:Quotient{<:UnivariatePolynomial{Q}}}} = G[q.val]

Quotient(g::G) where {Id,T,Q,G<:GaloisField{Id,T,Q}} = quotient(Q, g)
Quotient(::Type{G}) where {Id,T,Q,G<:GaloisField{Id,T,Q}} = Q
Polynomial(::Type{G}) where G<:GaloisField = Polynomial(Quotient(G))

promote_rule(G::Type{GaloisField{Id,T,Q}}, ::Type{<:Integer}) where {Id,T,Q} = G
_promote_rule(G::Type{GaloisField{Id,T,Q}}, ::Type{Q}) where {Id,T,Q} = G
_promote_rule(G::Type{<:GaloisField{Id,T,<:Quotient{<:UnivariatePolynomial{Q}}}}, ::Type{Q}) where {Id,T,Q} = G

function quotient(::Type{Q}, g::G) where {Id,T,Z,Q<:Quotient{UnivariatePolynomial{Z,:α}},G<:GaloisField{Id,T,Q}}
    et = gettypevar(G).exptable
    toquotient(et[g.val+1], Q)
end

(::Type{Q})(g::G) where {Id,T,Z,Q<:Quotient{UnivariatePolynomial{Z,:α}},G<:GaloisField{Id,T,Q}} = quotient(Q, g)

import Base.Broadcast: broadcastable
broadcastable(x::Type{<:GaloisField}) = collect(x)

function isless(a::G, b::G) where G<:GaloisField
    isless(Quotient(a), Quotient(b))
end

basetype(::Type{GaloisField{Id,T,Q}}) where {Id,T,Q} = Q
depth(G::Type{<:GaloisField}) = depth(basetype(G)) + 1
characteristic(G::Type{<:GaloisField}) = characteristic(basetype(G))
dimension(G::Type{<:GaloisField}) = dimension(basetype(G))
order(G::Type{<:GaloisField}) = order(basetype(G))
lognegone(G::Type{<:GaloisField}) = characteristic(G) == 2 ? 0 : (order(G) - 1) ÷ 2
modulus(G::Type{<:GaloisField}) = modulus(basetype(G))
issimpler(a::G, b::G) where G<:GaloisField = a.val < b.val

# multiplication using lookup tables
function *(a::G, b::G) where G<:GaloisField
    iszero(a.val) && return a
    iszero(b.val) && return b
    ord = order(G)
    # beware of overflows
    G(mod(a.val - 2 + b.val, ord - 1) + 1, NOCHECK)
end
function /(a::G, b::G) where G<:GaloisField
    iszero(a.val) && return a
    a.val == b.val && return one(G)
    iszero(b.val) && division_error()
    ord = order(G)
    # beware of overflows
    G(mod(a.val - 1 + ord - b.val, ord - 1) + 1, NOCHECK)
end
function inv(a::G) where G<:GaloisField
    ord = order(G)
    iszero(a) && division_error()
    nlog = mod(ord - a.val, ord - 1) + 1
    G(nlog, NOCHECK)
end

"""
    num_primitives(::Type{G})

Number of primitive elements (generators of multiplicative group) of GaloisField `G`.
"""
function num_primitives(::Type{G}) where G<:Union{ZZmod,GaloisField,QuotientRing}
    Primes.totient(order(G)-1)
end

"""
    isprimitive(g::G)

Return iff `g` is a primitive element of its ring (its powers are the complete multilicative subgroup of `G`)
"""

function isprimitive(g::G) where G<:Union{ZZmod,GaloisField,QuotientRing}
    iszero(g) && return false
    n = order(G) - 1
    for p in keys(factor(n))
        isone(g ^ (n ÷ p)) && return false
    end
    true
end

import Base: ^, log

function ^(a::G, x::Integer) where G<:GaloisField
    ord = order(G)
    if iszero(a)
        return x > 0 ? a : division_error()
    end
    T = typeof(a.val)
    v = (log(a))
    #xx = T(mod(x, ord - 1))
    nlog = T(mod(widemul(v, x), ord-1) + 1)
    G(nlog, NOCHECK)
end

"""
    log(g::GaloisField)

Return the integer `k in 0:order(G)-2` with `g == α ^ k`, where `α` is the generator of `G`
or `-1` if `g == 0`. For example `log(one(G) == 0`.
"""
log(a::G) where G<:GaloisField = a.val - 1

"""
    log_zech(k::Integer, G::Type{<:GaloisField})

If `α` is the generator of `G`, for `k >= 0` return the `log(α^k + 1)`, for `k < 0` return `0`.
In other words: `α ^ log_zech(k, G) == α ^ k + 1` for `k >= 0`.
"""
function log_zech(k, G::Type{<:GaloisField})
    ord = order(G)
    zt = gettypevar(G).zechtable
    log_zech(k, ord, zt)
end
function log_zech(k::Integer, ord::Integer, zt::AbstractVector)
    k < 0 && return 0
    zt[mod(k, ord - 1) + 2] - 1
end

"""
    generator(::Type{<:GaloisField})

Return the generator element of this implementation of Galois field.
"""
function generator(G::Type{<:GaloisField})
    G[generatornum(G)]
end

"""
   generatornum(::Type{<:GaloisField})

Return numeric representation of generator (in range `0:order(G)-1`). Ideally `== characteristic(G)`.
"""
generatornum(::Type{<:GaloisField{Id}}) where Id = Id[4]

division_error() = throw(ArgumentError("cannot invert zero"))

*(a::G, b::Integer) where G<:GaloisField = G[mod(b, characteristic(G))] * a
*(b::Integer, a::G) where G<:GaloisField =  a * b
==(a::G, b::G) where G<:GaloisField = a.val == b.val

function typedep(::Type{G}) where G<:GaloisField
    ord = order(G)
    zechtable = gettypevar(G).zechtable
    e = lognegone(G)
    ord, e, zechtable
end

function +(a::G, b::G) where G <:GaloisField
    iszero(a) && return b
    iszero(b) && return a
    a = a.val
    b = b.val
    if a < b
        a, b = b, a
    end
    ord, e, zechtable = typedep(G)
    k = a - b
    k == e && return zero(G)
    @inbounds G(mod1(zechtable[k+2] - 1 + b, ord - 1), NOCHECK)
end

function -(a::G, b::G) where G <:GaloisField
    iszero(b) && return a
    iszero(a) && return -b
    a = a.val
    b = b.val
    a == b && return zero(G)
    ord, e, zechtable = typedep(G)
    b = b > e ? b - e : b + e
    if a < b
        a, b = b, a
    end
    k = a - b
    G(mod1(zechtable[k+2] - 1 + b, ord - 1), NOCHECK)
end

function -(a::G) where G<:GaloisField
    characteristic(G) == 2 && return a
    iszero(a) && return a
    e = lognegone(G)
    a = a.val
    a = a > e ? a - e : a + e
    G(a, NOCHECK)
end

iszero(a::GaloisField) = iszero(a.val)
isone(a::GaloisField) = isone(a.val)
isunit(a::GaloisField) = !iszero(a)
issimple(a::GaloisField) = true
value(g::GaloisField) = value(toquotient(g))

zero(::Type{G}) where G<:GaloisField = G(0, NOCHECK)
one(::Type{G}) where G<:GaloisField = G(1, NOCHECK)

function rand(r::AbstractRNG, ::SamplerType{G}) where G<:GaloisField
    ord = order(G)
    G[rand(r, 0:ord-1)]
end

function Base.show(io::IO, g::Type{<:GaloisField})
    sc(f, g) = try f(g) catch; "?" end
    print(io, :GaloisField, '{', sc(characteristic, g), ',', sc(dimension, g), '}')
end

function tonumber(a::Quotient, p::Integer)
    s = 0
    for c = reverse(a.val.coeff)
        s = s * p + c.val
    end
    s
end

function toquotient(g::G) where {Id,T,Q,G<:GaloisField{Id,T,Q}}
    tvar = gettypevar(G)
    exptable = tvar.exptable
    toquotient(exptable[g.val+1], Q)
end

function toquotient(a::Integer, ::Type{Q}) where {Z,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P}}
    p = characteristic(Q)
    r = dimension(Q)
    ord = order(Q)
    Q(tocoeffs(a, p, r, ord, Z))
end
function tocoeffs(a::Integer, p::Integer, r::Integer, ord::Integer, Z::Type)
    c = zeros(Z, r)
    0 <= a < ord || throw(ArgumentError("index must be in 0:$(ord-1)"))
    b = a
    for i = 1:r
        iszero(b) && break
        b, c[i] = divrem(b, p)
    end
    c
end

"""
    GFImpl(p[, m=1]; mod=modulus, nr=0)

Return a representation for Galois Field `GFImpl(p, m)`. `p` must be a prime number and
`m` a positive integer.

If `m == 1`, `ZZmod{p}` is returned,
otherwise an irreducible polynomial `g ∈ ZZmod{p}[:x] and deg(g) == m` is generated and
`ZZmod{p}[:x]/(g)` is returned.

Optionally either the modulus polynomial `mod` or a integer search modifier `nr` may be given
to control the selection of the modulus polynomial.

Elements of the field can be created like
```´
    GF7 = GFImpl(7)
    GF7(5)  # or
    GF53 = GFImpl(5, 3)
    GF53([1,2,3])
```
"""
function GFImpl(p::Integer, m::Integer=1; nr::Integer=0, mod=nothing)
    isprime(p) || throw(ArgumentError("base $p must be prime"))
    m > 0 || throw(ArgumentError("exponent m=$m must be positive"))
    Z = ZZ / p
    m == 1 && mod === nothing && return Z
    mm = intpower(p, m) - 1
    list = factors(mm)
    P = Z[:α]
    if mod === nothing
        x = monom(P)
        nx = max(nr, 0)
        for gen = irreducibles(P, m)
            if order(x, gen, list) == mm
                nx == 0 && return P / gen
                nx -= 1
            end
        end
        throw(ArgumentError("no irreducible polynomial of degree $m found with generator p(γ) = γ (nr = $nr)"))
    else
        gen = P(Z.(mod.coeff))
        if isirreducible(gen)
            return P / gen
        end
        throw(ArgumentError("given polynomial $gen is not irreducible in $P"))
    end
end

function Base.show(io::IO, g::G) where G<:GaloisField

    m = dimension(G)
    p = characteristic(G)
    c = toquotient(g).val.coeff
    n = length(c)
    cc(i) = i > n ? 0 : c[i].val
    print(io, '{', cc(m))
    for k = m-1:-1:1
        print(io, ':', cc(k))
    end
    print(io, '%', p, '}')
end

"""
    normalmatrix(a::Q[, m])

Return a matrix of element type `ZZ/p`, whose colums are the coefficient
vectors of `a^(p^i) for i = 0:m-1`.

Here Q is a galois field of characteristic `p` and order `p^r`.

`m` m defaults to `r`.

If `normalmatrix(a))` is regular, the field elements `a^(p^i) for i = 0:r-1` form a
base of `Q` as a vector space over `ZZ/p` (a normal base).
"""
function normalmatrix(a::Q, m::Integer=0) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P}}
    p = characteristic(Q)
    r = dimension(Q)
    m = m <= 0 ? r : m
    M = Matrix{Z}(undef, r, m)
    for i = 0:m-1
        c = a.val.coeff
        k = length(c)
        for j = 1:r
            M[j,i+1] = j <= k ? c[j] : 0
        end
        a ^= p
    end
    M
end

"""
    normalmatrix(::Type{Q}[, m])

Return `normalmatrix(a, m)` for the first `a` in `Q` for which this ihas maximal rank.
"""
function normalmatrix(::Type{Q}, m::Integer=0) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P}}
    normalmatrix(normalbase(Q), m)
end

function normalbases(::Type{Q}) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P}}
    r = dimension(Q)
    Base.Iterators.filter(x->rank(normalmatrix(x, r)) == r, Q)
end
"""
    normalbase(::Type{Q})

Find the first `a in Q` for which `normalmatrix(a)` is regular.
"""
function normalbase(::Type{Q}) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P}}
    bases = normalbases(Q)
    isempty(bases) && throw(ArgumentError("quotient type with modulus $(modulus(Q)) has no normal bases - probably modulus is not an irreducible polynomial"))
    first(bases)
end

import Base: *

"""
    sized(a::Vector, r)

Return a vector of length `r`, which starts with `a` and is filled up with zeros if required.
"""
function sized(a::Vector{Z}, r::Integer) where Z
    n = length(a)
    n == r ? a : n < r ? vcat(a, zeros(Z, r - n)) : a[1:r]
end

mulsized(M::AbstractMatrix{Z}, a::Vector{Z}) where Z<:Ring = M * sized(a, size(M, 2))

function *(M::AbstractMatrix{Z}, a::Q) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P}}
    mulsized(M, a.val.coeff)
end

function monom(::Type{Q}, k::Integer) where {P<:UnivariatePolynomial,Q<:Quotient{P}}
    Q(monom(P, k))
end

"""
    homomorphism(Q, R [,nr=0])

Return a function `iso: Q -> R`, which describes an homomorphism between two Galois fields
`Q` and `R` of the same characteristic. If `Q == R` that are the Frobenius automorphisms,
if `order(Q) == order(R)` isomorphisms, in the case of `order(R) == order(S)^s` with s > 1
the (injective) monomorphisms.

The optional `nr ∈ 0:r-1` produces all possible monomorphisms (automorphisms) between `Q` and `R`.
In the automorphism case, `nr = 0` is the identity.co
"""
function homomorphism end
function homomorphism(f::Function, ::Type{G}, ::Type{H}) where {G,H}
    Hom{G,H}(f)
end
function _homomorphism(::Type{Q}, ::Type{R}) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P},R<:Quotient{P}}

    r = dimension(Q)
    s = dimension(R)
    p = characteristic(Q)
    pr = p^r
    mod(s, r) == 0 || throw(ArgumentError("dimension of Q ($r) must divide that of R ($s)"))
    f = normalbase(Q)
    M = normalmatrix(f, r)
    M1 = inv(M)
    k = 3 - (p % 2)
    h = M1 * f^k

    xr = monom(R, 1)
    L = ((xr^k)^pr for k = 0:s-1)
    S = hcat(collect(sized(x.val.coeff, s) for x in L)...)
    for i = 1:s
        S[i,i] -= one(Z)
    end
    K = Matrix(nullspace(S))

    for g0 in Q
        if R == Q
            g = f
            N = M
        else
            g = R(K * g0)
            deg(g.val) <= 0 && continue
            N = normalmatrix(g, r)
        end
        if sized((g^k).val.coeff, s) == N * h
            #if rank(N) != r
            #    throw(ErrorException("expected rank $r, but was $(rank(N)) $g"))
            #end
            return N, M1
        end
    end
    throw(ErrorException("no homomorphism found - not reachable"))
end

function _homomorphism(::Type{Z}, ::Type{R}) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},R<:Quotient{P}}
    1, 1
end

function homomorphism(iso::Function, nr::Integer=0)
    N = iso.A.N
    M1 = iso.A.M1
    Q = iso.A.Q
    R = iso.A.R
    _homomorphism(Q, R, N, M1, nr)
end

function homomorphism(::Type{Z}, ::Type{H}, nr::Integer=0) where {Z<:ZZmod,H<:GaloisField}
    Hom{Z,H}(x -> H(x))
end

function homomorphism(::Type{G}, ::Type{H}, nr::Integer=0) where {G<:GaloisField,H<:GaloisField}
    N, M1 = _homomorphism(basetype(G), basetype(H))
    Hom{G,H}(_homomorphism(G, H, N, M1, nr))
end

function homomorphism(::Type{Q}, ::Type{R}, nr::Integer=0) where {Z<:ZZmod,P<:UnivariatePolynomial{Z,:α},Q,R<:Quotient{P}}
    N, M1 = _homomorphism(Q, R)
    _homomorphism(Q, R, N, M1, nr)
end

function _homomorphism(::Type{Q}, ::Type{R}, N, M1, nr::Integer) where {Q,R}
    r = size(N, 2)
    nr = mod(nr, r)
    # cyclic permutation of columns of N
    if nr != 0
        N = hcat(N[:,nr+1:r], N[:,1:nr])
    end
    A = (T = N * M1, N = N, M1 = M1, Q = Q, R = R)
    quot(x) = x
    quot(x::GaloisField) = Quotient(x)
    iso(a::Q) = R(A.T * quot(a))
    iso
end

"""
    allzeros(p, vx)

Assume `p` is an irreducible polynomial over `ZZ/q`.
If `p(vx) == 0` for a galois field element `vx`,
find all zeros of `p` in the galois field, vx belongs to.
"""
function allzeros(p::P, vx::Q) where {P<:UnivariatePolynomial,Q<:Ring}
    q = characteristic(Q)
    r = deg(p)
    return (vx^q^k for k = 0:r-1)
    #=
    r = deg(p)
    m = dimension(Q)
    M = normalmatrix(normalbase(P/p), r)
    N = hcat( (sized((vx^k).val.coeff, m) for k = 0:r-1)...) *  M
    a = inv(M)[:,2]
    cp(N, k) = [N[:,k+1:r] N[:,1:k]] # cyclically permutating columns
    [ Q(cp(N, k) * a) for k in 0:r-1]
    =#
end

"""
    order(p, q, xlist)

Return first `n` in `xlist` where `rem(p^n, q) == 1`; otherwise 0.
"""
function order(p, q, xlist)
    for x in xlist
        rem(p^x, q) == 1 && return x
    end
    0
end

function order(x::G) where G<:GaloisField
    iszero(x) && return 0
    isone(x) && return 1
    ord = mult_order(G)
    ord ÷ gcd(x.val - 1, ord)
end

function logg0(G)
    p = characteristic(G)
    ord = order(G)
    zt = gettypevar(G).zechtable
    logg0(p, ord, zt)
end
function logg0(p::Integer, ord::Integer, zt)
    v = Vector{Int}(undef, p)
    v[1] = -1
    s = v[2] = 0
    for i = 3:length(v)
        s = log_zech(s, ord, zt)
        v[i] = s
    end
    v
end

loggx(a::AbstractVector, G::Type{<:GaloisField}) = loggx(a, order(G), gettypevar(G).zechtable, logg0(G))
function loggx(a::AbstractVector, ord::Integer, zt::AbstractVector, logt::AbstractVector)
    loggi(k) = logt[k+1]
    k = findfirst(!iszero, a)
    k === nothing && return -1
    accu = loggi(a[k]) + k - 1
    for i = k+1:length(a)
        ai = a[i]
        if !iszero(ai)
            lai = loggi(ai) + i - 1
            accu = log_zech(mod(accu - lai, ord-1), ord, zt) + lai
        end
    end
    mod(accu, ord-1)
end

function log_calc(k::Integer, G::Type{<:GaloisField})
    p, r, ord = characteristic(G), dimension(G), order(G)
    zt = gettypevar(G).zechtable
    log_calc(tocoeffs(k, p, r, ord, Int), ord, zt, logg0(p, ord, zt))
end
function log_calc(a::AbstractVector, ord::Integer, zt::AbstractVector, logt::AbstractVector)
    loggx(a, ord, zt, logt)
end


text:
function get_ae(obj::WrappingUnnormalizedProvider)
    return jcall(obj, "getAe", jdouble, ())
end

function get_max_degree(obj::WrappingUnnormalizedProvider)
    return jcall(obj, "getMaxDegree", jint, ())
end

function get_max_order(obj::WrappingUnnormalizedProvider)
    return jcall(obj, "getMaxOrder", jint, ())
end

function get_mu(obj::WrappingUnnormalizedProvider)
    return jcall(obj, "getMu", jdouble, ())
end

function get_offset(obj::WrappingUnnormalizedProvider, arg0::AbsoluteDate)
    return jcall(obj, "getOffset", jdouble, (AbsoluteDate,), arg0)
end

function get_reference_date(obj::WrappingUnnormalizedProvider)
    return jcall(obj, "getReferenceDate", AbsoluteDate, ())
end

function get_tide_system(obj::WrappingUnnormalizedProvider)
    return jcall(obj, "getTideSystem", TideSystem, ())
end

function on_date(obj::WrappingUnnormalizedProvider, arg0::AbsoluteDate)
    return jcall(obj, "onDate", UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics, (AbsoluteDate,), arg0)
end


text:
using MLBase: gridtune, Reverse
using Compat: @compat

function modeltune(dataset::Persa.CFDatasetAbstract, estfun::Function, params::@compat(Tuple{AbstractString, Any})...; verbose = true, k = 0.9, seed = NaN)
  if isnan(seed)
    ds_train, ds_test = get(HoldOut(dataset, k))
  else
    ds_train, ds_test = get(HoldOut(dataset, k, seed))
  end

  estfun2(args...) = estfun(ds_train, args...)

  evalfun2(model) = aval(model, ds_test).mae

  return r = gridtune(estfun2, evalfun2,
              params...;
              ord=Reverse,    # smaller msd value indicates better model
              verbose=verbose)   # show progress information
end

max_stars_repo_name: JuliaRecsys/ModelTuneCF

max_issues_repo_name: JuliaRecsys/ModelTuneCF

max_forks_repo_name: JuliaRecsys/ModelTuneCF


text:
# the channel example using a state-based approach
#
using DiscreteEvents, Printf, Random

abstract type Q end  # states
struct Idle <: Q end
struct Busy <: Q end
abstract type Σ end  # events
struct Arrive <: Σ end
struct Leave <: Σ end

mutable struct Server
    id::Int64
    name::AbstractString
    input::Channel
    output::Channel
    op     # operation to take
    state::Q
    token  # current token

    Server(id, name, input, output, op) = new(id, name, input, output, op, Idle(), nothing)
end

arrive(A::Server) = event!(SF(δ, A, A.state, Arrive()), SF(isready, A.input))

function δ(A::Server, ::Idle, ::Arrive)
    A.token = take!(A.input)
    @printf("%5.2f: %s %d took token %d\n", tau(), A.name, A.id, A.token)
    A.state=Busy()
    event!(SF(δ, A, A.state, Leave()), after, rand())
end

function δ(A::Server, ::Busy, ::Leave)
    put!(A.output, A.op(A.id,A.token))
    A.state=Idle()
    arrive(A)
end

δ(A::Server, q::Q, σ::Σ) =               # fallback transition
        println(stderr, "$(A.name) $(A.id) undefined transition $q, $σ")

reset!(𝐶)
Random.seed!(123)

ch1 = Channel(32)  # create two channels
ch2 = Channel(32)

s = shuffle(1:8)
for i in 1:2:8
    arrive(Server(s[i], "foo", ch1, ch2, +))
    arrive(Server(s[i+1], "bar", ch2, ch1, *))
end

put!(ch1, 1) # put first token into channel 1

run!(𝐶, 10)


text:
function sdf_capsule_capsule_constraint(xa, qa, ha, ra, xb, qb, hb, rb)
    sdf_capsule_capsule(xa, qa, ha, ra, xb, qb, hb, rb)
end

sdf_capsule_capsule_constraint(p1, o1, h1, r1, p2, o2, h2, r2)

function sdf_capsule_capsule_constraint_jacobian(relative, xa, qa, ha, ra, xb, qb, hb, rb;
    attjac=false)

    # get closest points
    pa, pb = closest_points_capsule_capsule(xa, qa, ha, ra, xb, qb, hb, rb)
    d = pb - pa

    # Jacobians
    if relative == :parent
        dx = FiniteDiff.finite_difference_jacobian(x -> closest_points_capsule_capsule(x, qa, ha, ra, xb, qb, hb, rb)[2], xa)
        dx -= FiniteDiff.finite_difference_jacobian(x -> closest_points_capsule_capsule(x, qa, ha, ra, xb, qb, hb, rb)[1], xa)
        dq = FiniteDiff.finite_difference_jacobian(q -> closest_points_capsule_capsule(xa, Quaternion(q..., false), ha, ra, xb, qb, hb, rb)[2], vector(qa))
        dq -= FiniteDiff.finite_difference_jacobian(q -> closest_points_capsule_capsule(xa, Quaternion(q..., false), ha, ra, xb, qb, hb, rb)[1], vector(qa))
        attjac && (Q *= LVᵀmat(qa))
    elseif relative == :child
        dx = FiniteDiff.finite_difference_jacobian(x -> closest_points_capsule_capsule(xa, qa, ha, ra, x, qb, hb, rb)[2], xb)
        dx -= FiniteDiff.finite_difference_jacobian(x -> closest_points_capsule_capsule(xa, qa, ha, ra, x, qb, hb, rb)[1], xb)
        dq = FiniteDiff.finite_difference_jacobian(q -> closest_points_capsule_capsule(xa, qa, ha, ra, xb, Quaternion(q..., false), hb, rb)[2], vector(qb))
        dq -= FiniteDiff.finite_difference_jacobian(q -> closest_points_capsule_capsule(xa, qa, ha, ra, xb, Quaternion(q..., false), hb, rb)[1], vector(qb))
        attjac && (Q *= LVᵀmat(qb))
    end

    return d' ./ norm(d) * [dx dq]
end

Jp = constraint_jacobian(:parent, p1, o1, h1, r1, p2, o2, h2, r2, attjac=false)
Jc = constraint_jacobian(:child, p1, o1, h1, r1, p2, o2, h2, r2, attjac=false)

using FiniteDiff
Jpx = FiniteDiff.finite_difference_jacobian(x -> constraint(x, o1, h1, r1, p2, o2, h2, r2), p1)
Jpq = FiniteDiff.finite_difference_jacobian(q -> constraint(p1, Quaternion(q..., false), h1, r1, p2, o2, h2, r2), vector(o1))

Jcx = FiniteDiff.finite_difference_jacobian(x -> constraint(p1, o1, h1, r1, x, o2, h2, r2), p2)
Jcq = FiniteDiff.finite_difference_jacobian(q -> constraint(p1, o1, h1, r1, p2, Quaternion(q..., false), h2, r2), vector(o2))

norm(Jp - [Jpx Jpq])
norm(Jc - [Jcx Jcq])

function force_mapping(relative, xa, qa, ha, ra, xb, qb, hb, rb)
    pa, pb = closest_points_capsule_capsule(xa, qa, ha, ra, xb, qb, hb, rb)

    na = (pb - pa) ./ norm(pb - pa)
    if norm(na - [1.0; 0.0; 0.0]) > 1.0e-8
        ta1 = RotX(0.5 * π) * na # add edge case
    else
        ta1 = RotY(0.5 * π) * na
    end
    ta2 = cross(na, ta1)

    if relative == :parent
        return [ta1 ta2 na]
    elseif relative == :child
        return -1.0 * [ta1 ta2 na]
    end
end

function force_mapping_jacobian(relative, jacobian, xa, qa, ha, ra, xb, qb, hb, rb, λ;
    attjac=false)
    if jacobian == :parent
        dx = FiniteDiff.finite_difference_jacobian(x -> force_mapping(relative, x, qa, ha, ra, xb, qb, hb, rb) * λ, xa)
        dq = FiniteDiff.finite_difference_jacobian(q -> force_mapping(relative, xa, Quaternion(q..., false), ha, ra, xb, qb, hb, rb) * λ, vector(qa))
        attjac && (dq *= LVᵀmat(qa))
    elseif jacobian == :child
        dx = FiniteDiff.finite_difference_jacobian(x -> force_mapping(relative, xa, qa, ha, ra, x, qb, hb, rb) * λ, xb)
        dq = FiniteDiff.finite_difference_jacobian(q -> force_mapping(relative, xa, qa, ha, ra, xb, Quaternion(q..., false), hb, rb) * λ, vector(qb))
        attjac && (dq *= LVᵀmat(qa))
    end
    return [dx dq]
end

force_mapping(:parent, p1, o1, h1, r1, p2, o2, h2, r2)
force_mapping(:child, p1, o1, h1, r1, p2, o2, h2, r2)

force_mapping_jacobian(:parent, :parent, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)
force_mapping_jacobian(:parent, :child, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)
force_mapping_jacobian(:child, :parent, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)
force_mapping_jacobian(:child, :child, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)

function torque_mapping(relative, xa, qa, ha, ra, xb, qb, hb, rb)
    pa, pb = closest_points_capsule_capsule(xa, qa, ha, ra, xb, qb, hb, rb)
    F = force_mapping(relative, xa, qa, ha, ra, xb, qb, hb, rb)

    if relative == :parent
        ra = pa - xa
        τa = rotation_matrix(inv(qa)) * skew(ra) * F
        return τa
    elseif relative == :child
        rb = pb - xb
        τb = rotation_matrix(inv(qb)) * skew(rb) * F
        return τb
    end
end

function torque_mapping_jacobian(relative, jacobian, xa, qa, ha, ra, xb, qb, hb, rb, λ;
    attjac=false)
    if jacobian == :parent
        dx = FiniteDiff.finite_difference_jacobian(x -> torque_mapping(relative, x, qa, ha, ra, xb, qb, hb, rb) * λ, xa)
        dq = FiniteDiff.finite_difference_jacobian(q -> torque_mapping(relative, xa, Quaternion(q..., false), ha, ra, xb, qb, hb, rb) * λ, vector(qa))
        attjac && (dq *= LVᵀmat(qa))
    elseif jacobian == :child
        dx = FiniteDiff.finite_difference_jacobian(x -> torque_mapping(relative, xa, qa, ha, ra, x, qb, hb, rb) * λ, xb)
        dq = FiniteDiff.finite_difference_jacobian(q -> torque_mapping(relative, xa, qa, ha, ra, xb, Quaternion(q..., false), hb, rb) * λ, vector(qb))
        attjac && (dq *= LVᵀmat(qa))
    end
    return [dx dq]
end

torque_mapping(:parent, p1, o1, h1, r1, p2, o2, h2, r2)
torque_mapping(:child, p1, o1, h1, r1, p2, o2, h2, r2)

torque_mapping_jacobian(:parent, :parent, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)
torque_mapping_jacobian(:parent, :child, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)
torque_mapping_jacobian(:child, :parent, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)
torque_mapping_jacobian(:child, :child, p1, o1, h1, r1, p2, o2, h2, r2, ones(3), attjac=false)

function impulse_map(relative, xa, qa, ha, ra, xb, qb, hb, rb)
    [
        force_mapping(relative, xa, qa, ha, ra, xb, qb, hb, rb);
        torque_mapping(relative, xa, qa, ha, ra, xb, qb, hb, rb);
    ]
end

function impulse_map_jacobian_configuration(relative, jacobian, xa, qa, ha, ra, xb, qb, hb, rb, λ)
    [
        force_mapping_jacobian(relative, jacobian, xa, qa, ha, ra, xb, qb, hb, rb, λ);
        torque_mapping_jacobian(relative, jacobian, xa, qa, ha, ra, xb, qb, hb, rb, λ);
    ]
end

function tangential_velocity(x3a, q3a, x2a, q2a, ha, ra, x3b, q3b, x2b, q2b, hb, rb, h)
    # contact point kinematics
    ka, kb = kinematics_capsule_capsule(x3a, q3a, ha, ra, x3b, q3b, hb, rb)
    # capsule a points
    p3a = x3a + vector_rotate(ka, q3a)
    p3b = x3b + vector_rotate(kb, q3b)
    # capsule b points
    p2a = x2a + vector_rotate(ka, q2a)
    p2b = x2b + vector_rotate(kb, q2b)
    # velocities
    va = (p3a - p2a) ./ h
    vb = (p3b - p2b) ./ h
    # relative velocity
    Δv = vb - va
    # projection matrix
    na = (p3b - p3a) ./ norm(p3b - p3a)
    if norm(na - [1.0; 0.0; 0.0]) > 1.0e-8
        ta1 = RotX(0.5 * π) * na # add edge case
    else
        ta1 = RotY(0.5 * π) * na
    end
    ta2 = cross(na, ta1)
    P = [ta1'; ta2']
    # tangential velocity
    return P * Δv
end

h = 0.1
x_shift = [0.0; 0.0; 1.0]
tangential_velocity(p1, o1, p1 + x_shift * h, o1, h1, r1, p2, o2, p2, o2, h2, r2, h)

## constraints
function constraint(model, s::AbstractVector{T}, γ::AbstractVector{T},
    x3a::AbstractVector{T}, q3a::Quaternion{T}, x2a::AbstractVector{T}, q2a::Quaternion{T},
    x3b::AbstractVector{T}, q3b::Quaternion{T}, x2b::AbstractVector{T}, q2b::Quaternion{T},
    h::T) where T

    # transforms the velocities of the origin of the link into velocities
    vp = tangential_velocity(x3a, q3a, x2a, q2a, h1, r1, x3b, q3b, x2b, q2b, h2, r2, h)

    SVector{4,T}(
        sdf_capsule_capsule(x3a, q3a, ha, ra, x3b, q3b, hb, rb),
        model.friction_coefficient * γ[1] - γ[2],
        (vp - s[@SVector [3,4]])...)
end


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

using ColorTypes, ImageView

export view

"""
    view(mri::MRI, plane::Char='a')

View an `MRI` structure in a slice viewer, along the specified `plane` ('a':
axial; 's': sagittal; 'c': coronal).
"""
function view(mri::MRI, plane::Char='a')

  # Find which axes of the volume correspond to the specified viewing plane
  ax = view_axes(mri.vox2ras, plane)

  (ax1, ax2) = abs.(ax)
  (flip1, flip2) = (ax .< 0)

  if mri.ispermuted
    (ax1 == 1) && (ax1 = 2)
    (ax1 == 2) && (ax1 = 1)
    (ax2 == 1) && (ax2 = 2)
    (ax2 == 2) && (ax2 = 1)
  end

  # Find the through-plane axis of the volume
  ax3 = setdiff(1:3, [ax1, ax2])[1]

  # Set names and colors of axis labels
  if plane == 'a'			# Axial plane: A->P, R->L
    label1 = ["A", "P"]
    color1 = julia_green

    label2 = ["R", "L"]
    color2 = julia_red
  elseif plane == 's'			# Sagittal plane: S->I, P->A
    label1 = ["S", "I"]
    color1 = julia_blue

    label2 = ["P", "A"]
    color2 = julia_green
  else					# Coronal plane: S->I, R->L
    label1 = ["S", "I"]
    color1 = julia_blue

    label2 = ["R", "L"]
    color2 = julia_red
  end

  # Convert volume to RGB/Gray array
  rgb = vol_to_rgb(mri.vol)

  # Set the size of the canvas that the volume will be displayed in
  csize = size(rgb)[[ax1, ax2]]
  csize = (300, Int.(round.(300 * csize[2] / csize[1])))

  # Display volume
  gui = ImageView.imshow(rgb, axes=(ax1,ax2), flipy=flip1, flipx=flip2,
                              canvassize=csize, name=mri.fspec)

  # Move display to middle slice
  push!(gui["roi"]["slicedata"].signals[1], Int(size(mri.vol, ax3)/2))

  # Add annotation for image axes
  fsize = round(5 * maximum(mri.volsize) / 256)

  ImageView.annotate!(gui,
            AnnotationText(size(mri.vol, ax2)*.5, size(mri.vol, ax1)*.05,
            label1[1], color=color1, fontsize=fsize))
  ImageView.annotate!(gui,
            AnnotationText(size(mri.vol, ax2)*.5, size(mri.vol, ax1)*.95,
            label1[2], color=color1, fontsize=fsize))
  ImageView.annotate!(gui,
            AnnotationText(size(mri.vol, ax2)*.05, size(mri.vol, ax1)*.5,
            label2[1], color=color2, fontsize=fsize))
  ImageView.annotate!(gui,
            AnnotationText(size(mri.vol, ax2)*.95, size(mri.vol, ax1)*.5,
            label2[2], color=color2, fontsize=fsize))

  # Add annotation for b-value and gradient vector
  blabel = Vector{String}(undef, 0)

  if !isempty(mri.bval)
    blabel = "b=" .* string.(Int.(round.(mri.bval)))
  end

  if !isempty(mri.bvec)
    blabel = blabel .* "\ng=[" .*
                       string.(round.(mri.bvec[:,1] * 100)/100) .* "," .*
                       string.(round.(mri.bvec[:,2] * 100)/100) .* "," .*
                       string.(round.(mri.bvec[:,3] * 100)/100) .* "]"
  end

  for ib = 1:length(blabel)
    ImageView.annotate!(gui,
              AnnotationText(size(mri.vol, ax2)*.1, size(mri.vol, ax1)*.1,
              blabel[ib], t = ib,
              color=RGB(1,1,1), fontsize=fsize, halign="left"))
  end

  return gui
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
struct ConvexPolygon{T,N,M} <: ParametricSurface{T,N}
    plane::Plane{T,N}
    origin::SVector{N,T}
    vertices::SMatrix{N,M,T}
end

# uv(surface::ParametricSurface{T,N}, p::SVector{N,T}) -> SVector{2,T}
# uvrange(surface::ParametricSurface{T,N}) -> Tuple{Tuple{T,T},Tuple{T,T}}
# point(surface::ParametricSurface{T,N}, u::T, v::T) -> SVector{N,T}
# partials(surface::ParametricSurface{T,N}, u::T, v::T) -> Tuple{SVector{N,T}, SVector{N,T}}
# normal(surface::ParametricSurface{T,N}, u::T, v::T) -> SVector{N,T}
# inside(surface::ParametricSurface{T,N}, p: :SVector{N,T}) -> Bool
# onsurface(surface::ParametricSurface{T,N}, p::SVector{N,T}) -> Bool
# surfaceintersection(surface::ParametricSurface{T,N}, AbstractRay::Ray{T,N}) -> Union{EmptyInterval{T},Interval{T},DisjointUnion{T}}
# interface(surface::ParametricSurface{T,N}) -> OpticalInterface{T}


text:
"""

General Framework for the data augmented Gaussian Processes

"""
module AugmentedGaussianProcesses

export AbstractGP, GP, VGP, SVGP
export Likelihood,  RegressionLikelihood, ClassificationLikelihood, MultiClassLikelihood
export GaussianLikelihood, StudentTLikelihood, LaplaceLikelihood
export LogisticLikelihood, BayesianSVM
export SoftMaxLikelihood, LogisticSoftMaxLikelihood
export PoissonLikelihood
export Inference, Analytic, AnalyticVI, AnalyticSVI, GibbsSampling, MCIntegrationVI, MCIntegrationSVI, QuadratureVI, QuadratureSVI
export NumericalVI, NumericalSVI

#Deprecated
export BatchGPRegression, SparseGPRegression, MultiClass, SparseMultiClass, BatchBSVM, SparseBSVM, BatchXGPC, SparseXGPC, BatchStudentT, SparseStudentT

#Useful functions and module
include("kernels/KernelModule.jl")
include("kmeans/KMeansModule.jl")
include("functions/PGSampler.jl")
#include("functions/PerturbativeCorrection.jl")
# include("functions/GPAnalysisTools.jl")
# include("functions/IO_model.jl")
#Custom modules
using .KernelModule
using .KMeansModule
using .PGSampler
# using .PerturbativeCorrection
# using .GPAnalysisTools
# using .IO_model
#General modules
# using MLKernels
using GradDescent
export Optimizer, Adam, VanillaGradDescent, ALRSVI, InverseDecay
using DataFrames
using Distributions
using LinearAlgebra
using StatsBase
using StatsFuns
using SpecialFunctions
using Expectations
using Random
import Base: convert, show, copy
#Exported modules
export KMeansModule
#Useful functions
export train!
export predict_f, predict_y, proba_y
# export getLog, getMultiClassLog
export Kernel,  Matern3_2Kernel, Matern5_2Kernel, RBFKernel, LaplaceKernel, SigmoidKernel, PolynomialKernel, ARDKernel
export kerneldiagmatrix, kerneldiagmatrix!, kernelmatrix, kernelmatrix!
export fstar, multiclasspredictproba, multiclasspredictprobamcmc, multiclasspredict, ELBO
export setvalue!,getvalue,setfixed!,setfree!,getvariance,getlengthscales,setoptimizer!
export opt_diag, opt_trace
export rand
export KMeansInducingPoints

# Main classes
abstract type Inference{T<:Real} end
abstract type Likelihood{T<:Real}  end

const LatentArray = Vector #For future optimization : How collection of latent GP parameters and local variables are stored
include("prior/meanprior.jl")

include("models/AbstractGP.jl")
include("models/GP.jl")
include("models/VGP.jl")
include("models/SVGP.jl")

include("inference/inference.jl")
include("likelihood/likelihood.jl")

include("functions/utils.jl")
include("functions/init.jl")
include("functions/KLdivergences.jl")
# include("functions/plotting.jl")
#Deprecated constructors
include("deprecated.jl")

#Training Functions
include("training.jl")
include("autotuning.jl")
include("predictions.jl")
end #End Module


text:
# Tests for Profit DEA Models
@testset "ProfitDEAModel" begin

    ## Test Profit DEA Model
    # Test with Zofio, Pastor and Aparicio (2013) data
    X = [1 1; 1 1; 0.75 1.5; 0.5 2; 0.5 2; 2 2; 2.75 3.5; 1.375 1.75]
    Y = [1 11; 5 3; 5 5; 2 9; 4 5; 4 2; 3 3; 4.5 3.5]
    P = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1]
    W = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1]

    GxGydollar = 1 ./ (sum(P, dims = 2) + sum(W, dims = 2))
    GxGydollar = repeat(GxGydollar, 1, 2)

    deaprofitdollar = deaprofit(X, Y, W, P, GxGydollar, GxGydollar)
    @test efficiency(deaprofitdollar, :Economic)   ≈ [2; 2; 0; 2; 2; 8; 12; 4] atol = 1e-3
    @test efficiency(deaprofitdollar, :Technical)  ≈ [0; 0; 0; 0; 0; 6; 12; 3] atol = 1e-3
    @test efficiency(deaprofitdollar, :Allocative) ≈ [2; 2; 0; 2; 2; 2; 0; 1] atol = 1e-3

    # Print
    show(IOBuffer(), deaprofitdollar)

    # Test errors
    @test_throws ErrorException deaprofit([1; 2 ; 3], [4 ; 5], [1; 1; 1], [4; 5], [1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws ErrorException deaprofit([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], [4; 5; 6], [1; 2; 3], [4; 5; 6]) # Different number of observation in input prices
    @test_throws ErrorException deaprofit([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6; 7], [1; 2; 3], [4; 5; 6]) # Different number of observation in output prices
    @test_throws ErrorException deaprofit([1 1; 2 2; 3 3], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [4; 5; 6], [1 1; 2 2; 3 3], [4; 5; 6]) # Different number of input prices and inputs
    @test_throws ErrorException deaprofit([1; 2; 3], [4 4; 5 5; 6 6], [1; 2; 3], [4 4 4; 5 5 5; 6 6 6], [1; 2; 3], [4 4; 5 5; 6 6]) # Different number of oputput prices and outputs
    @test_throws ErrorException deaprofit([1 1; 2 2; 3 3], [4; 5; 6], [1 1; 2 2; 3 3], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [4; 5; 6]) # Different size of inputs direction
    @test_throws ErrorException deaprofit([1; 2; 3], [4 4; 5 5; 6 6], [1; 2; 3], [4 4; 5 5; 6 6], [1; 2; 3], [4 4 4; 5 5 5; 6 6 6]) # Different size of inputs direction

end


text:
using Pkg

dependencies = [
    "MySQL",
    "DBInterface",
    "Bcrypt",
    "DotEnv"
]

Pkg.add(dependencies)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libjlnode_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libjlnode")
JLLWrappers.@generate_main_file("libjlnode", UUID("e3ef64f0-f261-50ad-b884-29c8a26457f8"))
end  # module libjlnode_jll