text:
using DynamicGrids, Dispersal, Test, Colors, FieldDefaults
import DynamicGrids: grid2image

init =  [1.0  1.0  0.5
         0.0  0.0  0.0
         0.0  0.0  0.0]

regionlookup = [1  2  3
                1  2  3
                2  2  0]

occurance = Bool[0  1  0
                 1  0  1
                 1  1  0]

start = 1
detectionthreshold = 0.0
framesperstep = 2
objective = RegionObjective(detectionthreshold, regionlookup, occurance, framesperstep, start)

@testset "region output works" begin
    rule = ExactExponentialGrowth(intrinsicrate = log(2.0))
    output = RegionOutput(init; objective=objective, tspan=1:5)
    sim!(output, rule)
end

@testset "image processor colors regions by fit" begin
    image = [RGB24(0.2,0.2,0.2)  RGB24(1.0,1.0,1.0)  RGB24(0.5,0.5,0.5)
             RGB24(1.0,0.0,1.0)  RGB24(1.0,1.0,0.0)  RGB24(1.0,1.0,0.0)
             RGB24(1.0,1.0,0.0)  RGB24(1.0,1.0,0.0)  RGB24(0.0,1.0,1.0)]

    truescheme = Greyscale()
    falsescheme = Greyscale(min=nothing, max=0.2)
    truezerocolor =  (1.0, 0.0, 1.0)
    falsezerocolor = (1.0, 1.0, 0.0)
    maskcolor =      (0.0, 1.0, 1.0)

    processor = ColorRegionFit(objective, truescheme, falsescheme,
                               truezerocolor, falsezerocolor, maskcolor)
    output = DynamicGrids.NoDisplayImageOutput(init; tspan=1:5)

    # This output needs an update anyway
    @test image == grid2image(processor, output, Ruleset(), (_default_=init,), 1)
end


text:
include("cp_layer/cp_layer.jl")

"""
    DefaultStateRepresentation{F, TS}

This is the default representation used by SeaPearl unless the user define his own.

It consists in a tripartite graph representation of the CP Model, with features associated with each node
and an index specifying the variable that should be branched on.
"""
mutable struct DefaultStateRepresentation{F, TS} <: FeaturizedStateRepresentation{F, TS}
    cplayergraph::CPLayerGraph
    nodeFeatures::Union{Nothing, AbstractMatrix{Float32}}
    globalFeatures::Union{Nothing, AbstractVector{Float32}}
    variableIdx::Union{Nothing, Int64}
    allValuesIdx::Union{Nothing, Vector{Int64}}
end

function DefaultStateRepresentation{F, TS}(model::CPModel; action_space=nothing) where {F, TS}
    g = CPLayerGraph(model)
    allValuesIdx = nothing
    if !isnothing(action_space)
        allValuesIdx = indexFromCpVertex.([g], ValueVertex.(action_space))
    end
    sr = DefaultStateRepresentation{F, TS}(g, nothing, nothing, nothing, allValuesIdx)
    sr.nodeFeatures = featurize(sr)
    sr.globalFeatures = global_featurize(sr)
    return sr
end

function DefaultTrajectoryState(sr::DefaultStateRepresentation{F, DefaultTrajectoryState}) where F
    if isnothing(sr.variableIdx)
        throw(ErrorException("Unable to build a DefaultTrajectoryState, when the branching variable is nothing."))
    end
    adj = Matrix(adjacency_matrix(sr.cplayergraph))
    fg = isnothing(sr.globalFeatures) ?
        FeaturedGraph(adj; nf=sr.nodeFeatures) : FeaturedGraph(adj; nf=sr.nodeFeatures, gf=sr.globalFeatures)
    return DefaultTrajectoryState(fg, sr.variableIdx, sr.allValuesIdx)
end

"""
    update_representation!(sr::DefaultStateRepresentation, model::CPModel, x::AbstractIntVar)

Update the StateRepesentation according to its Type and Featurization.
"""
function update_representation!(sr::DefaultStateRepresentation, model::CPModel, x::AbstractIntVar)
    update_features!(sr, model)
    sr.variableIdx = indexFromCpVertex(sr.cplayergraph, VariableVertex(x))
    return sr
end

struct DefaultFeaturization <: AbstractFeaturization end

"""
    featurize(sr::DefaultStateRepresentation{DefaultFeaturization, TS})

Create features for every node of the graph. Supposed to be overwritten.

Default behavior consists in a 3D One-hot vector that encodes whether the node represents a Constraint, a Variable or a Value.
"""
function featurize(sr::FeaturizedStateRepresentation{DefaultFeaturization, TS}) where TS
    g = sr.cplayergraph
    features = zeros(Float32, 3, nv(g))
    for i in 1:nv(g)
        cp_vertex = SeaPearl.cpVertexFromIndex(g, i)
        if isa(cp_vertex, ConstraintVertex)
            features[1, i] = 1.0f0
        end
        if isa(cp_vertex, VariableVertex)
            features[2, i] = 1.0f0
        end
        if isa(cp_vertex, ValueVertex)
            features[3, i] = 1.0f0
        end
    end
    features
end

"""
    feature_length(gen::AbstractModelGenerator, ::Type{FeaturizedStateRepresentation})

Returns the length of the feature vector, for the `DefaultFeaturization`.
"""
feature_length(::Type{<:FeaturizedStateRepresentation{DefaultFeaturization, TS}}) where TS = 3

DefaultStateRepresentation(m::CPModel) = DefaultStateRepresentation{DefaultFeaturization, DefaultTrajectoryState}(m::CPModel)

ext: jl


text:
using PkgBenchmark

current = BenchmarkConfig(id = "transpose_storage", juliacmd = `julia -O3`)
baseline = BenchmarkConfig(id = "master", juliacmd = `julia -O3`)
results = judge("YaoArrayRegister", current, baseline)
export_markdown("report.md", results)


text:
using FastTransforms, Compat
using Compat.Test

if VERSION < v"0.7-"
    @testset "Test vector field transforms" begin
        # f = (θ,φ) -> cospi(θ) + sinpi(θ)*(1+cospi(2θ))*sinpi(φ) + sinpi(θ)^5*(cospi(5φ)-sinpi(5φ))
        ∇θf = (θ,φ) -> -sinpi(θ) + (cospi(θ)*(1+cospi(2θ)) - 2*sinpi(θ)*sinpi(2θ))*sinpi(φ) + 5*sinpi(θ)^4*cospi(θ)*(cospi(5φ)-sinpi(5φ))
        ∇φf = (θ,φ) -> (1+cospi(2θ))*cospi(φ) - 5*sinpi(θ)^4*(sinpi(5φ)+cospi(5φ))

        n = 6
        θ = (0.5:n-0.5)/n
        φ = (0:2n-2)*2/(2n-1)
        ∇θF = [∇θf(θ,φ) for θ in θ, φ in φ]
        ∇φF = [∇φf(θ,φ) for θ in θ, φ in φ]
        V1 = zero(∇θF)
        V2 = zero(∇φF)
        Pa = FastTransforms.plan_analysis(∇θF)
        mul!(V1, V2, Pa, ∇θF, ∇φF)
        P = SlowSphericalHarmonicPlan(V1)

        U1 = zero(V1)
        U2 = zero(V2)
        At_mul_B!(U1, U2, P, V1, V2)

        W1 = zero(U1)
        W2 = zero(U2)

        mul!(W1, W2, P, U1, U2)

        Ps = FastTransforms.plan_synthesis(W1)

        G1 = zero(∇θF)
        G2 = zero(∇φF)

        mul!(G1, G2, Ps, W1, W2)

        @test vecnorm(∇θF - G1)/vecnorm(∇θF) < n*eps()
        @test vecnorm(∇φF - G2)/vecnorm(∇φF) < n*eps()

        y = (1.0, 2.0, 3.0)
        for k in (10, 20, 40)
            ∇θf = (θ,φ) -> -2k*sin(k*((sinpi(θ)*cospi(φ) - y[1])^2 + (sinpi(θ)*sinpi(φ) - y[2])^2 + (cospi(θ) - y[3])^2))*( (sinpi(θ)*cospi(φ) - y[1])*(cospi(θ)*cospi(φ)) + (sinpi(θ)*sinpi(φ) - y[2])*(cospi(θ)*sinpi(φ)) - (cospi(θ) - y[3])*sinpi(θ) )
            ∇φf = (θ,φ) -> -2k*sin(k*((sinpi(θ)*cospi(φ) - y[1])^2 + (sinpi(θ)*sinpi(φ) - y[2])^2 + (cospi(θ) - y[3])^2))*( (sinpi(θ)*cospi(φ) - y[1])*(-sinpi(φ)) + (sinpi(θ)*sinpi(φ) - y[2])*(cospi(φ)) )
            n = 12k

            θ = (0.5:n-0.5)/n
            φ = (0:2n-2)*2/(2n-1)
            ∇θF = [∇θf(θ,φ) for θ in θ, φ in φ]
            ∇φF = [∇φf(θ,φ) for θ in θ, φ in φ]
            V1 = zero(∇θF)
            V2 = zero(∇φF)
            Pa = FastTransforms.plan_analysis(∇θF)
            mul!(V1, V2, Pa, ∇θF, ∇φF)
            P = SlowSphericalHarmonicPlan(V1)

            U1 = zero(V1)
            U2 = zero(V2)
            At_mul_B!(U1, U2, P, V1, V2)

            W1 = zero(U1)
            W2 = zero(U2)

            mul!(W1, W2, P, U1, U2)

            Ps = FastTransforms.plan_synthesis(W1)

            G1 = zero(∇θF)
            G2 = zero(∇φF)

            mul!(G1, G2, Ps, W1, W2)

            @test vecnorm(∇θF - G1)/vecnorm(∇θF) < n*eps()
            @test vecnorm(∇φF - G2)/vecnorm(∇φF) < n*eps()
        end
    end
else
    @testset "Test vector field transforms" begin
        # f = (θ,φ) -> cospi(θ) + sinpi(θ)*(1+cospi(2θ))*sinpi(φ) + sinpi(θ)^5*(cospi(5φ)-sinpi(5φ))
        n = 6
        θ = (0.5:n-0.5)/n
        φ = (0:2n-2)*2/(2n-1)
        ∇θf = (θ,φ) -> -sinpi(θ) + (cospi(θ)*(1+cospi(2θ)) - 2*sinpi(θ)*sinpi(2θ))*sinpi(φ) + 5*sinpi(θ)^4*cospi(θ)*(cospi(5φ)-sinpi(5φ))
        ∇φf = (θ,φ) -> (1+cospi(2θ))*cospi(φ) - 5*sinpi(θ)^4*(sinpi(5φ)+cospi(5φ))

        ∇θF = [∇θf(θ,φ) for θ in θ, φ in φ]
        ∇φF = [∇φf(θ,φ) for θ in θ, φ in φ]
        V1 = zero(∇θF)
        V2 = zero(∇φF)
        Pa = FastTransforms.plan_analysis(∇θF)
        mul!(V1, V2, Pa, ∇θF, ∇φF)
        P = SlowSphericalHarmonicPlan(V1)

        U1 = zero(V1)
        U2 = zero(V2)
        mul!(U1, U2, transpose(P), V1, V2)

        W1 = zero(U1)
        W2 = zero(U2)

        mul!(W1, W2, P, U1, U2)

        Ps = FastTransforms.plan_synthesis(W1)

        G1 = zero(∇θF)
        G2 = zero(∇φF)

        mul!(G1, G2, Ps, W1, W2)

        @test norm(∇θF - G1)/norm(∇θF) < n*eps()
        @test norm(∇φF - G2)/norm(∇φF) < n*eps()

        for k in (10, 20, 40)
            y = (1.0, 2.0, 3.0)
            n = 12k
            θ = (0.5:n-0.5)/n
            φ = (0:2n-2)*2/(2n-1)
            ∇θf = (θ,φ) -> -2k*sin(k*((sinpi(θ)*cospi(φ) - y[1])^2 + (sinpi(θ)*sinpi(φ) - y[2])^2 + (cospi(θ) - y[3])^2))*( (sinpi(θ)*cospi(φ) - y[1])*(cospi(θ)*cospi(φ)) + (sinpi(θ)*sinpi(φ) - y[2])*(cospi(θ)*sinpi(φ)) - (cospi(θ) - y[3])*sinpi(θ) )
            ∇φf = (θ,φ) -> -2k*sin(k*((sinpi(θ)*cospi(φ) - y[1])^2 + (sinpi(θ)*sinpi(φ) - y[2])^2 + (cospi(θ) - y[3])^2))*( (sinpi(θ)*cospi(φ) - y[1])*(-sinpi(φ)) + (sinpi(θ)*sinpi(φ) - y[2])*(cospi(φ)) )

            ∇θF = [∇θf(θ,φ) for θ in θ, φ in φ]
            ∇φF = [∇φf(θ,φ) for θ in θ, φ in φ]
            V1 = zero(∇θF)
            V2 = zero(∇φF)
            Pa = FastTransforms.plan_analysis(∇θF)
            mul!(V1, V2, Pa, ∇θF, ∇φF)
            P = SlowSphericalHarmonicPlan(V1)

            U1 = zero(V1)
            U2 = zero(V2)
            FastTransforms.mul!(U1, U2, transpose(P), V1, V2)

            W1 = zero(U1)
            W2 = zero(U2)

            mul!(W1, W2, P, U1, U2)

            Ps = FastTransforms.plan_synthesis(W1)

            G1 = zero(∇θF)
            G2 = zero(∇φF)

            mul!(G1, G2, Ps, W1, W2)

            @test norm(∇θF - G1)/norm(∇θF) < n*eps()
            @test norm(∇φF - G2)/norm(∇φF) < n*eps()
        end
    end
end


text:
module QuizQuestions

using Markdown
using Mustache
using Random

include("question_types.jl")
include("html_templates.jl")
include("show_methods.jl")

export numericq, radioq, multiq, booleanq, yesnoq, stringq

end


text:
#=
see description in native.jl

allows objective to minimize frobenius norm or nuclear norm of residual matrix,
plus regularization; similar to
https://arxiv.org/ftp/arxiv/papers/1405/1405.1207.pdf
=#

using SparseArrays

struct MatrixRegressionJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    Y::Matrix{T}
    X::Matrix{T}
    nuc_obj::Bool # use nuclear norm loss, else squared loss
    lam_fro::Real # penalty on Frobenius norm
    lam_nuc::Real # penalty on nuclear norm
    lam_las::Real # penalty on l1 norm
    lam_glr::Real # penalty on penalty on row group l1 norm
    lam_glc::Real # penalty on penalty on column group l1 norm
end

function MatrixRegressionJuMP{Float64}(
    n::Int,
    m::Int,
    p::Int,
    args...;
    A_max_rank::Int = div(m, 2) + 1,
    A_sparsity::Real = max(0.2, inv(sqrt(m * p))),
    Y_noise::Real = 0.01,
    )
    @assert p >= m
    @assert 1 <= A_max_rank <= m
    @assert 0 < A_sparsity <= 1
    A_left = sprandn(p, A_max_rank, A_sparsity)
    A_right = sprandn(A_max_rank, m, A_sparsity)
    A = 10 * A_left * A_right
    X = randn(n, p)
    Y = X * A + Y_noise * randn(n, m)
    return MatrixRegressionJuMP{Float64}(Y, X, args...)
end

function MatrixRegressionJuMP{Float64}(n::Int, m::Int)
    @assert n >= m >= 1
    X = randn(n, m)
    Y = randn(n, m)
    return MatrixRegressionJuMP{Float64}(Y, X, true, 0.1, 0, 0, 0, 0)
end

function build(inst::MatrixRegressionJuMP{T}) where {T <: Float64}
    (Y, X) = (inst.Y, inst.X)
    @assert min(inst.lam_fro, inst.lam_nuc, inst.lam_las,
        inst.lam_glr, inst.lam_glc) >= 0
    (data_n, data_m) = size(Y)
    data_p = size(X, 2)
    @assert size(X, 1) == data_n
    @assert data_p >= data_m
    @assert !inst.nuc_obj || data_m <= data_n
    model = JuMP.Model()
    JuMP.@variable(model, A[1:data_p, 1:data_m])
    JuMP.@variable(model, loss)
    loss_mat = Y - X * A

    if inst.nuc_obj
        JuMP.@constraint(model, vcat(loss, vec(loss_mat')) in
            MOI.NormNuclearCone(data_m, data_n))
    else
        if data_n > data_p
            # dimension reduction via QR
            F = qr(X, ColumnNorm())
            loss_mat = (F.Q' * Y)[1:data_p, :] - F.R[1:data_p, 1:data_p] *
                F.P' * A
        end
        JuMP.@constraint(model, vcat(loss, 1, vec(loss_mat) / sqrt(data_n)) in
            JuMP.RotatedSecondOrderCone())
    end

    obj = one(T) * loss
    if !iszero(inst.lam_fro)
        JuMP.@variable(model, t_fro)
        JuMP.@constraint(model, vcat(t_fro, inst.lam_fro * vec(A)) in
            JuMP.SecondOrderCone())
        obj += t_fro
    end
    if !iszero(inst.lam_nuc)
        JuMP.@variable(model, t_nuc)
        JuMP.@constraint(model, vcat(t_nuc, inst.lam_nuc * vec(A)) in
            MOI.NormNuclearCone(data_p, data_m))
        obj += t_nuc
    end
    if !iszero(inst.lam_las)
        JuMP.@variable(model, t_las)
        JuMP.@constraint(model, vcat(t_las, inst.lam_las * vec(A)) in
            MOI.NormOneCone(data_p * data_m + 1))
        obj += t_las
    end
    if !iszero(inst.lam_glr)
        JuMP.@variable(model, t_glr[1:data_p])
        JuMP.@constraint(model, [i = 1:data_p], vcat(t_glr[i],
            inst.lam_glr * A[i, :]) in JuMP.SecondOrderCone())
        obj += sum(t_glr)
    end
    if !iszero(inst.lam_glc)
        JuMP.@variable(model, t_glc[1:data_m])
        JuMP.@constraint(model, [i = 1:data_m], vcat(t_glc[i],
            inst.lam_glc * A[:, i]) in JuMP.SecondOrderCone())
        obj += sum(t_glc)
    end
    JuMP.@objective(model, Min, obj)

    model.ext[:A_var] = A # save for use in tests

    return model
end

function test_extra(inst::MatrixRegressionJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective value is correct
    tol = eps(T)^0.2
    (Y, X) = (inst.Y, inst.X)
    A_opt = JuMP.value.(model.ext[:A_var])
    loss_mat = Y - X * A_opt
    if inst.nuc_obj
        loss = sum(svdvals(loss_mat))
    else
        loss = sum(abs2, loss_mat) / (2 * size(Y, 1))
    end
    obj_result = loss +
        inst.lam_fro * norm(A_opt, 2) +
        inst.lam_nuc * sum(svdvals(A_opt)) +
        inst.lam_las * norm(A_opt, 1) +
        inst.lam_glr * sum(norm, eachrow(A_opt)) +
        inst.lam_glc * sum(norm, eachcol(A_opt))
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    return
end


text:
__precompile__()

module DiffEqDiffTools

using LinearAlgebra

include("function_wrappers.jl")
include("finitediff.jl")
include("derivatives.jl")
include("gradients.jl")
include("jacobians.jl")

end # module


text:
using Pkg

dependencies = [
    "GraphRecipes",
    "LightGraphs",
    "SimpleWeightedGraphs",
    "GraphPlot",
    "MetaGraphs",
    "Documenter",
    "StatsBase",
    "LinearAlgebra",
    "Statistics",
    "Colors",
    "Plots",
    "LaTeXStrings",
    "Cairo",
    "Compose",
    "SparseArrays",
    "QuadGK",
    "SatelliteToolbox",
    "GR"]

Pkg.add(dependencies)
Pkg.update(dependencies)

max_stars_repo_name: drpeterrohde/QuNet

max_issues_repo_name: SamuelMarks/QuNet

max_forks_repo_name: SamuelMarks/QuNet


text:
module Functionals

export F_ss, J_T_ss, chi_ss!, F_sm, J_T_sm, chi_sm!, F_re, J_T_re, chi_re!
export grad_J_T_sm!

import ..WeightedObjective

using LinearAlgebra

"""Average complex overlap of the target states with forward-propagated states.

```julia
f_tau(ϕ, objectives; τ=nothing)
```
"""
function f_tau(ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k in 1:N]
    end
    f::ComplexF64 = 0
    for k = 1:N
        obj = objectives[k]
        w = isa(obj, WeightedObjective) ? obj.weight : 1.0
        f += w * τ[k]
    end
    return f / N
end

"""State-to-state phase-insensitive fidelity.

```julia
F_ss(ϕ, objectives; τ=nothing)
```
"""
function F_ss(ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k in 1:N]
    end
    F::ComplexF64  = f_tau(ϕ, objectives; τ=abs.(τ).^2)
    @assert imag(F) < 1e-10
    return real(F)
end

"""State-to-state phase-insensitive functional.

```julia
J_T_ss(ϕ, objectives; τ=nothing)
```
"""
function J_T_ss(ϕ, objectives; τ=nothing)
    return 1.0 - F_ss(ϕ, objectives; τ=τ)
end

"""Krotov-states χ for functional [`J_T_ss`](@ref).

```julia
chi_ss!(χ, ϕ, objectives; τ=nothing)
```
"""
function chi_ss!(χ, ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k in 1:N]
    end
    for k = 1:N
        obj = objectives[k]
        ϕₖ_tgt = obj.target_state
        copyto!(χ[k], ϕₖ_tgt)
        w = isa(obj, WeightedObjective) ? obj.weight : 1.0
        lmul!((τ[k] * w)/N, χ[k])
    end
end

"""Square-modulus fidelity.

```julia
F_sm(ϕ, objectives; τ=nothing)
```
"""
function F_sm(ϕ, objectives; τ=nothing)
    return abs(f_tau(ϕ, objectives; τ=τ))^2
end

"""Square-modulus functional.

```julia
J_T_sm(ϕ, objectives; τ=nothing)
```
"""
function J_T_sm(ϕ, objectives; τ=nothing)
    return 1.0 - F_sm(ϕ, objectives; τ=τ)
end

"""Gradient for [`J_T_sm`](@ref)."""
function grad_J_T_sm!(G, τ, ∇τ)
    N = length(τ) # number of objectives
    L, N_T = size(∇τ[1])  # number of controls/time intervals
    G′ = reshape(G, L, N_T)  # writing to G′ modifies G
    for l = 1:L
        for n = 1:N_T
            G′[l, n] = real(sum([conj(τ[k′]) * ∇τ[k][l, n]
                                 for k′=1:N for k=1:N]))
        end
    end
    lmul!(-2/N^2, G)
    return G
end

"""Krotov-states χ for functional [`J_T_sm`](@ref).

```julia
chi_sm!(χ, ϕ, objectives; τ=nothing)
```
"""
function chi_sm!(χ, ϕ, objectives; τ=nothing)

    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k in 1:N]
    end

    w = ones(N)
    for k = 1:N
        obj = objectives[k]
        w[k] = isa(obj, WeightedObjective) ? obj.weight : 1.0
    end

    for k = 1:N
        obj = objectives[k]
        ϕₖ_tgt = obj.target_state
        copyto!(χ[k], ϕₖ_tgt)
        lmul!(w[k] * sum(w .* τ) / N^2, χ[k])
    end

end

"""Real-part fidelity.

```julia
F_re(ϕ, objectives; τ=nothing)
```
"""
function F_re(ϕ, objectives; τ=nothing)
    return real(f_tau(ϕ, objectives; τ=τ))
end

"""Real-part functional.

```julia
J_T_re(ϕ, objectives; τ=nothing)
```
"""
function J_T_re(ϕ, objectives; τ=nothing)
    return 1.0 - F_re(ϕ, objectives; τ=τ)
end

"""Krotov-states χ for functional [`J_T_re`](@ref).

```julia
chi_re!(χ, ϕ, objectives; τ=nothing)
```
"""
function chi_re!(χ, ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k in 1:N]
    end
    for k = 1:N
        obj = objectives[k]
        ϕₖ_tgt = obj.target_state
        copyto!(χ[k], ϕₖ_tgt)
        w = isa(obj, WeightedObjective) ? obj.weight : 1.0
        lmul!(w/(2N), χ[k])
    end
end

end


text:
@testset "action_selectors" begin
    include("alternate_selector.jl")
    include("epsilon_greedy_selector.jl")
    include("usb_selector.jl")
    include("weighted_selector.jl")
end


text:
"""
    PTModel(logdensity, logreference)
`logdensity` is the function that implements the log of the joint probability
density, i.e. the unormalized posterior.
`logreference` is the reference distribution you wish to sampler from.

**Note:**
The reference distribution is forced to be easy to sample from in order to
induce faster mixing in the lower temperature chains.
"""
struct PTModel{L,P<:Distribution} <: AbstractMCMC.AbstractModel
    logdensity::L
    reference::P
end


text:
# -------- Explicitly Typed values for automated json encoding --------------- #

struct Typed{T}
  typ :: String
  value :: T

  # Overwrite default constructor
  Typed(value) = new{typeof(value)}(typeof(value) |> string |> lowercase, value)
end

# -------- Individual Layers that compose an Image --------------------------- #

abstract type Layer end

const Point = Tuple{Float32, Float32}

struct Lines <: Layer
  data :: Vector{Tuple{Point, Point}}
  name :: String
  style :: String
end

Lines(data; name = "lines", style = "black") = Lines(data, name, style)

function Lines(xs, ys; kwargs...)
  data = map(xs, ys) do x, y
    ((x[1], y[1]), (x[2], y[2]))
  end
  Lines(data; kwargs...)
end

struct Heatmap <: Layer
  data :: Vector{Float32}
  name :: String
  style :: String
  min :: Float32
  max :: Float32
end

function Heatmap( data
                ; min = minimum(data)
                , max = maximum(data)
                , name = "heatmap"
                , style = "grey" )
  Heatmap(data, name, style, min, max)
end

struct Tokens <: Layer
  data :: Vector{String}
  name :: String
  style :: String
end

function Tokens(game :: Union{MNKGame, MetaTac}; name = "tokens", style = "black")
  tokens = Dict(1  => "X", -1 => "O", 0  => "")
  Tokens(map(x -> tokens[x], game.board), name, style)
end

struct Actions <: Layer
  data :: Vector{Int}
  name :: String
end

function Actions(game :: Union{MNKGame, MetaTac}; name = "actions")
  actions = legal_actions(game)
  action_data = -ones(Int, policy_length(game))
  action_data[actions] = actions
  Actions(action_data, name)
end

# -------- An image as a stack of layers with some geometry ------------------ #

struct Image
  layers
  width :: Int
  height :: Int
  name :: String
  value :: Union{Float32, Nothing}

  # Overwrite default constructor
  Image(layers :: Vector{Layer}, args...) = new(Typed.(layers), args...)
end

function Image( game :: Union{MNKGame, MetaTac}
              , policy :: Vector{Float32}
              , value :: Float32 )

  m, n = size(game)[1:2]

  heatmap = Heatmap(policy, min = 0., max = 1.)
  name = typeof(game) |> string |> lowercase

  layers = [heatmap, Tokens(game), Actions(game)]

  if isa(game, MetaTac)
    x = Point[(3, 3), (6, 6), (0, 9), (0, 9)]
    y = Point[(0, 9), (0, 9), (3, 3), (6, 6)]
    push!(layers, Lines(x, y))
  end

  Image(layers, m, n, name, value)
end

max_stars_repo_name: roSievers/jtac


text:
include("../src/MatrixReader.jl")

function fakedata(fakecount, blocksize, expsize; maxval=100)
    testdata = zeros(Int, (blocksize, expsize))
    for _ in 1:fakecount
        i = rand(1:blocksize)
        j = rand(1:expsize)
        v = rand(1:maxval)
        testdata[i,j] = v
    end
    testdata
end

# Testing with 10000 random data with random sparcity, dimensions.

println("For UInt8 version")

for itr in 1:10000
    expsize = rand(20:100)

    ######################
    # Write in fake data #
    ######################
    mw = MatrixWriter("test.data", expsize, UInt8)

    # data1
    t1 = fakedata(rand(1:200), 100, expsize, maxval=100)
    writeMatrix(mw, t1)
    # data2
    t2 = fakedata(rand(1:200), 100, expsize, maxval=100)
    writeMatrix(mw, t2)
    close(mw)

    #####################
    # Read in fake data #
    #####################
    mr = MatrixReader("test.data", 25, buffsize=rand(10:50))
    new = advance!(mr)
    m1 = mr.data
    new = advance!(mr)
    m2 = mr.data
    new = advance!(mr)
    m3 = mr.data
    new = advance!(mr)
    m4 = mr.data
    new = advance!(mr)
    m5 = mr.data
    new = advance!(mr)
    m6 = mr.data
    new = advance!(mr)
    m7 = mr.data
    new = advance!(mr)
    m8 = mr.data
    close(mr)

    #################
    # Test equality #
    #################
    x1 = vcat(t1, t2)
    new = vcat(m1, m2, m3, m4, m5, m6, m7, m8)
    @assert new == x1
    if itr%1000 == 0
        println(itr," test passed. ", size(x1), ":", sum(x1))
    end
end

println("For UInt16 version")

for itr in 1:10000
    expsize = rand(20:100)

    ######################
    # Write in fake data #
    ######################
    mw = MatrixWriter("test.data", expsize, UInt16)

    # data1
    t1 = fakedata(rand(1:200), 100, expsize, maxval=60000)
    writeMatrix(mw, t1)
    # data2
    t2 = fakedata(rand(1:200), 100, expsize, maxval=60000)
    writeMatrix(mw, t2)
    close(mw)

    #####################
    # Read in fake data #
    #####################
    mr = MatrixReader("test.data", 25, buffsize=rand(10:50))
    new = advance!(mr)
    m1 = mr.data
    new = advance!(mr)
    m2 = mr.data
    new = advance!(mr)
    m3 = mr.data
    new = advance!(mr)
    m4 = mr.data
    new = advance!(mr)
    m5 = mr.data
    new = advance!(mr)
    m6 = mr.data
    new = advance!(mr)
    m7 = mr.data
    new = advance!(mr)
    m8 = mr.data
    close(mr)

    #################
    # Test equality #
    #################
    x1 = vcat(t1, t2)
    new = vcat(m1, m2, m3, m4, m5, m6, m7, m8)
    @assert new == x1
    if itr%1000 == 0
        println(itr," test passed. ", size(x1), ":", sum(x1))
    end
end


text:
using SpinSymmetry
using Test

@testset "SpinSymmetry.jl" begin
    include("basis.jl")
    include("symmetries.jl")
end


text:
export @rx,
  react,
  dispatch,
  Event,
  ValueEvent,
  CompletedEvent,
  ErrorEvent,
  Observer,
  Observable,
  Reactor,
  collector,
  onEvent,
  subscribe!,
  chain!,
  notify!,
  Collector

using Distributed

# / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / /
#
# Core objects for implementing a reactive programming model
# (inspired by ReactiveX)
#
# / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / /

"""
An [`Event`](@ref) containing a value to deliver to an [`Observable`](@ref).
"""
struct ValueEvent{T}
  value::T
end

"""
An [`Event`](@ref) signifying no more events are available from the originating
[`Observable`](@ref).
"""
struct CompletedEvent end

"""
An [`Event`](@ref) signifying the originating [`Observable`](@ref) encountered an error;
no more events will be delivered after events of tyie.
"""
struct ErrorEvent{E}
  error::E
end

"""
An `Event` is any of the following:
* A _value_, encapsulated as [`ValueEvent`](@ref)
* A _completion_ event, encapsulated as [`CompletedEvent`](@ref)
* An _error_, encapsulated as [`ErrorEvent`](@ref)
"""
Event = Union{ValueEvent,CompletedEvent,ErrorEvent}

"""
An `Observer` is a receiver of [`Event`](@ref)s via [`onEvent`](@ref)
"""
abstract type Observer end

"""
    onEvent(observer, event)

Deliver an [`Event`](@ref) to the `Observer`
"""
function onEvent(observer, event) end

"""
`Observers` is a collection of subscribed [`Observer`](@ref)s
"""
Observers = Vector{Observer}

"""
An `Observable` is a source of [`Event`](@ref)s, with [`Observer`](@ref)s
[`subscribe!`](@ref)ing to the `Observable` in order to receive those events.

See [Observable in ReactiveX documentation](http://reactivex.io/documentation/observable.html)
"""
mutable struct Observable
  observers::Observers
end

Observable() = Observable([])

"""
    notify!(observers, event)

Notify [`Observer`](@ref) by invoking the block on each
"""
function notify!(observers, event)
  for observer in observers
    onEvent(observer, event)
  end
end

"""
    error!(observers, err)

Emit an [`ErrorEvent`](@ref) to [`Observers`](@ref)
"""
function error!(observers, err)
  evt = ErrorEvent(err)
  notify!(observers, evt)
end

"""
    complete!(observers)

Emit a [`CompletedEvent`](@ref) to [`Observers`](@ref)
"""
function complete!(observers)
  evt = CompletedEvent()
  notify!(observers, evt)
end

"""
    subscribe!(observable::Observable, observer)

Subscribe an [`Observer`](@ref) to the given [`Observable`](@ref), return the [`Observer`](@ref)
"""
function subscribe!(observable::Observable, observer)
  register!(observable.observers, observer)
end

"""
    register!(observers::Observers, observer)
Add an [`Observer`](@ref) to an existing set of [`Observers`](@ref), return the [`Observer`](@ref)
"""
function register!(observers::Observers, observer)
  push!(observers, observer)
  observer
end

"""
    chain!(observable, observer)

Subscribe the [`Observer`](@ref) to the [`Observable`](@ref) using [`subscribe!`](@ref),
and return the [`Observable`](@ref).
"""
function chain!(observable, observer)
  subscribe!(observable, observer)
  observable
end

"""
    Reactor() = Reactor(notify!)
    Reactor(fn) = Reactor(fn, [])

A `Reactor` is an [`Observer`](@ref) that is also useful for building [`Observer`](@ref)s that
are [`Observable`](@ref) as well.

A `Reactor` is the analogue of a [`Subject`](http://reactivex.io/documentation/subject.html) in the ReactiveX model.
"""
mutable struct Reactor <: Observer
  """
  Function with arguments `(observers::Observers, value)` to handle
  each `ValueEvent` for the `Reactor`.
  """
  fn::Any
  """
  The `Observers` subscribed to this `Reactor`
  """
  observers::Observers
end

Reactor() = Reactor(notify!)
Reactor(fn) = Reactor(fn, [])

"""
    dispatch(fn::Function)::Reactor = Reactor(fn)

Return a [`Reactor`](@ref) around the provided function which takes [`Observers`](@ref)
and an [`Event`](@ref), taking action as appropriate, and producing an [`Observable`](@ref)
that can also be an [`Observer`](@ref) of other [`Observable`](@ref)s.
"""
dispatch(fn::Function)::Reactor = Reactor(fn)

"""
    react(fn::Function)::Reactor
Return a [`Reactor`](@ref) around the provided function which takes [`Observers`](@ref)
and a value, taking action as appropriate, and producing an [`Observable`](@ref)
that can also be an [`Observer`](@ref) of other [`Observable`](@ref)s. [`CompletedEvent`](@ref)s
and [`ErrorEvent`](@ref)s are passed on to [`Observers`](@ref), while [`ValueEvent`](@ref)s
are passed to the supplied function. The supplied function takes [`Observers`](@ref) and a value
as arguments.
"""
react(fn::Function)::Reactor =
  dispatch() do observers::Observers, event::Event
    if isa(event, ValueEvent)
      fn(observers, event.value)
    else
      notify!(observers, event)
    end
  end

function onEvent(reactor::Reactor, event::Event)
  reactor.fn(reactor.observers, event)
end

function subscribe!(reactor::Reactor, observer)
  register!(reactor.observers, observer)
end

"""
An [`Observer`](@ref) that accumulates events on a `Channel`, which then may be
retrieved for iteration
"""
struct Collector <: Observer
  events::RemoteChannel
end

"""
    collector(buffer = 0)
Return a [`Collector`](@ref) to accumulate observed [`Event`](@ref)s. If set to a value
greater than 0, then events will be buffered with a `Channel` with the specified buffer size,
throttling incoming events until the channel has empty slots.
"""
function collector(buffer = 0)
  Collector(RemoteChannel( () -> Channel(buffer)) )
end

function onEvent(collector::Collector, event::Event)
  put!(collector.events, event)
  if !isa(event, ValueEvent)
    close(collector.events)
  end
end

function Base.iterate(collector::Collector, _s = nothing)
  evt = try
    # We have to catch an exception in case the channel i closed
    take!(collector.events)
  catch e
    if isa(e, InvalidStateException)
      return nothing
    else
      throw(e)
    end
  end
  if isa(evt, CompletedEvent)
    nothing
  elseif isa(evt, ValueEvent)
    return (evt.value, collector)
  elseif isa(evt, ErrorEvent)
    throw(evt.error)
  else
    error("Unknown event $show(evt)")
  end
end

function Base.IteratorSize(collector::Collector)
  Base.SizeUnknown()
end

"""
    rx(blk,source=nothing; buffer=0)

Given a do-block where each statement is an [`Observable`](@ref),
[`subscribe!`](@ref) to each in sequence to the one proceeding using [`chain!`](@ref).
Return an object that one can use to iterate over the events from
the last [`Observable`](@ref) in the block.

If source (a positional argument) is provided and is not `nothing`, then it is used as the
initial  [`Observable`](@ref) in the pipeline; this can be useful if the result of one
pipeline needs to be fed into another pipeline as the initial [`Observable`](@ref).

The value of `buffer` (a keyword argument is passed to [`collector`](@ref) in order to
create a potentially buffered `Channel` for the [`Collector`](@ref) at
the end of the pipeline.

Example:

```@example
results = @rx() do
count
filter(:even)
end
printlin([result for result in results])
```
"""
macro rx(args...)
  # the blk is actually an expression of type :-> (closure),
  # and its first arg is the args list (should be () ), followed by
  # the array of statements in the closure's block
  blk, posargs, kwargs = macroArguments(args)
  # use the buffer keyword to choose a buffer size for the collector
  buffer = get(kwargs, :buffer, 0)
  steps = Vector(reverse(blk.args[2].args))
  # if there is a 2nd argument to the macro, regard it
  # as an initial source for the pipeline
  if length(posargs) >= 1
    push!(steps, posargs[1])
  end
  pipes = map(steps) do p
    if isa(p, LineNumberNode)
      p
    else
      step = repr(p)
      :( it = chain!($(esc(p)), it) )
    end
  end
  return quote
    let evts = collector($(esc(buffer)))
      fn = function(evts)
        try
          it = evts
          $(pipes...)
        catch e
          @error "Encountered error in pipeline; $e"
        end
      end
      @anywhere fn(evts)
      evts
    end
  end
end

"""
    generate_value!(value, observer)

Treat a value as an [`Observable`](@ref), emitting the value followed by a [`CompletedEvent`](@ref)
"""
function generate_value!(value, observer)
  begin
    try
      evt = ValueEvent(value)
      onEvent(observer, evt)
      onEvent(observer, CompletedEvent())
    catch e
      println(stderr, "Caught error $e")
      for (exc, bt) in Base.catch_stack()
        showerror(stdout, exc, bt)
        println()
      end
      onEvent(observer, ErrorEvent(e))
    end
  end
end

"""
    generate_iterable!(iterable, observer)

Treat an iteraable as an [`Observable`](@ref), emitting each of its elements in turn
in a [`ValueEvent`](@ref), and concluding with a [`CompletedEvent`](@ref).
"""
function generate_iterable!(iterable, observer)
  begin
    try
      for item in iterable
        evt = ValueEvent(item)
        onEvent(observer, evt)
      end
      onEvent(observer, CompletedEvent())
    catch e
      println(stderr, "Caught error $e")
      for (exc, bt) in Base.catch_stack()
        showerror(stdout, exc, bt)
        println()
      end
      onEvent(observer, ErrorEvent(e))
    end
  end
end

function generate!(observable, observer)
  if applicable(iterate, observable)
    generate_iterable!(observable, observer)
  else
    generate_value!(observable, observer)
  end
end

function generate!(value::String, observer)
  generate_value!(value, observer)
end

"""
Treat any value as an [`Observable`](@ref), and subscribe the [`Observer`](@ref) to it; the value
will be emitted once with a [`ValueEvent`](@ref), then a [`CompletedEvent`](@ref) will be emitted. If the
value is iterable (e.g., `applicable(iterate, value)` returns `true`), then emit
over each value returned by the iterable and pass on to the [`Observer`](@ref), concluding
with a [`CompletedEvent`](@ref). Note that strings are treated as an atomic value and not
an iterable.

```jldoctest
# works for booleans
julia> evts = @rx(() -> true)
Collector(Channel{Any}(sz_max:32,sz_curr:2))

julia> collect(evts)
1-element Array{Any,1}:
 true
```

```jldoctest
# works for integers
julia> evts = @rx( ()-> 1 )
Collector(Channel{Any}(sz_max:32,sz_curr:2))

julia> collect(evts)
1-element Array{Any,1}:
 1
```

```jldoctest
# works for strings
julia> evts = @rx( ()-> "foo" )
Collector(Channel{Any}(sz_max:32,sz_curr:2))

julia> collect(evts)
1-element Array{Any,1}:
 "foo"
```

```jldoctest
# works for iterables
julia> evts = @rx( ()-> [1 2 3] )
Collector(Channel{Any}(sz_max:32,sz_curr:4))

julia> collect(evts)
3-element Array{Any,1}:
 1
 2
 3
```
"""
function subscribe!(observable, observer)
  generate!(observable, observer)
  observer
end


text:
# This file was generated, do not modify it. # hide
Xc = coerce(X, :Proline=>Continuous, :Magnesium=>Continuous);


text:
export OptAlgo, PopulationBasedOptAlgo, best

abstract OptAlgo{C}

abstract PopulationBasedOptAlgo{C} <: OptAlgo{C}

best(algo::PopulationBasedOptAlgo) = algo.popMan.best

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
function f753(N::Int, K = 0, n₃ = 0, n₅ = 0, n₇ = 0)
    if K > N
        return 0
    else
        result = (n₃ > 0 && n₅ > 0 && n₇ > 0) +
            f753(N, 10K + 3, n₃ + 1, n₅, n₇) +
            f753(N, 10K + 5, n₃, n₅ + 1, n₇) +
            f753(N, 10K + 7, n₃, n₅, n₇ + 1)
        return result
    end
end

function main()
    N = parse(Int, readline())
    println(f753(N))
end

main()


text:
module SudokuSolver

include("constants.jl")
include("types.jl")
include("operations.jl")
include("process solution.jl")

# Solve the input Sudoku puzzle.
# All solutions found are printed to STDOUT, and returned.
function solveSudoku(input::Matrix{Number})
    gs = GameState()
    numbersFilled = populateGame!(gs, input)
    placeNextNumber!(gs, numbersFilled)

    @printf "\n\nTotal COUNT = %d\n" getTotalCount(gs.levelCount)

    return gs.solutions
end

# Populate the initial state of the game using the input Sudoku puzzle.
# Return the number of cells filled up by the input.
function populateGame!(gs::GameState, input::Matrix{Number})
    numbersFilled = 0

    for i in 1:9, j in 1:9
        if input[i, j] != 0
            cell = Cell(i, j)
            setCell!(gs, cell, input[i, j])
            numbersFilled += 1
            swapSequenceEntries!(gs, numbersFilled, findfirst(gs.sequence, cell))
        end
    end

    return numbersFilled
end

# Place the next number.
# Choose a cell and run through all the possible numbers we can place in it.
# A solution is found if we filled up all the cells successfully.
function placeNextNumber!(gs::GameState, numbersFilled::Int)
    if numbersFilled ≥ 81
        # Solution found.
        processSolution!(gs)
        return
    end

    currentIndex = numbersFilled + 1
    gs.levelCount[currentIndex] += 1

    nextCellIndex = chooseNextCell(gs, currentIndex)
    swapSequenceEntries!(gs, currentIndex, nextCellIndex)

    cell = gs.sequence[currentIndex]

    leftoverNumbers = getLeftoverNumbers(gs, cell)
    for number in leftoverNumbers
        setCell!(gs, cell, number)
        placeNextNumber!(gs, numbersFilled + 1)
        clearCell!(gs, cell)
    end

    swapSequenceEntries!(gs, currentIndex, nextCellIndex)
end

# Choose the next cell to fill in number.
# Not all cell are equally good,
# we choose one that has the least possible numbers to fill in, for efficiency sake.
# This is to reduce branching, and increase chance of early termination.
function chooseNextCell(gs::GameState, currentIndex::Int)
    nextIndex = currentIndex
    MinTotalChoices = 10 #Max bitCount of any cell is 9.

    for T in currentIndex:81
        cell = gs.sequence[T]
        leftoverNumbers = getLeftoverNumbers(gs, cell)
        totalChoices = getSize(leftoverNumbers)

        if totalChoices < MinTotalChoices
            MinTotalChoices = totalChoices
            nextIndex = T
        end
    end

    return nextIndex
end

end #SudokuSolver

max_stars_repo_name: Teo-ShaoWei/SudokuSolverSG

max_issues_repo_name: Teo-ShaoWei/SudokuSolverSG

max_forks_repo_name: Teo-ShaoWei/SudokuSolverSG


text:
"""
Some tests check that KGE or KGE' is ≈ -0.41
when Q_sim is equal to μ_obs. See Knoben et al., (2019).

In other words, values > -0.41 indicates that the model
performs "better" than simply taking the mean of observations.

References
----------
1. Knoben, W.J.M., Freer, J.E., Woods, R.A., 2019.
    Technical note: Inherent benchmark or not?
    Comparing Nash-Sutcliffe and Kling-Gupta efficiency scores (preprint).
    Catchment hydrology/Modelling approaches.
    https://doi.org/10.5194/hess-2019-327
"""

using Statistics
using Test
using Streamfall

test_seq = [1,10,100,5,5,4,10,99]
b_seq = repeat([mean(test_seq)], length(test_seq))
c_seq = repeat([0.0], length(test_seq))

@testset "Kling-Gupta" begin
    @test Streamfall.KGE(test_seq, test_seq) == 1.0
    @test Streamfall.KGE(test_seq, b_seq) ≈ -0.4142135623730
    @test Streamfall.KGE(test_seq, c_seq) < 0.0
end

@testset "Modified Kling-Gupta" begin
    @test Streamfall.mKGE(test_seq, test_seq) == 1.0
    @test Streamfall.mKGE(test_seq, b_seq) ≈ -0.4142135623730
    @test Streamfall.mKGE(test_seq, c_seq) < 0.0
end

@testset "non-parametric Kling-Gupta" begin
    @test Streamfall.npKGE(test_seq, test_seq) == 1.0
    @test Streamfall.npKGE([10.0, 10.0, 10.0], [10.0, 10.0, 10.0]) == 1.0
    @test Streamfall.npKGE([0.0, 0.0, 0.0], [10.0, 0.0, 0.0]) ≈ -1.386303510546059
    @test Streamfall.npKGE(test_seq, c_seq) < 0.0

    @test !isnan(Streamfall.npKGE(repeat([0], 100), repeat([0.01], 100)))
end


text:
#==============================================================================#
# IAM/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSCore
using AWSIAM
using Test

AWSCore.set_debug_level(1)

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config(region="ap-southeast-2")

#-------------------------------------------------------------------------------
# IAM tests
#-------------------------------------------------------------------------------

@test occursin(r"^arn:aws:iam::[0-9]+:[^:]+$", iam_whoami(aws))

#==============================================================================#
# End of file.
#==============================================================================#


text:
using Revise
module BayesMNIST

using Flux, Flux.Data.MNIST, Statistics, LinearAlgebra, StatsBase
using Flux: onehotbatch, onecold, crossentropy, throttle, @epochs
using Flux: softmax, gradient#, params
using DynamicHMC, LogDensityProblems, TransformVariables, DynamicHMC.Diagnostics, Parameters
using CUDAapi, Random, Distributions, DataFrames, Zygote#, Plots
import ForwardDiff
if has_cuda()
    @info "CUDA is on"
    import CuArrays
    CuArrays.allowscalar(false)
end
using Revise

#load the data
#the purpose is to explore flux in a Bayesian context

Random.seed!(11)

const D=784
const K=10

#Lasso Logistic Categorical regression
struct LLCRegression{TY, TX, Tαλ, Tθλ, Tατ, Tθτ}
  Y::TY
  X::TX

  αλ::Tαλ
  θλ::Tθλ

  ατβ::Tατ
  θτβ::Tθτ

end

###Dimensionality
#=
We have N datapoints, D features, K classes
Y ∈

=#

function (problem::LLCRegression)(θ)
    #unpack parameters and data
    @unpack Y, X, αλ, θλ, ατβ, θτβ = problem
    @unpack βs, τβs, λ²s = θ
    Ncovariates, Npoints = size(X)
    Nclasses = length(λ²s)

    @assert size(βs,1) == Ncovariates
    @assert size(βs,2) == length(τβs) == length(λ²s) == Nclasses

    #p(λ²ᵢ|αλ,θλ)
    logλ²pri = sum([logpdf(Gamma(αλ, θλ), λ²s[i]) for i ∈ 1:Nclasses])

    logτβpri = sum([logpdf(Gamma(ατβ, θτβ), τβs[i]) for i ∈ 1:Nclasses])

    #p(βᵢ|λ²ᵢ, τβᵢ, βᵢ)
    logβpri = sum([logpdf(Laplace(0, 1/τβs[i]), λ²s[i]^0.5*sum((abs).(βs[:,i]))) for i ∈ 1:Nclasses])

    μs = softmax(X' * βs, dims=2)
    #loglik = sum(logpdf(Multinomial(Npoints, μs[n, :])))
    loglik = sum([logpdf(Multinomial(1, μs[n, :]), Y[n]) for n ∈ 1:Npoints])

    return logλ²pri + logτβpri + logβpri + loglik
end

#hyperpriors
const Π = Dict{Symbol, Float64}(
  :αλ=>1.0f0,
  :θλ=>1.0f0,

  :ατβ=>1.0f0,
  :θτβ=>1.0f0
)

function LLCtransform(p::LLCRegression)
    return as((βs=as(Array, size(p.X,1), K),
      τβs=as(Array, asℝ₊, K),
      λ²s=as(Array, asℝ₊, K)))
end

function dynamicmnist(::Type{V} = AbstractVector, ::Type{T}=Float64) where {V,T}
  imgs = MNIST.images()
  labels = MNIST.labels()

  CuArrays.allowscalar(false)

  #make each image into a feature vector
  #concatenate each feature vector into a matrix
  preX =hcat(float.(reshape.(imgs,:))...)

  #add an intercept
  preX = vcat(preX, ones(1,size(preX,2)))[:,1:1000]
  X = preX# |> gpu

  #creates a one-hot matrix of 0 vectors

  Y = Matrix{Float64}(onehotbatch(labels[1:1000], 0:9))

  #error("stop")
  #make the problem
  p = LLCRegression(Y, X, Π[:αλ], Π[:θλ], Π[:ατβ], Π[:θτβ])

  #transform the support
  t = LLCtransform(p)
  P = TransformedLogDensity(t, p)

  #differentiate via AD
  ∇P =  ADgradient(:ForwardDiff, P)

  @info "beginning MCMC"
  testβs = rand(D+1,K)# |> gpu
  testτβs = rand(K)
  testλ²s = rand(K)
  #p((βs=testβs, τβs=testτβs, λ²s=rand(K)))
  results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 10, reporter = LogProgressReport()#,
    #initialization = (ϵ = 0.1, ),
    #warmup_stages = fixed_stepsize_warmup_stages()
    )

  posterior = transform.(t, results.chain)

  β_posterior = first.(transform.(t, results.chain))
  display(mean(β_posterior))

  λ_posterior = last.(transform.(t, results.chain))
  display(mean(λ_posterior) .^ 0.5)

  println("sizeX: $(size(X))")
  println("sizeY: $(size(Y))")

end

end # module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function verify_gcc(gcc)
    try
        return success(`$gcc --version`)
    catch
        return false
    end
end

if Sys.iswindows()
    using WinRPM
end

function build()
    gccpath = ""
    if isfile("deps.jl")
        include("deps.jl")
        if verify_gcc(gcc)
            @info "GCC already installed and package already built"
            return
        else
            rm("deps.jl")
        end
    end

    if haskey(ENV, "CC")
        if !verify_gcc(`$(ENV["CC"]) -v`)
            error("Using compiler override from environment variable CC = $(ENV["CC"]), but unable to run `$(ENV["CC"]) -v`.")
        end
        gccpath = ENV["CC"]
        @info "Using `$gccpath` as C compiler from environment variable CC"
    end

    @info "Looking for GCC"

    gccargs = ``
    if verify_gcc("cc")
        gccpath = "cc"
        @info "Using `cc` as C compiler"
    elseif Sys.iswindows()
        sysroot = joinpath(WinRPM.installdir, "usr", "$(Sys.ARCH)-w64-mingw32", "sys-root")
        gccpath = joinpath(sysroot, "mingw", "bin", "gcc.exe")
        if !isfile(gccpath)
            @info "Could not find GCC, installing using WinRPM"
            WinRPM.install("gcc", yes = true)
        end
        if !isfile(gccpath)
            error("Couldn't install gcc via WinRPM")
        end
        gccargs = `$gccargs --sysroot $sysroot`
        @info "Using `gcc` from WinRPM as C compiler"
    elseif Sys.isunix() && verify_gcc("gcc")
        gccpath = "gcc"
        @info "Using `gcc` as C compiler"
    end

    if isempty(gccpath)
        error("""
            Please make sure to provide a working gcc in your path!
            You may need to install GCC.
            """
        )
        Sys.isapple() && @info "You can install GCC using Homebrew by `brew install gcc`"
        Sys.islinux() && @info """
                                You can install GCC using `sudo apt-get install gcc` on Debian,
                                or use your distro's package manager.
                                """
    end
    open("deps.jl", "w") do io
        println(io, "const gcc = ", repr(`$gccpath $gccargs`))
    end
end

build()


text:
module Examples

# dummy module to appease Documenter

end

max_stars_repo_name: JuliaImageRecon/Examples

max_issues_repo_name: JuliaImageRecon/Examples

max_forks_repo_name: JuliaImageRecon/Examples


text:
using ReservoirComputing, Random

const res_size=20
const ts = 0.:0.1:50.0
const data = sin.(ts)
const train_len = 400
const input_data = reduce(hcat, data[1:train_len-1])
const target_data = reduce(hcat, data[2:train_len])
const predict_len = 100
const test = reduce(hcat, data[train_len+1:train_len+predict_len])
const training_method = StandardRidge(10e-6)

nlas = [NLADefault(), NLAT1(), NLAT2(), NLAT3()]

for n in nlas
    Random.seed!(77)
    esn = ESN(input_data;
        reservoir=RandSparseReservoir(res_size, 1.2, 0.1),
        nla_type = n)
    output_layer = train(esn, target_data, training_method)
    output = esn(Generative(predict_len), output_layer)
    @test maximum(abs.(test .- output)) ./ maximum(abs.(test)) < 0.1
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
module SpWord2Vec
println(@__FILE__)
include("config_global.jl")
include("extend.jl")
include("SpSkipGramA.jl")
include("SpSkipGramB.jl")
include("SpSkipGramC.jl")
include("SpSkipGramD.jl")
include("SpSkipGramE.jl")
include("SpSkipGramBa.jl")
include("SpSkipGramBb.jl")
end # module


text:
using Documenter, Graphene

makedocs(;
    modules=[Graphene],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/chenspc/Graphene.jl/blob/{commit}{path}#L{line}",
    sitename="Graphene.jl",
    authors="Chen Huang",
    assets=String[],
)

deploydocs(;
    repo="github.com/chenspc/Graphene.jl",
)


text:
using Base.Iterators
using DataStructures
using Dates
using JSON

struct DatasetBuildError <: Exception
    error_message::String
end

struct DataSet
    dimensions::OrderedDict
    variables::OrderedDict
    attributes::OrderedDict
    encoding::Dict
end

#  TODO: missing arguments:
#    - grib_errors
#    - index_path
#    - filter_by_keys
function DataSet(
        path::String;
        read_keys::Array{String,1}=String[],
        kwargs...
    )::DataSet

    index_keys = sort([ALL_KEYS..., read_keys...])
    index = FileIndex(path, index_keys)

    return DataSet(
        build_dataset_components(index; read_keys=read_keys, kwargs...)...
    )
end

#  TODO: build_array
struct OnDiskArray
    grib_path::String
    size::Tuple
    offsets::OrderedDict
    message_lengths::Array{Int, 1} #  TODO: Fix seek offset issue
    missing_value::Any
    geo_ndim::Int
    dtype::Type
end

expand_key(key, shape) = Tuple((1:l)[k] for (k, l) in zip(key, shape))

Base.size(A::OnDiskArray) = A.size

Base.axes(A::OnDiskArray) = Tuple(Base.OneTo(i) for i in size(A))
Base.axes(A::OnDiskArray, d::Int) = axes(A)[d]

Base.convert(::Type{T}, A::OnDiskArray) where T <: Array = A[repeat([Colon()], length(size(A)))...]

#  TODO: Use propper `to_indices`, add boundscheck
function Base.getindex(obj::OnDiskArray, key...)
    expanded_keys = expand_key(key, size(obj))
    #  Geograpyh dims (e.g. lat, lon) are on the end and need to be loaded fully
    #  so only look at the other dimensions
    header_items = expanded_keys[1:end-obj.geo_ndim]
    array_field_shape = (
        (length(l) for l in header_items)..., size(obj)[end-obj.geo_ndim+1:end]...
    )
    array_field = Array{obj.dtype}(undef, array_field_shape...)

    geo_ndim_idx = repeat([Colon()], obj.geo_ndim)

    GribFile(obj.grib_path) do file
        message_length_cumsum = cumsum(obj.message_lengths)
        for (header_indexes, offset) in pairs(obj.offsets)
            array_field_indexes = collect(flatten([
                findall(it .== ix)
                for (it, ix)
                in zip(header_items, header_indexes)
            ]))

            if length(array_field_indexes) != length(header_indexes)
                #  If the index (e.g. [10, 4, 1]) is not in the requested header
                #  range (e.g [1:10, 1:4, 2]) then findall will return fewer
                #  items than required (e.g 2 instead of 3). Skip these cases
                continue
            end

            offset_message_index = findfirst(message_length_cumsum .> offset) - 1
            seek(file, offset_message_index)
            message = Message(file)
            values = message["values"]
            array_field[array_field_indexes..., geo_ndim_idx...] = values
        end
    end

    #  TODO: Weird 'correction, not sure if this the right approach. In the case
    #  where the key is like [:, :, *2*, 120, 61] you might have an array field of
    #  shape (10, 4, 1, 120, 61), which correctly means that only the 2nd layer
    #  was loaded. However this means that the index *2* should now be 1
    corrected_key = collect(Any, deepcopy(key))
    corrected_key[collect(array_field_shape) .== 1] .= 1
    #  TODO: Skipped some sections of the python equivalent code as I don't get
    #  what they're for. Should check this out later
    #  TODO: Add in missing value substitution
    return getindex(array_field, corrected_key...)
end

#  TODO: Use parametric struct instead of any
struct Variable
    dimensions::Tuple{Vararg{String}}
    data::Union{Number, Array, OnDiskArray}
    attributes::Dict{String, Any}
end

function Base.:(==)(a::Variable, b::Variable)
    attributes = a.attributes == b.attributes
    dimensions = a.dimensions == b.dimensions
    data = a.data == b.data

    return attributes && dimensions && data
end

#  TODO: Implement filter_by_keys
function enforce_unique_attributes(
        header_values::OrderedDict{String, T} where T <: Array,
        attribute_keys::Array
    )
    attributes = OrderedDict()
    for key in attribute_keys
        values = header_values[key]

        if length(values) > 1
            throw(DatasetBuildError(
                "Attributes are not unique for " *
                "$key: $(values)"
            ))
        end

        value = values[1]

        if !ismissing(value) && !(value in ["undef", "unknown"])
            attributes["GRIB_" * key] = value
        end
    end

    return attributes
end

#  TODO: Implement filter_by_keys
function enforce_unique_attributes(index::FileIndex, attribute_keys::Array)
    attributes = enforce_unique_attributes(
        index.header_values, attribute_keys
    )

    return attributes
end

function build_geography_coordinates(
        index, encode_cf, errors
    )

    first_message = first(index)
    geo_coord_vars = OrderedDict()
    grid_type = cfgrib.getone(index, "gridType")

    if "geography" in encode_cf && grid_type in GRID_TYPES_DIMENSION_COORDS

        column_major = getone(index, "jPointsAreConsecutive") != 0
        #  TODO: column/row major has always confused me, not sure if this
        #  is the correct approach here. Idea is taken from how GRIB.jl
        #  handles reading data with `codes_grib_get_data`:
        #  https://github.com/weech/GRIB.jl/blob/5710a1f462e888ad38f6e3b282df3fb953478d1b/src/message.jl#L355
        if column_major
            geo_dims = ("latitude", "longitude")
            geo_shape = (getone(index, "Ny"), getone(index, "Nx"))
        else
            geo_dims = ("longitude", "latitude")
            geo_shape = (getone(index, "Nx"), getone(index, "Ny"))
        end
        latitudes = first_message["distinctLatitudes"]
        geo_coord_vars["latitude"] = Variable(
            ("latitude",), latitudes, cfgrib.COORD_ATTRS["latitude"]
        )

        if latitudes[1] > latitudes[end]
            geo_coord_vars["latitude"].attributes["stored_direction"] =
                "decreasing"
        end

        geo_coord_vars["longitude"] = Variable(
            ("longitude",),
            first_message["distinctLongitudes"],
            cfgrib.COORD_ATTRS["longitude"],
        )
    elseif "geography" in encode_cf && grid_type in GRID_TYPES_2D_NON_DIMENSION_COORDS
        throw("unimplemented")
    else
        geo_dims = ("values", )
        geo_shape = (getone(index, "numberOfPoints"), )
        try
            latitude = first_message["latitudes"]
            geo_coord_vars["latitude"] = Variable(
                ("values", ),
                latitude,
                COORD_ATTRS["latitude"]
            )

            longitude = first_message["longitudes"]
            geo_coord_vars["longitude"] = Variable(
                ("values", ),
                longitude,
                COORD_ATTRS["longitude"]
            )
        catch e
            rethrow(e)
        end
    end

    return geo_dims, geo_shape, geo_coord_vars
end

function encode_cf_first(
        data_var_attrs::OrderedDict,
        encode_cf::Tuple{Vararg{String}}=("parameter", "time"),
        time_dims::Tuple{Vararg{String}}=("time", "step")
    )

    #  NOTE: marking value as `const` just means it cannot be reassigned, the
    #  value can still be mutated/appended to, so be carfeul `append!`ing to
    #  the constants
    coords_map = deepcopy(cfgrib.ENSEMBLE_KEYS)
    param_id = get(data_var_attrs, "GRIB_paramId", missing)
    data_var_attrs["long_name"] = "original GRIB paramId: $(param_id)"
    data_var_attrs["units"] = "1"

    if "parameter" in encode_cf
        if haskey(data_var_attrs, "GRIB_cfName")
            data_var_attrs["standard_name"] = data_var_attrs["GRIB_cfName"]
        end

        if haskey(data_var_attrs, "GRIB_name")
            data_var_attrs["long_name"] = data_var_attrs["GRIB_name"]
        end

        if haskey(data_var_attrs, "GRIB_units")
            data_var_attrs["units"] = data_var_attrs["GRIB_units"]
        end
    end

    if "time" in encode_cf
        if issubset(time_dims, cfgrib.ALL_REF_TIME_KEYS)
            append!(coords_map, time_dims)
        else
            throw("time_dims $(time_dims) is not a subset of " *
                  "$(cfgrib.ALL_REF_TIME_KEYS)"
            )
        end
    else
        append!(coords_map, cfgrib.DATA_TIME_KEYS)
    end

    append!(coords_map, cfgrib.VERTICAL_KEYS)
    append!(coords_map, cfgrib.SPECTRA_KEYS)

    return coords_map
end

#  TODO: Add filter_by_keys
function build_variable_components(
        index; encode_cf=(),
        errors="warn", squeeze=true, read_keys=String[],
        time_dims=("time", "step")
    )
    data_var_attrs_keys = cfgrib.DATA_ATTRIBUTES_KEYS
    data_var_attrs_keys = [
        data_var_attrs_keys;
        get(cfgrib.GRID_TYPE_MAP, index["gridType"][1], [])
    ]
    data_var_attrs_keys = [data_var_attrs_keys; read_keys]

    data_var_attrs = enforce_unique_attributes(index, data_var_attrs_keys)

    coords_map = encode_cf_first(data_var_attrs, encode_cf, time_dims)

    coord_name_key_map = Dict()
    coord_vars = OrderedDict()

    for coord_key in coords_map
        values = index[coord_key]
        if length(values) == 1 && ismissing(values[1])
            #  TODO: Add logging
            #  @warn "Missing from GRIB Stream $(coord_key)"
            continue
        end

        coord_name = coord_key

        if ("vertical" in encode_cf && coord_key == "level"
                && haskey(data_var_attrs, "GRIB_typeOfLevel"))
            coord_name = data_var_attrs["GRIB_typeOfLevel"]
            coord_name_key_map[coord_name] = coord_key
        end

        attributes = Dict(
            "long_name" => "original GRIB coordinate for key:" *
                           "$(coord_key)($(coord_name))",
            "units"     => "1",
        )

        merge!(attributes, copy(get(cfgrib.COORD_ATTRS, coord_name, Dict())))

        data = sort(
            values,
            rev=get(attributes, "stored_direction", "none") == "decreasing"
        )
        dimensions = (coord_name, )

        if squeeze && length(values) == 1
            data = data[1]
            #  Should single values be in an array as well?
            # typeof(data) == Array ? nothing : data = [data]
            dimensions = ()
        end

        coord_vars[coord_name] = Variable(dimensions, data, attributes)
    end

    header_dimensions = Tuple(
        d for (d, c)
        in pairs(coord_vars)
        if !squeeze || length(c.data) > 1
    )
    #  Loses information on which shape belongs to which dimension
    #  doesn't seem to matter though
    header_shape = Iterators.flatten(
        Tuple(size(coord_vars[d].data) for d in header_dimensions)
    )

    geo_dims, geo_shape, geo_coord_vars = build_geography_coordinates(
        index, encode_cf, errors
    )

    dimensions = (header_dimensions..., geo_dims...)
    shape = (header_shape..., geo_shape...)

    merge!(coord_vars, geo_coord_vars)

    offsets = OrderedDict{NTuple{length(header_dimensions), Int64}, Int}()
    for (header_values, offset) in index.offsets
        header_indexes = Array{Int}(undef, length(header_dimensions))
        for (i, dim) in enumerate(header_dimensions)
            coord_name = get(coord_name_key_map, dim, dim)
            coord_idx = findfirst(index.index_keys .== coord_name)
            header_value = header_values[coord_idx]
            header_indexes[i] = findfirst(coord_vars[dim].data .== header_value)
        end

        offsets[Tuple(header_indexes)] = offset
    end

    data = OnDiskArray(
        index.grib_path,
        shape,
        offsets,
        index.message_lengths,  #  TODO: Fix seek offset issue
        missing,
        length(geo_dims),
        Float32
    )

    if haskey(coord_vars, "time") && haskey(coord_vars, "step")
        # add the 'valid_time' secondary coordinate
        dims, time_data = build_valid_time(
            coord_vars["time"].data,
            coord_vars["step"].data
        )
        attrs = cfgrib.COORD_ATTRS["valid_time"]
        coord_vars["valid_time"] = Variable(dims, time_data, attrs)
    end

    data_var_attrs["coordinates"] = join(keys(coord_vars), " ")
    data_var = Variable(dimensions, data, data_var_attrs)
    dims = OrderedDict(
        (d => s)
        for (d, s)
        in zip(dimensions, size(data_var.data))
    )

    return dims, data_var, coord_vars
end

#  TODO: logging, filter_by_keys
function build_dataset_attributes(index; encoding)
    attributes = enforce_unique_attributes(index, GLOBAL_ATTRIBUTES_KEYS)
    attributes["Conventions"] = "CF-1.7"

    if "GRIB_centreDescription" in keys(attributes)
        attributes["institution"] = attributes["GRIB_centreDescription"]
    end

    attributes_namespace = Dict(
        "cfgrib_version" => cfgrib_jl_version,  # TODO: Package versions are experimental, this should be changed: https://julialang.github.io/Pkg.jl/dev/api/#Pkg.dependencies
        "cfgrib_open_kwargs" => JSON.json(encoding),
        "eccodes_version" => "missing",  # TODO: Not sure how to get this
        "timestamp" => string(Dates.now()),
    )

    history_in = (
        "timestamp GRIB to CDM+CF via " *
        "cfgrib-cfgrib_version/ecCodes-eccodes_version with cfgrib_open_kwargs"
    )

    [history_in=replace(history_in, p) for p in attributes_namespace]
    #  TODO: Fix quotes, should probably still be double quotes not single
    history_in = replace(history_in, "\"" => "'")
    attributes["history"] = history_in

    return attributes
end

  # TODO: Add filter_by_keys
function build_dataset_components(
        index; errors="warn",
        encode_cf=("parameter", "time", "geography", "vertical"),
        squeeze=true, read_keys=String[], time_dims=("time", "step")
    )

    dimensions = OrderedDict()
    variables = OrderedDict()
    for param_id in index["paramId"]
        var_index = filter(index, paramId=param_id)
        dims, data_var, coord_vars = build_variable_components(
            var_index;
            encode_cf=encode_cf,
            errors=errors,
            squeeze=squeeze,
            read_keys=read_keys,
            time_dims=time_dims,
        )

        short_name = get(data_var.attributes, "GRIB_shortName", "paramId$(param_id)")
        var_name = get(data_var.attributes, "GRIB_cfVarName", "unknown")

        if ("parameter" in encode_cf) && !(var_name == "unknown") && !ismissing(var_name)
            short_name = var_name
        end

        merge!(variables, coord_vars)
        merge!(variables, Dict(short_name => data_var))
        merge!(dimensions, dims)
    end

    encoding = Dict(
        "source" => index.grib_path,
        "filter_by_keys" => "not_implemented",  # TODO: Add filter_by_keys
        "encode_cf" => encode_cf
    )
    attributes = build_dataset_attributes(index, encoding=encoding)

    return dimensions, variables, attributes, encoding
end


text:
function barzanRegressTreeLearn(responseVars, features)
	#% learns a compact regression tree
	features=features''
	model = build_tree(vec(responseVars),features,5)
	#model = RegressionTree.fit(features, responseVars, 'MergeLeaves', 'on', 'MinLeaf', 1, 'Prune', 'off', 'MinParent', 10);

	#%save('RegresstionTree_Model', 'model')
	#%load('RegresstionTree_Model', 'model')
	#%disp 'RTREE PAUSE'
	##%pause
	#%model = compact(tree);
	return model
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
function draw(r)
    k = 1
    isExit = !isborder(r, Nord)
    while !isExit
        for _ = 1:k
            move!(r, West)
            if (!isborder(r, Nord)) break end
        end
        if (!isborder(r, Nord)) break end
        k+=1
        for _ = 1:k
            move!(r, Ost)
            if (!isborder(r, Nord)) break end
        end
        if (!isborder(r, Nord)) break end
        k+=1
    end
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
using SymbolicRegression
using Test
using JSON3

X = 2 .* randn(Float32, 2, 1000)
y = 3 * cos.(X[2, :]) + X[1, :] .^ 2 .- 2

options = SymbolicRegression.Options(;
    binary_operators=(+, *, /, -),
    unary_operators=(cos,),
    recorder=true,
    recorder_file="pysr_recorder.json",
    crossoverProbability=0.0,  # required for recording, as not set up to track crossovers.
    npopulations=2,
    npop=100,
    maxsize=20,
    complexity_of_operators=[cos => 2],
)

hall_of_fame = EquationSearch(X, y; niterations=5, options=options, multithreading=true)

data = open(options.recorder_file, "r") do io
    JSON3.read(io)
end

@test haskey(data, :options)
@test haskey(data, :out1_pop1)
@test haskey(data, :out1_pop2)
@test haskey(data, :mutations)

# Test that "Options" is part of the string in `data.options`:
@test contains(data.options, "Options")
@test length(data.mutations) > 1000

# Check whether 10 random elements have the right properties:
for (i, key) in enumerate(keys(data.mutations))
    @test haskey(data.mutations[key], :events)
    @test haskey(data.mutations[key], :score)
    @test haskey(data.mutations[key], :tree)
    @test haskey(data.mutations[key], :loss)
    @test haskey(data.mutations[key], :parent)
    if i > 10
        break
    end
end


text:
using Luxor

demo = Movie(1000, 1000, "comput-graph")

backdrop(scene, framenumber) = background("none")

function frame(scene, framenumber)
    radius = 40
    node_color = "#6ca1f7"

    sethue("black")
    p1 = Point(200, -200)
    p2 = Point(-200, 200)
    temps = Point[between(p1, p2, i) for i in LinRange(0, 1, 4)]
    height = abs(temps[3].y - temps[4].y)
    width = abs(temps[3].x - temps[4].x)
    p_A = Point(temps[2].x, temps[4].y)
    p_x = Point(temps[3].x, temps[4].y + height)
    p_b = Point(temps[1].x, temps[4].y + height)
    p_c = Point(p_b.x + 2 * width, p_b.y)
    p_z = Point(temps[2].x + 2 * width, temps[4].y)

    arrow_configs = (arrowheadlength=13, arrowheadangle=pi/8, linewidth=1.5)
    for (prev_p, next_p) in zip(temps[1:end-1], temps[2:end])
        _, a_p1, _ = intersectionlinecircle(prev_p, next_p, prev_p, radius)
        _, _, a_p2 = intersectionlinecircle(prev_p, next_p, next_p, radius)
        arrow(a_p2, a_p1; arrow_configs...)
    end

    _, _, a_p1 = intersectionlinecircle(p_A, temps[3], temps[3], radius)
    _, a_p2, _ = intersectionlinecircle(p_A, temps[3], p_A, radius)
    arrow(a_p2, a_p1; arrow_configs...)

    _, _, a_p1 = intersectionlinecircle(p_x, temps[4], temps[4], radius)
    _, a_p2, _ = intersectionlinecircle(p_x, temps[4], p_x, radius)
    arrow(a_p2, a_p1; arrow_configs...)

    _, _, a_p1 = intersectionlinecircle(p_x, temps[2], temps[2], radius)
    _, a_p2, _ = intersectionlinecircle(p_x, temps[2], p_x, radius)
    arrow(a_p2, a_p1; arrow_configs...)

    _, _, a_p1 = intersectionlinecircle(p_z, temps[1], temps[1], radius)
    _, a_p2, _ = intersectionlinecircle(p_z, temps[1], p_z, radius)
    arrow(a_p2, a_p1; arrow_configs...)

    _, _, a_p1 = intersectionlinecircle(p_c, temps[1], temps[1], radius)
    _, a_p2, _ = intersectionlinecircle(p_c, temps[1], p_c, radius)
    arrow(a_p2, a_p1; arrow_configs...)

    _, a_p1, _ = intersectionlinecircle(p_x, p_z, p_x, radius)
    _, _, a_p2 = intersectionlinecircle(p_x, p_z, p_z, radius)
    arrow(a_p1, a_p2; arrow_configs...)

    _, a_p1, _ = intersectionlinecircle(p_b, p_z, p_b, radius)
    _, _, a_p2 = intersectionlinecircle(p_b, p_z, p_z, radius)
    arrow(a_p1, a_p2; arrow_configs...)

    # line(p_x, temps[4])
    # line(p_A, temps[3])
    strokepath()

    xaxis = Point[]

    expr = String[
        "f(u) = transpose(u)",
        "f(u, v) = uv",
        "f(M, v) = Mv",
        "f(x1, x2, x3) = x1 + x2 + x3"
    ]

    for (p, ex) in zip(temps, reverse(expr))
        sethue(node_color)
        circle(p, radius, :fill)
        fontsize(20)
        sethue("black")
        text(ex, p - (radius + 5, 0), halign=:right, valign=:middle)
    end

    sethue(node_color)
    circle(p_A, radius, :fill)
    circle(p_x, radius, :fill)
    circle(p_b, radius, :fill)
    circle(p_z, radius, :fill)
    circle(p_c, radius, :fill)
    # circle(xaxis[])

    sethue("black")
    text("f(u, v) = uv", p_z - (radius + 5, 0), halign=:right, valign=:middle)

    fontsize(30)
    configs = (halign=:center, valign=:middle)
    text("A", p_A; configs...)
    text("x", p_x; configs...)
    text("b", p_b; configs...)
    text("c", p_c; configs...)
end

function move_frame(scene, framenumber)
    node_color = "#aaf76c"
    radius = 40
    p1 = Point(200, -200)
    p2 = Point(-200, 200)
    temps = Point[between(p1, p2, i) for i in LinRange(0, 1, 4)]
    height = abs(temps[3].y - temps[4].y)
    width = abs(temps[3].x - temps[4].x)
    p_A = Point(temps[2].x, temps[4].y)
    p_x = Point(temps[3].x, temps[4].y + height)
    p_b = Point(temps[1].x, temps[4].y + height)
    p_c = Point(p_b.x + 2 * width, p_b.y)
    p_z = Point(temps[2].x + 2 * width, temps[4].y)

    sethue(node_color)
    setopacity(0.8)
    if framenumber == 1
        circle(p_x, radius, :fill)
        circle(p_A, radius, :fill)
        circle(p_b, radius, :fill)
        circle(p_c, radius, :fill)
    elseif framenumber == 2
        circle(temps[4], radius, :fill)
        circle(temps[3], radius, :fill)
        circle(p_z, radius, :fill)
    elseif framenumber == 3
        circle(temps[2], radius, :fill)
    elseif framenumber == 4
        circle(temps[1], radius, :fill)
    end

    sethue("black")
    fontsize(30)
    configs = (halign=:center, valign=:middle)
    text("A", p_A; configs...)
    text("x", p_x; configs...)
    text("b", p_b; configs...)
    text("c", p_c; configs...)
end

animate(demo, [
    # Scene(demo, backdrop, 0:200),
    Scene(demo, frame, 0:4)
    Scene(demo, move_frame, 0:4)
    ],
    framerate=2,
    creategif=true,
    pathname="test/comput-graph.gif"
)


text:
using PlutoUI
using Test

# has to be outside of the begin block for julia 1.0 compat
struct Uhm end

@testset "DisplayTricks" begin

    x = "asdf"

    # Show

    @test repr(MIME"derp1"(), Show(MIME"derp1"(), x)) == codeunits(x)

    Base.istextmime(::MIME"derp2") = true

    @test repr(MIME"derp2"(), Show(MIME"derp2"(), codeunits(x))) == x

    s = Show(MIME"cool1"(), x)

    @test Base.showable(MIME"text/plain"(), x)
    @test !Base.showable(MIME"cool1"(), x)
    # everything is plaintext showable
    @test Base.showable(MIME"text/plain"(), s)
    @test Base.showable(MIME"cool1"(), s)

    # AsMIME

    mimer = as_mime(MIME"cool2"())
    am = mimer(s)

    @test Base.showable(MIME"cool1"(), s)
    @test !Base.showable(MIME"cool2"(), s)
    @test Base.showable(MIME"text/plain"(), s)

    @test !Base.showable(MIME"cool1"(), am)
    @test Base.showable(MIME"cool2"(), am)
    @test Base.showable(MIME"text/plain"(), am)

    # WithIOContext

    @test WithIOContext(s, :a => 1, :b => 2) ==
        WithIOContext(s, :a => 1; b = 2) ==
        WithIOContext(s, a = 1; b = 2) ==
        WithIOContext(s; a = 1, b = 2)

    @test repr(MIME"text/plain"(), WithIOContext(1.345234523452, compact=true)) == "1.34523"
    @test repr(MIME"text/plain"(), WithIOContext(1.345234523452, compact=false)) == "1.345234523452"

    m = MIME"hello"()
    Base.istextmime(::MIME"hello") = true

    Base.show(io::IO, ::MIME"hello", ::Uhm) = print(io, get(io, :prop, missing))

    u = Uhm()

    @test repr(m, WithIOContext(u, prop=1)) == "1"
end


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The transportation problem

# Allocation of passenger cars to trains to minimize cars required or car-miles
# run. Based on:
#
# Fourer, D.M. Gay and Brian W. Kernighan, A Modeling Language for Mathematical
# Programming, https://ampl.com/REFS/amplmod.ps.gz Appendix D.
#
# Originally contributed by Louis Luangkesorn, January 30, 2015.

using JuMP
import GLPK
import Test

function example_transp()
	ORIG = ["GARY", "CLEV", "PITT"]
	DEST = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"]
	supply = [1_400, 2_600, 2_900]
	demand = [900, 1_200, 600, 400, 1_700, 1_100, 1_000]
	Test.@test sum(supply) == sum(demand)
	cost = [
		39   14   11   14   16   82    8;
		27    9   12    9   26   95   17;
		24   14   17   13   28   99   20
	]
	model = Model(GLPK.Optimizer)
	@variable(model, trans[1:length(ORIG), 1:length(DEST)] >= 0)
	@objective(
		model,
		Min,
		sum(
			cost[i, j] * trans[i, j]
			for i in 1:length(ORIG), j in 1:length(DEST)
		)
	)
	@constraints(model, begin
		[i in 1:length(ORIG)], sum(trans[i, :]) == supply[i]
		[j in 1:length(DEST)], sum(trans[:, j]) == demand[j]
	end)
	optimize!(model)
	Test.@test termination_status(model) == MOI.OPTIMAL
	Test.@test primal_status(model) == MOI.FEASIBLE_POINT
	Test.@test objective_value(model) == 196200.0
	println("The optimal solution is:")
	println(value.(trans))
	return
end

example_transp()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# using Stars

@testset "Projection string checking" begin
    ga = rast("data/MOD13A2_Henan_2015_2020_10km.tif")
    LON, LAT = st_coords(ga)

    @test size(LON) == size(LAT)
    @test size(LON) == (38, 30)

    LON2, LAT2 = st_coordinates(ga)
    @test maximum(LON - LON2) <= 1e-7
    @test maximum(LAT - LAT2) <= 1e-7

    lon2 = st_coordinates(ga, :x)
    lat2 = st_coordinates(ga, :y)

    lon, lat = st_dim(ga)

    @test maximum(abs.(lon - lon2)) <= 1e-7
    @test maximum(abs.(lat - lat2)) <= 1e-7
end


text:
function fread_note!(S::GphysData, N::Array{Int64,1}, method::String, fmt::String, filestr::String, opts::String)
  src_str = string(timestamp(), " ¦ +source ¦ ", method, "(S, \"", fmt, "\", \"", filestr, "\", ", opts, ")" )
  for i in N
    push!(S.notes[i], src_str)
  end
  return nothing
end

# Five arguments: S, N, method, fname, opts
function fwrite_note!(S::GphysData, i::Int64, method::String, fname::String, opts::String)
  wstr = string(timestamp(), " ¦ write ¦ ", method, "(S", opts, ") ¦ wrote to file ", fname)
  push!(S.notes[i], wstr)
  return nothing
end

function fwrite_note!(C::GphysChannel, method::String, fname::String, opts::String)
  wstr = string(timestamp(), " ¦ write ¦ ", method, "(C", opts, ") ¦ wrote to file ", fname)
  push!(C.notes, wstr)
  return nothing
end

proc_note!(S::GphysData, nn::Array{Int64, 1}, method::String, desc::String) = note!(S, nn, string("processing ¦ ", method, " ¦ ", desc))

proc_note!(S::GphysData, i::Int64, method::String, desc::String) = note!(S, i, string("processing ¦ ", method, " ¦ ", desc))

proc_note!(C::GphysChannel, method::String, desc::String) = note!(C, string("processing ¦ ", method, " ¦ ", desc))


text:
# wrapper functions for measuring velocities with EchelleCCFs
import EchelleCCFs

"""
    calc_ccf(lambdas, intensities, lines, depths, resolution; normalize=true)

Compute the cross correlation function from a spectrum (lambdas and intensities) with a mask computed from a line list (lines).

# Arguments
- `lambdas::AbstractArray{Float64,1}`: List of wavelengths in spectrum.
- `intensities::AbstractArray{Float64,1}`: List of intensities in spectrum.
- `lines::AbstractArray{Float64,1}`: List of line centers for use in CCF mask.
- `depths::AbstractArray{Float64,1}`: List of line depths for use as weights in CCF mask.
- `resolution::Float64`: Spectral resolution of spcectrum.
"""
function calc_ccf(lambdas::AA{T,1}, intensities::AA{T,1}, lines::AA{T,1},
                  depths::AA{T,1}, resolution::T; normalize::Bool=true) where T<:AF
    # make a line list
    line_list = EchelleCCFs.BasicLineList(lines, depths)

    # set mask properties
    speed_of_light = c_ms::T
    mask_width = speed_of_light/resolution
    mask_shape = EchelleCCFs.TopHatCCFMask(mask_width)
    midpoint = 0.0  # just by eye so things aren't wildly shifted
    max_bc = 0.0    #  ~ 2pi AU/year  in m/s
    Δv_step = 400   # m/s arbitrary probably smaller/slower than you need
    Δv_max = 30e3   # m/s arbitrary

    # make ccf_plan
    ccf_plan = EchelleCCFs.BasicCCFPlan(line_list=line_list, mask_shape=mask_shape,
                                        step=Δv_step, max=Δv_max, midpoint=midpoint,
                                        range_no_mask_change=max_bc)

    # calculate the ccf
    v_grid = EchelleCCFs.calc_ccf_v_grid(ccf_plan)
    ccf = EchelleCCFs.ccf_1D(lambdas, intensities, ccf_plan)

    # normalize if normalize==true
    if normalize
        ccf ./= maximum(ccf)
    end
    return v_grid, ccf
end

function calc_ccf(lambdas::AA{T,1}, intensities::AA{T,2}, lines::AA{T,1},
                  depths::AA{T,1}, resolution::T; normalize::Bool=true) where T<:AF
    func = x -> calc_ccf(lambdas, x, lines, depths, resolution, normalize=normalize)
    out = map(func, intensities[:,i] for i in 1:size(intensities, 2))
    return out[1][1], cat([x[2] for x in out]..., dims=2)
end

function calc_ccf(lambdas::AA{T,1}, intensities::AA{T,1}, spec::SpecParams{T}; normalize::Bool=true) where T<:AF
    return calc_ccf(lambdas, intensities, spec.lines, spec.depths, spec.resolution, normalize=normalize)
end

function calc_ccf(lambdas::AA{T,1}, intensities::AA{T,2}, spec::SpecParams{T}; normalize::Bool=true) where T<:AF
    func = x -> calc_ccf(lambdas, x, spec, normalize=normalize)
    out = map(func, intensities[:,i] for i in 1:size(intensities, 2))
    return out[1][1], cat([x[2] for x in out]..., dims=2)
end

"""
    calc_rvs_from_ccf(v_grid, ccf)

Calculate apparent radial velocity from a CCF and velocity grid.

# Arguments
- `v_grid::AbstractArray{Float64,1}`: List of velocities returned by calc_ccf.
- `ccf::AbstractArray{Float64,1}`: CCF values returned by calc_ccf.
"""
function calc_rvs_from_ccf(v_grid::AA{T,1}, ccf::AA{T,1}; frac_of_width_to_fit::T=0.65) where T<:AF
    mrv = EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=frac_of_width_to_fit)
    return mrv(v_grid, ccf)
end

function calc_rvs_from_ccf(v_grid::AA{T,1}, ccf::AA{T,2}; frac_of_width_to_fit::T=0.65) where T<:AF
    func = x -> calc_rvs_from_ccf(v_grid, x, frac_of_width_to_fit=frac_of_width_to_fit)
    out = map(func, ccf[:,i] for i in 1:size(ccf, 2))
    return cat([x[1] for x in out]..., dims=1), cat([x[2] for x in out]..., dims=1)
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
using Polynomials

"""
    poly_chebyshev(n, [::Type{T}=Int, [var=:x]])

Create a Chebyshev polynomial of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_chebyshev(4)
Polynomials.Polynomial(1 - 8*x^2 + 8*x^4)

julia> p2 = poly_chebyshev(4, BigInt)
Polynomials.Polynomial(1 - 8*x^2 + 8*x^4)

julia> p3 = poly_chebyshev(4, Float64, :y)
Polynomials.Polynomial(1.0 - 8.0*y^2 + 8.0*y^4)

```

"""
function poly_chebyshev(n, ::Type{T}=Int, var=:x) where {T<:Number}

    if n==0
        return Polynomial{T}([one(T)], var)
    elseif n==1
        return Polynomial{T}([zero(T), one(T)], var)
    end

    px = Polynomial{T}([zero(T), one(T)], var)
    T0 = Polynomial{T}([one(T)], var)
    T1 = px

    for i = 2:n
        T2 = 2*px*T1 - T0
        T0 = T1
        T1 = T2
    end

    return T1

end

"""
    poly_dchebyshev(n, [::Type{T}=Int, [var=:x]])

Create the derivative of Chebyshev polynomial of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_dchebyshev(4)
Polynomials.Polynomial(-16*x + 32*x^3)

julia> p2 = poly_dchebyshev(4, BigInt)
Polynomials.Polynomial(-16*x + 32*x^3)

julia> p3 = poly_dchebyshev(4, Float64, :y)
Polynomials.Polynomial(-16.0*y + 32.0*y^3)
```
"""
poly_dchebyshev(n, ::Type{T}=Int, var=:x) where {T<:Number} = derivative(poly_chebyshev(n,T, var))

"""
    poly_chebyshev2(n, [::Type{T}=Int, [var=:x]])

Create a Chebyshev of the second kind  of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_chebyshev2(4)
Polynomials.Polynomial(1 - 12*x^2 + 16*x^4)

julia> p2 = poly_chebyshev2(4, BigInt)
Polynomials.Polynomial(1 - 12*x^2 + 16*x^4)

julia> p3 = poly_chebyshev2(4, Float64, :y)
Polynomials.Polynomial(1.0 - 12.0*y^2 + 16.0*y^4)
```

"""
function poly_chebyshev2(n, ::Type{T}=Int, var=:x) where {T<:Number}

    if n==0
        return Polynomial{T}([one(T)], var)
    elseif n==1
        return Polynomial{T}([zero(T), 2*one(T)], var)
    end

    px = Polynomial{T}([zero(T), one(T)], var)
    U0 = Polynomial{T}([one(T)], var)
    U1 = Polynomial{T}([zero(T), 2*one(T)], var)

    for i = 2:n
        U2 = 2*px*U1 - U0
        U0 = U1
        U1 = U2
    end

    return U1

end

"""
    poly_dchebyshev2(n, [::Type{T}=Int, [var=:x]])

Create a derivative of Chebyshev polynomial of the second kind  of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_dchebyshev2(4)
Polynomials.Polynomial(-24*x + 64*x^3)

julia> p2 = poly_dchebyshev2(4, BigInt)
Polynomials.Polynomial(-24*x + 64*x^3)

julia> p3 = poly_dchebyshev2(4, Float64, :y)
Polynomials.Polynomial(-24.0*y + 64.0*y^3)
```

"""
poly_dchebyshev2(n, ::Type{T}=Int, var=:x) where {T<:Number} = derivative(poly_chebyshev2(n,T, var))

"""
    poly_legendre(n, [::Type{T}=Int, [var=:x]])

Create a Legendre polynomial of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_legendre(4)
Polynomials.Polynomial(0.375 - 3.75*x^2 + 4.375*x^4)

julia> p2 = poly_legendre(4, Rational{Int})
Polynomials.Polynomial(3//8 - 15//4*x^2 + 35//8*x^4)

julia> p3 = poly_legendre(4, BigFloat, :y)
Polynomials.Polynomial(0.375 - 3.75*y^2 + 4.375*y^4)
```

"""
function poly_legendre(n, ::Type{T}=Float64, var=:x) where {T<:Number}

    if n==0
        return Polynomial{T}([one(T)], var)
    elseif n==1
        return Polynomial{T}([zero(T), one(T)], var)
    end

    px = Polynomial{T}([zero(T), one(T)], var)
    p0 = Polynomial{T}([one(T)], var)
    p1 = px

    for i = 2:n
        p2 = ( (2i-1)*px*p1 - (i-1)*p0 ) / i
        p0 = p1
        p1 = p2
    end

    return p1
end

"""
    poly_dlegendre(n, [::Type{T}=Int, [var=:x]])

Create a derivative of Legendre polynomial of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_dlegendre(4)
Polynomials.Polynomial(-7.5*x + 17.5*x^3)

julia> p2 = poly_dlegendre(4, Rational{Int})
Polynomials.Polynomial(-15//2*x + 35//2*x^3)

julia> p3 = poly_dlegendre(4, BigFloat, :y)
Polynomials.Polynomial(-7.5*y + 17.5*y^3)
```

"""
poly_dlegendre(n, ::Type{T}=Int, var=:x) where {T<:Number} = derivative(poly_legendre(n,T, var))

"""
    poly_jacobi(n, [::Type{T}=Int, [var=:x]])

Create a Jacobi polynomial of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `a` Weight α of the Jacobi polynomial
 * `b` Weight β of the Jacobi polynomial
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_jacobi(4, 1, 1)
Polynomials.Polynomial(0.625 - 8.75*x^2 + 13.125*x^4)

julia> p2 = poly_jacobi(4, 1, 1, Rational{Int})
Polynomials.Polynomial(5//8 - 35//4*x^2 + 105//8*x^4)

julia> p3 = poly_jacobi(4, 1, 1, BigFloat, :y)
Polynomials.Polynomial(0.625 - 8.75*y^2 + 13.125*y^4)
```
"""
function poly_jacobi(n, a, b, ::Type{T}=Float64, var=:x) where {T<:Number}

    ox = one(T)
    zx = zero(T)
    if n==0
        return Polynomial{T}([one(T)], var)
    elseif n==1
        return Polynomial{T}([(a-b)/2, (a+b+2)/2], var)
    end
    p0 = Polynomial{T}([one(T)], var)
    p1 = Polynomial{T}([(a-b)/2, (a+b+2)/2], var)
    px = Polynomial{T}([zero(T), one(T)], var)

    for i = 1:(n-1)
	a1 = 2*(i+1)*(i+a+b+1)*(2*i+a+b);
	a2 = (2*i+a+b+1)*(a*a-b*b);
	a3 = (2*i+a+b)*(2*i+a+b+1)*(2*i+a+b+2);
	a4 = 2*(i+a)*(i+b)*(2*i+a+b+2);
	p2 = ox/a1*( (a2 + a3*px)*p1 - a4*p0);

        p0 = p1
        p1 = p2
    end

    return p1
end

"""
    poly_djacobi(n, [::Type{T}=Int, [var=:x]])

Create a derivative of Jacobi polynomial of order n `Polynomial` object.

Instead of computing the polynomial at a specific point,
use `Polynomials` package to compute the actual polynomial

Parameters:

 * `n` Order of polynomials
 * `a` Weight α of the Jacobi polynomial
 * `b` Weight β of the Jacobi polynomial
 * `T` Type of polynomial coefficients
 * `var` symbol to be used as variable by the polynomial

# Examples
```julia-repl
julia> p = poly_djacobi(4, 1, 1)
Polynomials.Polynomial(-17.5*x + 52.5*x^3)

julia> p2 = poly_djacobi(4, 1, 1, Rational{Int})
Polynomials.Polynomial(-35//2*x + 105//2*x^3)

julia> p3 = poly_djacobi(4, 1, 1, BigFloat, :y)
Polynomials.Polynomial(-17.5*y + 52.5*y^3)
```
"""
poly_djacobi(n, a, b, ::Type{T}=Float64, var=:x) where {T<:Number} = derivative(poly_jacobi(n,a,b,T,var))


text:
"""
    Encryptor

An `Encryptor` can be used to encrypt a `Plaintext` instance, yielding a `Ciphertext` instance.

See also: [`Plaintext`](@ref), [`Ciphertext`](@ref)
"""
mutable struct Encryptor <: SEALObject
  handle::Ptr{Cvoid}

  function Encryptor(context::SEALContext, public_key::PublicKey, secret_key::SecretKey)
    handleref = Ref{Ptr{Cvoid}}(0)
    retval = ccall((:Encryptor_Create, libsealc), Clong,
                   (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   context, public_key, secret_key, handleref)
    @check_return_value retval
    return Encryptor(handleref[])
  end

  function Encryptor(context::SEALContext, public_key::PublicKey)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Encryptor_Create, libsealc), Clong,
                   (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   context, public_key, C_NULL, handleref)
    @check_return_value retval
    return Encryptor(handleref[])
  end

  function Encryptor(context::SEALContext, secret_key::SecretKey)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Encryptor_Create, libsealc), Clong,
                   (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   context, C_NULL, secret_key, handleref)
    @check_return_value retval
    return Encryptor(handleref[])
  end

  function Encryptor(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::Encryptor)
  if isallocated(object)
    @check_return_value ccall((:Encryptor_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function set_secret_key!(encryptor::Encryptor, secret_key::SecretKey)
  retval = ccall((:Encryptor_SetSecretKey, libsealc), Clong,
                  (Ptr{Cvoid}, Ptr{Cvoid}),
                  encryptor, secret_key)
  @check_return_value retval
  return nothing
end

function encrypt!(destination::Ciphertext, plain::Plaintext, encryptor::Encryptor)
  retval = ccall((:Encryptor_Encrypt, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                 encryptor, plain, destination, C_NULL)
  @check_return_value retval
  return destination
end

function encrypt_symmetric!(destination::Ciphertext, plain::Plaintext, encryptor::Encryptor)
  retval = ccall((:Encryptor_EncryptSymmetric, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, UInt8, Ptr{Cvoid}, Ptr{Cvoid}),
                 encryptor, plain, false, destination, C_NULL)
  @check_return_value retval
  return destination
end
function encrypt_symmetric(plain::Plaintext, encryptor::Encryptor)
  destination = Ciphertext()
  retval = ccall((:Encryptor_EncryptSymmetric, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, UInt8, Ptr{Cvoid}, Ptr{Cvoid}),
                 encryptor, plain, true, destination, C_NULL)
  @check_return_value retval
  return destination
end