text:
using Test

using CompScienceMeshes
using BEAST
using LinearAlgebra

sol = 1.0
G = BEAST.HH3DSingleLayerTDBIO(sol)

S = readmesh(joinpath(@__DIR__, "assets", "sphere2.in"))
S1 = meshrectangle(1.0, 1.0, 0.25, 3)
S2 = CompScienceMeshes.translate(S1, point(1.0,1.0,1.0))

@show numcells(S1)
@show numcells(S2)

X1 = lagrangecxd0(S1)
x1 = refspace(X1)

X2 = lagrangecxd0(S2)
x2 = refspace(X2)

Δt = 0.2
Nt = 300
ΔR = sol*Δt
P = timebasiscxd0(Δt, Nt)
H = timebasisc0d1(Δt, Nt)

Q = convolve(P, H)
q = refspace(Q)

# Overwrite the default quadrature strategy for
# this operator/space combination
τ1 = chart(S1, first(Iterators.drop(cells(S1),0)))
τ2 = chart(S2, first(Iterators.drop(cells(S2),0)))
ι = BEAST.ring(first(BEAST.rings(τ1, τ2, ΔR)), ΔR)

struct DoubleQuadTimeDomainRule end

function momintegrals!(z, op::typeof(G),
    g::typeof(x1), f::typeof(x2), T::typeof(q),
    τ::typeof(τ1), σ::typeof(τ2), ι::typeof(ι),
    qr::DoubleQuadTimeDomainRule)

    qp_τ = quadpoints(g, [τ], (3,))[1]
    qp_σ = quadpoints(f, [σ], (4,))[1]

    for wpv_τ in qp_τ
        x = wpv_τ.point
        gx = wpv_τ.value[1]
        for wpv_σ in qp_σ
            y = wpv_σ.point
            fy = wpv_σ.value[1]

            R = norm(cartesian(x)-cartesian(y))
            first(ι) <= R <= last(ι) || continue
            maxdegree = numfunctions(T)-1
            TR = [(-R)^d for d in 0:maxdegree]
            dxdy = wpv_τ.weight * wpv_σ.weight

            for i in 1 : numfunctions(g)
                for j in 1 : numfunctions(f)
                    for k in 1 : numfunctions(T)
                        z[i,j,k] += dxdy * gx[i] * fy[j] * TR[k] / (4π*R)
end end end end end end

# Compare results for a single monomial
z1 = zeros(numfunctions(x1), numfunctions(x2), numfunctions(q))
for r in BEAST.rings(τ1, τ2, ΔR)
    ι = BEAST.ring(r, ΔR)
    momintegrals!(z1, G, x1, x2, q, τ1, τ2, ι, DoubleQuadTimeDomainRule())
end

qd = quaddata(G, x1, x2, q, [τ1], [τ2], nothing)
z2 = zeros(numfunctions(x1), numfunctions(x2), numfunctions(q))
for r in BEAST.rings(τ1, τ2, ΔR)
    ι = BEAST.ring(r, ΔR)
    quad_rule = quadrule(G, x1, x2, q, 1, τ1, 1, τ2, r, ι, qd)
    BEAST.momintegrals!(z2, G, x1, x2, q, τ1, τ2, ι, quad_rule)
end

@test z1≈z2 rtol=1e-6

# Compare results for a single basis function
CSM = CompScienceMeshes
distance = norm(cartesian(CSM.center(τ1)) - cartesian(CSM.center(τ2)))
k = floor(Int, distance/Δt/sol)

timead = BEAST.temporalassemblydata(Q, kmax=k)


text:
import .ExpressionExplorer: UsingsImports, SymbolsState

"A container for the result of parsing the cell code, with some extra metadata."
Base.@kwdef struct ExprAnalysisCache
    code::String=""
    parsedcode::Expr=Expr(:toplevel, LineNumberNode(1), Expr(:block))
    module_usings_imports::UsingsImports = UsingsImports()
    function_wrapped::Bool=false
    forced_expr_id::Union{PlutoRunner.ObjectID,Nothing}=nothing
end

ExprAnalysisCache(notebook, cell::Cell) = let
    parsedcode=parse_custom(notebook, cell)
    ExprAnalysisCache(
        code=cell.code,
        parsedcode=parsedcode,
        module_usings_imports=ExpressionExplorer.compute_usings_imports(parsedcode),
        function_wrapped=ExpressionExplorer.can_be_function_wrapped(parsedcode),
    )
end

function ExprAnalysisCache(old_cache::ExprAnalysisCache; new_properties...)
    properties = Dict{Symbol,Any}(field => getproperty(old_cache, field) for field in fieldnames(ExprAnalysisCache))
    merge!(properties, Dict{Symbol,Any}(new_properties))
    ExprAnalysisCache(;properties...)
end

"The (information needed to create the) dependency graph of a notebook. Cells are linked by the names of globals that they define and reference. 🕸"
Base.@kwdef struct NotebookTopology
    nodes::ImmutableDefaultDict{Cell,ReactiveNode}=ImmutableDefaultDict{Cell,ReactiveNode}(ReactiveNode)
    codes::ImmutableDefaultDict{Cell,ExprAnalysisCache}=ImmutableDefaultDict{Cell,ExprAnalysisCache}(ExprAnalysisCache)
    cell_order::ImmutableVector{Cell}=ImmutableVector{Cell}()

    unresolved_cells::ImmutableSet{Cell} = ImmutableSet{Cell}()
end

# BIG TODO HERE: CELL ORDER
all_cells(topology::NotebookTopology) = topology.cell_order.c

is_resolved(topology::NotebookTopology) = isempty(topology.unresolved_cells)

function set_unresolved(topology::NotebookTopology, unresolved_cells::Vector{Cell})
    codes = Dict{Cell,ExprAnalysisCache}(
        cell => ExprAnalysisCache(topology.codes[cell]; function_wrapped=false, forced_expr_id=nothing)
        for cell in unresolved_cells
    )
    NotebookTopology(
        nodes=topology.nodes,
        codes=merge(topology.codes, codes),
        unresolved_cells=union(topology.unresolved_cells, unresolved_cells),
        cell_order = topology.cell_order,
    )
end

ext: jl


text:
@testset "FieldVector" begin
    @testset "Immutable Point3D" begin
        eval(quote
            immutable Point3D <: FieldVector{Float64}
                x::Float64
                y::Float64
                z::Float64
            end

            @inline Point3D(xyz::NTuple{3,Float64}) = Point3D(xyz[1], xyz[2], xyz[3])
        end)

        p = Point3D(1.0, 2.0, 3.0)

        @test_throws Exception p[2] = 4.0

        @test size(p) === (3,)
        @test length(p) === 3
        @test eltype(p) === Float64

        @test (p[1], p[2], p[3]) === (p.x, p.y, p.z)

        m = @SMatrix [2.0 0.0 0.0;
                      0.0 2.0 0.0;
                      0.0 0.0 2.0]

        @test m*p === Point3D(2.0, 4.0, 6.0)

        @test similar_type(Point3D) == Point3D
        @test similar_type(Point3D, Float64) == Point3D
        @test similar_type(Point3D, Float32) == SVector{3,Float32}
        @test similar_type(Point3D, (4,)) == SVector{4,Float64}
        @test similar_type(Point3D, Float32, (4,)) == SVector{4,Float32}
    end

    @testset "Mutable Point2D" begin
        eval(quote
            type Point2D{T} <: FieldVector{T}
                x::T
                y::T
            end

            @inline Point2D(xy::NTuple{2}) = Point2D(xy[1], xy[2])
        end)

        p = Point2D(0.0, 0.0)
        p[1] = 1.0
        p[2] = 2.0

        @test size(p) === (2,)
        @test length(p) === 2
        @test eltype(p) === Float64

        @test (p[1], p[2]) === (p.x, p.y)
        @test (p[1], p[2]) === (1.0, 2.0)

        m = @SMatrix [2.0 0.0;
                      0.0 2.0]

        @test (m*p)::Point2D == Point2D(2.0, 4.0)

        @test similar_type(Point2D{Float64}) == Point2D{Float64}
        @test similar_type(Point2D{Float64}, Float32) == Point2D{Float32}
        @test similar_type(Point2D{Float64}, (4,)) == SVector{4,Float64}
        @test similar_type(Point2D{Float64}, Float32, (4,)) == SVector{4,Float32}
    end
end


text:
using Test

@testset "./sr/LocalPolytope/generaotrs.jl" begin

using BellScenario

@testset "generator_vertex()" begin
    @testset "simple 3-2-3 example" begin
        scenario = LocalSignaling(3,3,2)

        D = DeterministicStrategy([0 0 0;1 0 1;0 1 0],scenario)

        generator = LocalPolytope.generator_vertex(D, scenario)

        @test generator == [1 1 0;0 0 1;0 0 0]
        @test generator isa DeterministicStrategy
    end

    @testset "simple 6-4-4 example" begin
        scenario = LocalSignaling(6,4,4)
        D = DeterministicStrategy([1 0 1 1 0 0;0 0 0 0 1 0;0 1 0 0 0 0;0 0 0 0 0 1], scenario)

        generator = LocalPolytope.generator_vertex(D, scenario)

        @test generator == [1 1 1 0 0 0;0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1]
        @test D isa DeterministicStrategy
    end
end

@testset "_sort_perm_groups()" begin
    @test LocalPolytope._sort_perm_groups(
            [(1,1),(3,2),(4,3),(2,4),(5,5)],
            [[1],[2,3],[4,5]]
        ) == [(1,1),(4,3),(3,2),(5,5),(2,4)]
end

@testset "_perm_increase_lexico_score()" begin
    game = [1 1 1 0;0 0 2 0;0 2 0 0;2 0 0 0]
    allowed_col_perms = [collect(1:4)]
    target_row = 1

    canonical_game = LocalPolytope._perm_increase_lexico_score(game, target_row, allowed_col_perms)

    @test canonical_game == [2 0 0 0;1 1 1 0;0 2 0 0;0 0 2 0]
end

@testset "generator_facet()" begin
    @testset "success guessing game" begin
        BG = BellGame([0 1 0 0 0 0;0 0 0 0 1 0;1 0 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 0 1],2)
        scenario = LocalSignaling(6,6,2)

        @test LocalPolytope.generator_facet(BG, scenario) == [1 0 0 0 0 0;0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1]
    end

    @testset "ambiguous guessing games" begin
        BG = BellGame([2 0 0;0 2 0;0 0 2;1 1 1],4)
        scenario = LocalSignaling(3,4,2)

        @test LocalPolytope.generator_facet(BG, scenario) == [2 0 0;1 1 1;0 2 0;0 0 2]

        BG = BellGame([1 1 1 1 0;0 0 3 0 0;0 0 0 3 0;0 1 0 0 2;3 0 0 0 0], 6)
        scenario = LocalSignaling(5,5,2)

        @test LocalPolytope.generator_facet(BG,  scenario) == [3 0 0 0 0;1 1 1 1 0;0 3 0 0 0;0 0 3 0 0;0 0 0 1 2]
    end

    @testset "error guessing games" begin
        BG = BellGame([1 0 0 0;0 1 0 1;0 0 1 1;0 1 1 0], 3)
        scenario = LocalSignaling(4,4,2)

        @test LocalPolytope.generator_facet(BG, scenario) == [1 1 0 0;1 0 1 0;0 1 1 0;0 0 0 1]
    end

    @testset "coarse-grained output" begin
        BG = BellGame([0 1 0;1 0 0;0 0 1;0 1 0],2)
        scenario = LocalSignaling(3,4,2)

        @time @test LocalPolytope.generator_facet(BG, scenario) == [1 0 0;1 0 0;0 1 0;0 0 1]
    end

    @testset "7x7 case raw permutations fails here" begin
        BG = BellGame([
                1 1 1 1 1 1 0 0;0 0 0 0 0 0 0 1;
                0 1 1 1 1 1 1 0;1 0 1 1 1 1 1 0;
                1 1 1 1 0 1 1 0;1 1 0 1 1 1 1 0;
                1 1 1 0 1 1 1 0;1 1 1 1 1 0 1 0;
            ], 8)

        scenario = LocalSignaling(8,8,2)

        @test LocalPolytope.generator_facet(BG, scenario) == [
            1 1 1 1 1 1 0 0;1 1 1 1 1 0 1 0;
            1 1 1 1 0 1 1 0;1 1 1 0 1 1 1 0;
            1 1 0 1 1 1 1 0;1 0 1 1 1 1 1 0;
            0 1 1 1 1 1 1 0;0 0 0 0 0 0 0 1;
        ]
    end

    @testset "testing cases with duplicate columns" begin
        scenario = LocalSignaling(5,5,2)
        BG1 = BellGame([1 0 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 1 0;0 0 0 0 1], 2)
        BG2 = BellGame([1 0 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 1;0 0 0 0 1], 2)
        BG3 = BellGame([0 0 0 0 1;0 0 1 0 0;0 0 1 0 0;0 0 0 1 0;0 1 0 0 0], 2)

        match = [1 0 0 0 0;1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0]

        @test LocalPolytope.generator_facet(BG1, scenario) == match
        @test LocalPolytope.generator_facet(BG2, scenario) == match
        @test LocalPolytope.generator_facet(BG3, scenario) == match
    end

    @testset "test case where values are re-mapped" begin
        BG = BellGame([4 2 2 0 0;
                4 0 0 2 2; 3 0 3 3 0;
                0 3 0 3 0; 0 0 4 0 2
            ], 12)
        scenario = LocalSignaling(5,5,2)

        @test LocalPolytope.generator_facet(BG, scenario) == [4 2 2 0 0;4 0 0 2 2;3 3 0 3 0;0 3 0 0 3;0 0 2 4 0]
    end

    @testset "historic failure no longer fails" begin
        BG = BellGame([1 1 0 0 0; 1 0 1 0 0; 0 1 1 0 0; 0 0 0 0 1; 0 0 0 0 1], 3)
        scenario = LocalSignaling(5,5,2)

        @test LocalPolytope.generator_facet(BG,scenario)  == [1 1 0 0 0; 1 0 1 0 0; 0 1 1 0 0; 0 0 0 1 0; 0 0 0 1 0]
    end
end

end


text:
ispid(pid::AbstractString) = !isempty(readall(pipeline(`ps aux`, `awk '{print $2 }'`)) |> txt -> search(txt, pid))
function cleanpid(file::AbstractString)
	name, ext = splitext(file)
	if ext == ".pid"
		if ispid(name)
			run(`kill $name`)
		end
		rm(file)
	end
end
closeall(_) = map(cleanpid,readdir())
function startit(_)
	closeall(_)
	@async run(`x-terminal-emulator -e 'julia --startup-file=no Spectropolarimetry.jl'`)
end

it = Input(3)
nr = Input(1)
refresh = Input{Any}(leftbutton)
l315c0 = Input{Any}(leftbutton)
l0c0 = Input{Any}(leftbutton)
l45c0 = Input{Any}(leftbutton)
l90c0 = Input{Any}(leftbutton)
l0c315 = Input{Any}(leftbutton)
l0c45 = Input{Any}(leftbutton)
dark = Input{Any}(leftbutton)
quit = Input{Any}(leftbutton)

lift(startit, refresh)

consume(it, typ=Any, init=empty) do i
	open("integrationtime","w") do o
		print(o,1000i)
	end
end
consume(nr, typ=Any, init=empty) do i
	open("nreplicates","w") do o
		print(o,i)
	end
end

consume(l315c0, typ=Any, init=empty) do i
	open("rotation","w") do o
		print(o,"l315c0")
	end
end
consume(l0c0, typ=Any, init=empty) do i
	open("rotation","w") do o
		print(o,"l0c0")
	end
end
consume(l45c0, typ=Any, init=empty) do i
	open("rotation","w") do o
		print(o,"l45c0")
	end
end
consume(l90c0, typ=Any, init=empty) do i
	open("rotation","w") do o
		print(o,"l90c0")
	end
end
consume(l0c315, typ=Any, init=empty) do i
	open("rotation","w") do o
		print(o,"l0c315")
	end
end
consume(l0c45, typ=Any, init=empty) do i
	open("rotation","w") do o
		print(o,"l0c45")
	end
end
consume(dark, typ=Any, init=empty) do i
	open("rotation","w") do o
		print(o,"dark")
	end
end
lift(closeall, quit, typ=Any, init=empty)

function main(window)
	push!(window.assets, "widgets")

	vbox(
	hbox("Number of scans to average", slider(1:100, value=1) >>> nr),
	hbox("Integration time (milli seconds)", flex(slider(3:60_000, value=3) >>> it)),
	hbox(button("l315c0") >>> l315c0 |> fillcolor("#eeb"), button("l0c0") >>> l0c0 |> fillcolor("#eeb"), button("l45c0") >>> l45c0 |> fillcolor("#eeb"), button("l90c0") >>> l90c0 |> fillcolor("#eeb"), button("l0c315") >>> l0c315 |> fillcolor("#859"), button("l0c45") >>> l0c45 |> fillcolor("#859"), button("dark") >>> dark |> fillcolor("#875"), button("Quit") >>> quit |> fillcolor("#f75"),button("Refresh") >>> refresh |> fillcolor("#f75"))
	)
end


text:
testfile = ARGS[1]

expandtabs(str) = replace(str, r"\t", "  ")
totab4(str) = " " ^ (div(length(expandtabs(str)), 2) * 4)

function transform_line(line::AbstractString)
    m = match(r"^(\s*)facts\(\"(.*)\"\)\s+do(\s*)\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@testset \"$(m.captures[2])\" begin\n"
    end

    m = match(r"^(\s*)context\(\"(.*)\"\)\s+do\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@testset \"$(m.captures[2])\" begin\n"
    end

    m = match(r"^(\s*)@fact_throws\s+(.*)\s+([A-Za-z]+)\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test_throws $(m.captures[3]) $(m.captures[2])\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+is([a-z]+)\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test is$(m.captures[3])($(m.captures[2]))\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+not\(is([a-z]+)\)\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test !is$(m.captures[3])($(m.captures[2]))\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+less_than_or_equal\((.*)\)\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test $(m.captures[2]) <= $(m.captures[3])\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+greater_than_or_equal\((.*)\)\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test $(m.captures[2]) >= $(m.captures[3])\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+less_than\((.*)\)\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test $(m.captures[2]) < $(m.captures[3])\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+greater_than\((.*)\)\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test $(m.captures[2]) > $(m.captures[3])\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+true\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test $(m.captures[2])\n"
    end

    m = match(r"^(\s*)@fact\s+(.*)\s+-->\s+(.*)\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))@test $(m.captures[2]) == $(m.captures[3])\n"
    end

    m = match(r"^(\s*)(.*)\s*\n", line)
    if m != nothing
        return "$(totab4(m.captures[1]))$(m.captures[2])\n"
    end

    line
end

newlines = open(testfile, "r") do fh
    map(transform_line, readlines(fh))
end

# Delete so only max one empty line at end
while length(newlines) >= 2 && strip(newlines[end]) == "" && strip(newlines[end-1]) == ""
    deleteat!(newlines, length(newlines))
end

altname = tempname()
open(altname, "w") do fh
    print(fh, join(newlines, ""))
end

cmd = `diff $(testfile) $(altname)`
try
    run(cmd)
catch _e
end

# Remove file unless we should write it
if length(ARGS) < 2
    rm(altname)
elseif length(ARGS) >= 2 && in(lowercase(ARGS[2]), ["o", "-o", "--overwrite", "overwrite"])
    mv(altname, testfile; remove_destination = true)
else length(ARGS) >= 2 && in(lowercase(ARGS[2]), ["s", "-s", "--save", "write"])
    println("Altered file saved to: ", altname)
end


text:
using Distributions, StatsBase

function imcmc_multi_insert_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::SisMcmcInsertDeleteCenter,
    ind::AbstractVector{Int},
    V::UnitRange,
    K_in_ub::Int
    )

    prop_pointers = mcmc.prop_pointers
    ν_td = mcmc.ν_td
    N = length(S_curr)
    len_dist = mcmc.len_dist

    # Make insertion Categorical distribution ...
    # p_ins =
    mean_len = floor(Int,mean(length.(S)))
    log_ratio = 0.0
    for i in ind
        tmp = popfirst!(prop_pointers)
        m = rand(len_dist)
        resize!(tmp, m)
        sample!(p_ins, tmp)
        insert!(S_prop, i, tmp)
        log_ratio += - logpdf(len_dist, m) + m*log(length(V)) - Inf * (m > K_in_ub)
    end
    log_ratio += log(ν_td) - log(min(ν_td,N))
    return log_ratio

end


text:
module SeqEvol

export evolMSA, extract_params

using DelimitedFiles
using FastaIO
using LinearAlgebra
using StatsBase
using GZip

include("read_write.jl")
include("evol.jl")

end # module

max_stars_repo_name: matteobisardi/SeqEvol

max_issues_repo_name: matteobisardi/SeqEvol

max_forks_repo_name: matteobisardi/SeqEvol


text:
using Plots, Distributions

plot(
    Uniform(0, 1),
    label="Uniform(0, 1)",
    xlabel="Success Probability q",
    ylabel="Probability Density",
    xlims=(0, 1),
    ylims=(0, 2.8)
)
plot!(
    Beta(4, 6),
    label="Beta(α=4, β=6)",
    linestyle=:dash
)


text:
using Shamir, Polynomials, Test

@info "Testing creation of shares"
prod_coeffs = [1234, 166, 94]
n = 6
prod_poly = Poly(prod_coeffs)
shares = Shamir.construct_shares(n, prod_poly)
@test shares == [[1, 1494], [2, 1942], [3, 2578], [4, 3402], [5, 4414], [6, 5614]]

@info "Testing Recover secret"
prod_coeffs = [1234, 166, 94]
n = 6 #total number of parties
k = 3 #min num of shares
p = 1613 #field
prod_poly = Poly(prod_coeffs)
shares = Shamir.construct_shares(n, prod_poly)
secret = Shamir.recover_secret(shares, n, k, p)
@test secret == 1234.0

@info "All tests completed"


text:
using ArgoData
using Test

@testset "ArgoData.jl" begin
    # Write your own tests here.
end


text:
module miniWeather

using ArgParse
using NetCDF
using ProgressMeter
import Base.Threads.@spawn
#using LoopVectorization: @avx

include("const.jl")
include("utils.jl")
include("Initialize.jl")
include("Timestep.jl")

export Grid, Model, init, run!, reductions

end


text:
function posteriortemperature(lpdf::LineDensity, Trange::AbstractRange{Float64}, Srange::AbstractRange{Float64})

    @showprogress map(P -> logl(lpdf, P[2] * observedunitgalaxyvector(P[1])), Iterators.product(Trange, Srange))

end


text:
# %% imports
using Random
using LinearAlgebra
using Statistics
using Plots
using YAML
# %% classes

cfg = YAML.load(open("/Users/titou/Documents/PFE/RL-STDP/Julia/CartPole/CMAES/cfg_cmaes.yml"))

mutable struct CMAES
    N::Int
    μ::Int
    λ::Int
    τ::Float64
    τ_c::Float64
    τ_σ::Float64
    population::Array{Array{Float64}}
    offspring::Array{Array{Float64}}
    F_μ::Array{Float64}
    F_λ::Array{Float64}
    C::Array{Float64}
    s::Array{Float64}
    s_σ::Array{Float64}
    σ::Float64
    E::Array{Float64}
    W::Array{Float64}
    x::Array{Float64}
    cfg::Dict
end

function CMAES(init::Array{Float64} ;N=2, μ=1, λ=10, τ=sqrt(N), τ_c=N^2, τ_σ=sqrt(N), cfg::Dict = cfg)
    x = init
    x = clamp.(x, -1.0, 1.0)
    population = fill(x, µ)
    offspring = Array{Array{Float64}}(undef, λ)
    F_µ = Inf .* ones(µ)
    F_λ = Inf .* ones(λ)
    C = Array(Diagonal{Float64}(I, N))
    s = zeros(N)
    s_σ = zeros(N)
    σ = 1.0
    E = zeros(N, λ)
    W = zeros(N, λ);
    CMAES(N, μ, λ, τ, τ_c, τ_σ, population, offspring, F_µ, F_λ, C, s, s_σ, σ, E, W, x, cfg)
end

function step!(c::CMAES; obj=objective, visualize=false, anim=Nothing)
    # L1
    rng_offspring = Random.MersenneTwister(1234)
    sqrt_c = cholesky((c.C + c.C') / 2.0).U
    for i in 1:c.λ
        c.E[:,i] = randn(rng_offspring, c.N)
        c.W[:,i] = c.σ * (sqrt_c * c.E[:,i])
        c.offspring[i] = c.x + c.W[:,i]
        c.offspring[i] = clamp.(c.offspring[i], -1.0, 1.0)
        c.F_λ[i] = obj(c.offspring[i])
    end
    # Select new parent population
    idx = sortperm(c.F_λ)[1:c.μ]
    for i in 1:c.μ
        c.population[i] = c.offspring[idx[i]]
        c.F_μ[i] = c.F_λ[idx[i]]
    end
    # L2
    w = vec(mean(c.W[:,idx], dims=2))
    c.x += w
    # L3
    c.s = (1.0 - 1.0/c.τ)*c.s + (sqrt(c.μ/c.τ * (2.0 - 1.0/c.τ))/c.σ)*w
    # L4
    c.C = (1.0 - 1.0/c.τ_c).*c.C + (c.s./c.τ_c)*c.s'
    # L5
    ɛ = vec(mean(c.E[:,idx], dims=2))
    c.s_σ = (1.0 - 1.0/c.τ_σ)*c.s_σ + sqrt(c.μ/c.τ_σ*(2.0 - 1.0/c.τ_σ))*ɛ
    # L6
    c.σ = c.σ*exp(((c.s_σ'*c.s_σ)[1] - c.N)/(2*c.N*sqrt(c.N)))
    if visualize
        plot(xs, ys, fz, st=:contour)
        scatter!([c.offspring[i][1] for i in 1:c.λ], [c.offspring[i][2] for i in 1:c.λ],
            xlims=(-5, 5), ylims=(-5, 5), legend=:none)
        scatter!([c.x[1]], [c.x[2]], color=:black, marker=:rect,
            xlims=(-5, 5), ylims=(-5, 5), legend=:none)
        frame(anim)
    end
    c
end

function plot_obj()
    c = CMAES()
    println("x initial: ", c.x)
    anim = Animation()
    for i in 1:100
        v = mod(i, 1) == 0
        step!(c, visualize=v, anim=anim)
    end
    println("x final: ", c.x)
    gif(anim)
end

# %% test

# solution = [3.5, -0.2]
# sphere(x::Array{Float64}) = sum((x .- solution).^2)
#
# xs = -5.0:0.1:5.0
# ys = -5.0:0.1:5.0
#
# objective = sphere # sphere, himmelblau, styblinski_tang, rastrigin
# fz(x, y) = objective([x, y])
# println(solution) # optimal for sphere and rastrigin
# plot_obj()

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
# __precompile__()

module Sundials

import Reexport
Reexport.@reexport using DiffEqBase
import DataStructures
import Logging
import DiffEqBase
import SparseArrays
import LinearAlgebra

import Libdl
using CEnum

const warnkeywords = (
    :d_discontinuities,
    :isoutofdomain,
    :unstable_check,
    :calck,
    :internalnorm,
    :gamma,
    :beta1,
    :beta2,
    :qmax,
    :qmin,
    :qoldinit,
)

warnlist = Set(warnkeywords)
warnida = union(warnlist, Set((:dtmin,)))

using Sundials_jll

export solve,
    SundialsODEAlgorithm, SundialsDAEAlgorithm, ARKODE, CVODE_BDF, CVODE_Adams, IDA, KINSOL

# some definitions from the system C headers wrapped into the types_and_consts.jl
const DBL_MAX = prevfloat(Inf)
const DBL_MIN = nextfloat(-Inf)
const DBL_EPSILON = eps(Cdouble)
const int64_t = Int64
const MPI_DOUBLE = nothing
const MPI_INT64_T = nothing
macro SUNDIALS_F77_FUNC(name, NAME)
    Symbol(string(name) * "_64_")
end

# Issue 328: these symbols exist in other libraries
const libsundials_sundials = libsundials_cvode
const libsundials_sunlinsol = libsundials_cvode
const libsundials_sunnonlinsol = libsundials_cvode
# sunmatrix has been renamed to sunmatrix[dense/sparse/band]
# const libsundials_sunmatrix = libsundials_cvode

const SPARSE_SOLVERS = (:KLU,)

include("../lib/libsundials_common.jl")
include("types_and_consts_additions.jl")

include("handle.jl")
include("nvector_wrapper.jl")

include("../lib/libsundials_api.jl")

for ff in names(@__MODULE__; all=true)
    fname = string(ff)
    if occursin("SetLinearSolver", fname) &&
        !occursin("#", fname) && # filter out compiler generated names
        !occursin("Dls", fname) && !occursin("Spils", fname) # filter out old names
        @eval $ff(mem, LS::SUNLinearSolver, A::Ptr, args...) = $ff(mem, LS, convert(SUNMatrix, A), args...)
    end
end

include("simple.jl")
include("common_interface/function_types.jl")
include("common_interface/verbosity.jl")
include("common_interface/algorithms.jl")
include("common_interface/integrator_types.jl")
include("common_interface/integrator_utils.jl")
include("common_interface/solve.jl")

##################################################################
# Deprecations
##################################################################

end # module


text:
using EngEconomics, Roots, Plots

# Given
aInit = -9000
aAnnual = 5600
aScrap = 1000
aLife = 2

bInit = -5700
bAnnual = 2100
bScrap = 2700
bLife = 3

cInit = -4800
cAnnual = 1200
cScrap = 0
# cLife = infinity

MARR = 0.05

# Find
aAW = aInit * capitalRecoveryFactor(MARR, aLife) + aAnnual + aScrap * sinkingFundFactor(MARR, aLife)
bAW = bInit * capitalRecoveryFactor(MARR, bLife) + bAnnual + bScrap * sinkingFundFactor(MARR, bLife)
cAW = cInit * MARR + cAnnual


text:
module ACEbonds

import ACE
import ACE: B1pMultiplier, AbstractState

include("envelopes.jl")

include("bselectors.jl")

include("cylindrical.jl")

include("iterator.jl")

include("bondpot.jl")

end


text:
module Feature
    export @activate, is_activated
    separate_module_features = Dict()
    function is_activated(symbol, mod)
        get(separate_module_features, symbol) do
            nothing
        end |>
        function (modules)
            if modules === nothing
                false
            else
                mod in modules
            end
        end
    end
    macro activate(symbols...)
        if :TypeLevel in symbols
            modules = get(separate_module_features, :TypeLevel) do
                    separate_module_features[:TypeLevel] = []
            end
            if !(__module__ in modules)
                push!(modules, __module__)
            end
        end
    end
end


text:
module TestStacking

using Test
using MLJBase
using MLJModelInterface
using ..Models
using Random
using StableRNGs

import Distributions

rng = StableRNGs.StableRNG(1234)

function model_evaluation(models::NamedTuple, X, y; measure=rmse)
    cv = CV(;nfolds=3)
    results = []
    for model in models
        mach = machine(model, X, y)
        ev = evaluate!(mach; resampling=cv, verbosity=0, measure=measure, check_measure=false)
        push!(results, ev.measurement[1])
    end
    results
end

function test_internal_evaluation(internalreport, std_evaluation, modelnames)
    for model in modelnames
        model_ev = internalreport[model]
        std_ev = std_evaluation[model]
        @test model_ev isa PerformanceEvaluation
        @test model_ev.per_fold == std_ev.per_fold
        @test model_ev.measurement == std_ev.measurement
        @test model_ev.per_observation[1] === std_ev.per_observation[1] === missing
        @test model_ev.per_observation[2] == std_ev.per_observation[2]
        @test model_ev.operation == std_ev.operation
        @test model_ev.report_per_fold == std_ev.report_per_fold
        @test model_ev.train_test_rows == std_ev.train_test_rows
    end
end

@testset "Testing Stack on Continuous target" begin
    X, y = make_regression(500, 5; rng=rng)

    @testset "Testing Deterministic Stack" begin
    # Testing performance

    # The dataset is a simple regression model with intercept
    # No model in the stack can recover the true model on its own
    # Indeed, FooBarRegressor has no intercept
    # By combining models, the stack can generalize better than any submodel
    # And optimize the rmse

    models = (constant=DeterministicConstantRegressor(),
                decisiontree=DecisionTreeRegressor(),
                ridge_lambda=FooBarRegressor(;lambda=0.1),
                ridge=FooBarRegressor(;lambda=0))

    mystack = Stack(;metalearner=FooBarRegressor(),
                    resampling=CV(;nfolds=3),
                    models...)

    results = model_evaluation((stack=mystack, models...), X, y)
    @test argmin(results) == 1

    # Mixing ProbabilisticModels amd Deterministic models as members of the stack
    models = (constant=ConstantRegressor(),
                decisiontree=DecisionTreeRegressor(),
                ridge_lambda=FooBarRegressor(;lambda=0.1),
                ridge=FooBarRegressor(;lambda=0))

    mystack = Stack(;metalearner=FooBarRegressor(),
                    resampling=CV(;nfolds=3),
                    models...)
    # Testing attribute access of the stack
    @test propertynames(mystack) == (:metalearner, :resampling, :measures, :cache, :acceleration,
        :constant, :decisiontree, :ridge_lambda, :ridge)

    @test mystack.decisiontree isa DecisionTreeRegressor

    @test target_scitype(mystack) == target_scitype(FooBarRegressor())
    @test input_scitype(mystack) == input_scitype(FooBarRegressor())

    # Testing fitted_params results are easily accessible for each
    # submodel. They are in order of the cross validation procedure.
    # Here 3-folds then 3 machines + the final fit
    mach = machine(mystack, X, y)
    fit!(mach, verbosity=0)
    fp = fitted_params(mach)
    @test nrows(getfield(fp, :constant)) == 4
    @test nrows(getfield(fp, :decisiontree)) == 4
    @test nrows(getfield(fp, :ridge)) == 4
    @test nrows(getfield(fp, :ridge_lambda)) == 4

    # The metalearner has been fit and has one coefficient
    # for each model in the library (No intercept)
    @test fp.metalearner isa Vector{Float64}
    @test nrows(fp.metalearner) == 4

    # Testing prediction is Deterministic
    @test predict(mach) isa Vector{Float64}
    end

    @testset "Testing ProbabilisticStack" begin
        models = (constant=ConstantRegressor(),
                    decisiontree=DecisionTreeRegressor(),
                    ridge_lambda=FooBarRegressor(;lambda=0.1),
                    ridge=FooBarRegressor(;lambda=0))

        # The type of the stack is determined by the type of the metalearner
        metalearner = ConstantRegressor(;distribution_type=Distributions.Cauchy)
        mystack = Stack(;metalearner=metalearner,
                    resampling=CV(;nfolds=3),
                    models...)

        @test target_scitype(mystack) == target_scitype(metalearner)
        @test input_scitype(mystack) == input_scitype(FooBarRegressor())

        mach = machine(mystack, X, y)
        fit!(mach, verbosity=0)
        @test predict(mach) isa Vector{Distributions.Cauchy{Float64}}

    end

end

@testset "Testing ProbabilisticStack on Finite target" begin
    X, y = make_blobs(;rng=rng, shuffle=false)

    models = (constant=ConstantClassifier(),
                decisiontree=DecisionTreeClassifier(),
                knn=KNNClassifier())

    mystack = Stack(;metalearner=DecisionTreeClassifier(),
                    resampling=StratifiedCV(;nfolds=3),
                    models...)

    # Check input and target scitypes
    @test target_scitype(mystack) == target_scitype(DecisionTreeClassifier())
    # Here the greatest lower bound is the scitype of the knn
    @test input_scitype(mystack) == input_scitype(KNNClassifier())

    mach = machine(mystack, X, y)
    fit!(mach, verbosity=0)
    @test predict(mach) isa Vector{<:MLJBase.UnivariateFinite}

end

@testset "Stack constructor valid argument checks" begin
    # metalearner should have target_scitype:
    # Union{AbstractArray{<:Continuous}, AbstractArray{<:Finite}}
    @test_throws ArgumentError Stack(;metalearner=Standardizer(),
                        constant=ConstantClassifier())

    @test_throws ArgumentError Stack(;constant=KNNRegressor())
end

@testset "Misc" begin
    # Test setproperty! behaviour
    models = (constant=DeterministicConstantRegressor(),
                decisiontree=DecisionTreeRegressor(),
                ridge_lambda=FooBarRegressor(;lambda=0.1))

    mystack = Stack(;metalearner=FooBarRegressor(),
                    measures=rmse,
                    resampling=CV(;nfolds=3),
                    models...)
    @test mystack.ridge_lambda.lambda == 0.1
    @test mystack.metalearner isa FooBarRegressor
    @test mystack.resampling isa CV

    mystack.ridge_lambda = FooBarRegressor(;lambda=0.2)
    @test mystack.ridge_lambda.lambda == 0.2

    mystack.metalearner = DecisionTreeRegressor()
    @test mystack.metalearner isa DecisionTreeRegressor

    mystack.resampling = StratifiedCV()
    @test mystack.resampling isa StratifiedCV

    # Test measures accepts a single measure
    @test mystack.measures == [rmse,]

    # using inner constructor accepts :resampling and :metalearner
    # as modelnames
    modelnames = (:resampling, :metalearner)
    models = [DeterministicConstantRegressor(), FooBarRegressor(;lambda=0)]
    metalearner = DeterministicConstantRegressor()
    resampling = CV()
    cache = true
    acceleration = CPU1()

    MLJBase.DeterministicStack(modelnames, models, metalearner, resampling, nothing, cache, acceleration)

    # Test input_target_scitypes with non matching target_scitypes
    models = [KNNRegressor()]
    metalearner = KNNClassifier()
    inp_scitype, tg_scitype = MLJBase.input_target_scitypes(models, metalearner)
    @test tg_scitype == Unknown
    @test inp_scitype == Table{<:AbstractVector{<:Continuous}}

    # Test input_target_scitypes with non matching target_scitypes
    models = [ConstantClassifier(), DecisionTreeClassifier()]
    metalearner = KNNClassifier()
    inp_scitype, tg_scitype = MLJBase.input_target_scitypes(models, metalearner)
    @test tg_scitype == AbstractVector{<:Finite}
    @test inp_scitype == Table{<:Union{AbstractVector{<:Continuous},
                                       AbstractVector{<:Count},
                                       AbstractVector{<:OrderedFactor}}}

    # Changing a model to a non compatible target scitype
    initial_stack = Stack(;metalearner=FooBarRegressor(),
                            resampling=CV(;nfolds=3),
                            constant = DeterministicConstantRegressor(),
                            fb=FooBarRegressor())
    initial_stack.constant = ConstantClassifier()
    @test_throws DomainError clean!(initial_stack)

    # Test check_stack_measures with
    # probabilistic measure and deterministic model
    measures =[log_loss]
    stack = Stack(;metalearner=FooBarRegressor(),
                    resampling=CV(;nfolds=3),
                    measure=measures,
                    constant=ConstantRegressor(),
                    fb=FooBarRegressor())
    X, y = make_regression()

    @test_logs((:error, r"Problem fitting"),
               (:info, r"Running type"),
               (:info, r"Type checks okay"),
               @test_throws ArgumentError fit!(machine(stack, X, y), verbosity=0))

    @test_throws ArgumentError MLJBase.check_stack_measures(stack, 0, measures, y)

    # This will not raise
    stack.measures = nothing
    fit!(machine(stack, X, y), verbosity=0)
end

@testset  "function oos_set" begin
    X = (x = Float64[1, 1, 2, 2, 3, 3],)
    y = coerce(['a', 'b', 'b', 'c', 'a', 'a'], Multiclass)
    n = nrows(y)

    model1 = KNNClassifier(K=2)
    model2 = ConstantClassifier()
    judge = KNNClassifier(K=3)

    stack = Stack(metalearner=judge,
                model1=model1,
                model2=model2,
                resampling=CV(;nfolds=3, shuffle=true, rng=StableRNG(123)))

    Xs = source(X)
    ys = source(y)

    ttp = MLJBase.train_test_pairs(stack.resampling, 1:n, X, y)
    Zval, yval, folds_evaluations = MLJBase.oos_set(stack, Xs, ys, ttp)

    # No internal measure has been provided so the resulting
    # folds_evaluations contain nothing
    @test all(x === nothing for x in folds_evaluations)

    # To be accessed, the machines need to be trained
    fit!(Zval, verbosity=0)
    # Each model in the library should output a 3-dim vector to be concatenated
    # resulting in the table of shape (nrows, 6) here nrows=6
    # for future use by the metalearner
    sc = schema(Zval())
    @test nrows(Zval()) == 6
    @test sc.names == (:x1, :x2, :x3, :x4, :x5, :x6)

    # The lines of yval should match the reordering indexes
    # of the original y (reordering given by the folds node)
    reordering = vcat([x[2] for x in ttp]...)
    @test yval() == y[reordering]
    # And the same is true for Zval, let's check this for model1's output
    # on the first fold, ie (2 first rows, 3 first columns)
    # First we need to train the model
    trainingrows = ttp[1][1]
    Xtrain = selectrows(X, trainingrows)
    ytrain = selectrows(y, trainingrows)
    mach = machine(model1, Xtrain, ytrain)
    fit!(mach, verbosity=0)

    # Then predict on the validation rows
    Xpred = selectrows(X, ttp[1][2])
    Zval_expected_dist = predict(mach, Xpred)
    # This is a distribution, we need to apply the appropriate transformation
    Zval_expected = pdf(Zval_expected_dist, levels(first(Zval_expected_dist)))
    @test matrix(Zval())[1:2, 1:3] == Zval_expected

end

@testset "An integration test for stacked classification" begin

    # We train a stack by hand and compare with the canned version
    # `Stack(...)`. There are two base learners, with 3-fold
    # cross-validation used to construct the out-of-sample base model
    # predictions.

    probs(y) = pdf(y, levels(first(y)))

    # data:
    N = 200
    X = (x = rand(rng, 3N), )
    y = coerce(rand("abc", 3N), Multiclass)

    # row splits:
    test1 = 1:N
    test2 = (N + 1):2N
    test3 = (2N + 1):3N
    train1 = (N + 1):3N
    train2 = vcat(1:N, (2N + 1):3N)
    train3 = 1:2N

    # base `model1`:
    model1 = KNNClassifier(K=2)
    mach1 = machine(model1, X, y)
    fit!(mach1, rows=train1, verbosity=0)
    y11 = predict(mach1, rows=test1) |> probs
    mach1 = machine(model1, X, y)
    fit!(mach1, rows=train2, verbosity=0)
    y12 = predict(mach1, rows=test2) |> probs
    mach1 = machine(model1, X, y)
    fit!(mach1, rows=train3, verbosity=0)
    y13 = predict(mach1, rows=test3) |> probs
    y1_oos = vcat(y11, y12, y13)
    mach1_full = machine(model1, X, y)
    fit!(mach1_full, verbosity=0)
    y1 = predict(mach1_full, X) |> probs

    # base `model2`:
    model2 = DecisionTreeClassifier()
    mach2 = machine(model2, X, y)
    fit!(mach2, rows=train1, verbosity=0)
    y21 = predict(mach2, rows=test1) |> probs
    mach2 = machine(model2, X, y)
    fit!(mach2, rows=train2, verbosity=0)
    y22 = predict(mach2, rows=test2) |> probs
    mach2 = machine(model2, X, y)
    fit!(mach2, rows=train3, verbosity=0)
    y23 = predict(mach2, rows=test3) |> probs
    y2_oos = vcat(y21, y22, y23)
    mach2_full = machine(model2, X, y)
    fit!(mach2_full, verbosity=0)
    y2 = predict(mach2_full, X) |> probs

    # metalearner (`judge`):
    X_oos = MLJBase.table(hcat(y1_oos, y2_oos))
    judge = KNNClassifier(K=3)
    m_judge = machine(judge, X_oos, y)
    fit!(m_judge, verbosity=0)
    X_judge = MLJBase.table(hcat(y1, y2))
    yhat_matrix = predict(m_judge, X_judge) |> probs

    # alternatively, use stack:
    stack = Stack(metalearner=judge,
                model1=model1,
                model2=model2,
                resampling=CV(nfolds=3))
    mach = machine(stack, X, y)
    fit!(mach, verbosity=0)
    yhat_matrix_stack = predict(mach, X) |> probs

    # compare:
    @test yhat_matrix_stack ≈ yhat_matrix

end

@testset "Test store_for_evaluation" begin
    X, y = make_blobs(;rng=rng, shuffle=false)
    Xs, ys = source(X), source(y)
    mach = machine(KNNClassifier(), Xs, ys)
    fit!(mach, verbosity=0)
    measures = [accuracy, log_loss]
    mach_, Xtest, ytest = MLJBase.store_for_evaluation(mach, Xs, ys, measures)()
    @test Xtest == X
    @test ytest == y
    @test mach_ == mach
    # check fallback
    @test MLJBase.store_for_evaluation(mach, Xs, ys, nothing) === nothing
end

@testset "Test internal_stack_report" begin
    n = 500
    X, y = make_regression(n, 5; rng=rng)
    resampling = CV(;nfolds=2)
    measures = [rms, l2]
    constant = ConstantRegressor()
    ridge = FooBarRegressor()
    mystack = Stack(;metalearner=ridge,
                    resampling=resampling,
                    measures=measures,
                    constant=constant,
                    ridge=ridge)

    std_evaluation = (
        constant=evaluate(constant, X, y, resampling=resampling, measures=measures, verbosity=0),
        ridge=evaluate(ridge, X, y, resampling=resampling, measures=measures, verbosity=0)
    )
    ttp = MLJBase.train_test_pairs(resampling, 1:nrows(y), X, y)
    # Testing internal_stack_report default with nothing
    @test MLJBase.internal_stack_report(mystack, 0, ttp, nothing, nothing) == NamedTuple{}()

    # Simulate the evaluation nodes which consist of
    # - The fold machine
    # - Xtest
    # - ytest
    evaluation_nodes = []
    for (train, test) in MLJBase.train_test_pairs(resampling, 1:n, y)
        for model in getfield(mystack, :models)
            mach = machine(model, X, y)
            fit!(mach, verbosity=0, rows=train)
            Xtest = selectrows(X, test)
            ytest = selectrows(y, test)
            push!(evaluation_nodes, source((mach, Xtest, ytest)))
        end
    end

    internalreport = MLJBase.internal_stack_report(
        mystack,
        0,
        ttp,
        evaluation_nodes...
    ).report.cv_report()

    test_internal_evaluation(internalreport, std_evaluation, (:constant, :ridge))

    test_internal_evaluation(internalreport, std_evaluation, (:constant, :ridge))
    @test std_evaluation.constant.fitted_params_per_fold ==
        internalreport.constant.fitted_params_per_fold
    @test std_evaluation.ridge.fitted_params_per_fold ==
        internalreport.ridge.fitted_params_per_fold

end

@testset "Test internal evaluation of the stack in regression mode" begin
    X, y = make_regression(500, 5; rng=rng)
    resampling = CV(;nfolds=3)
    measures = [rms, l2]
    constant = ConstantRegressor()
    ridge = FooBarRegressor()
    mystack = Stack(;metalearner=FooBarRegressor(),
                    resampling=resampling,
                    measure=measures,
                    ridge=ridge,
                    constant=constant)

    mach = machine(mystack, X, y)
    fit!(mach, verbosity=0)
    internalreport = report(mach).cv_report
    # evaluate decisiontree and ridge out of stack and check results match
    std_evaluation = (
        constant = evaluate(constant, X, y,
                            measure=measures,
                            resampling=resampling,
                            verbosity=0),
        ridge = evaluate(ridge, X, y, measure=measures, resampling=resampling, verbosity=0)
        )

    test_internal_evaluation(internalreport, std_evaluation, (:constant, :ridge))
    @test std_evaluation.constant.fitted_params_per_fold ==
        internalreport.constant.fitted_params_per_fold
    @test std_evaluation.ridge.fitted_params_per_fold ==
        internalreport.ridge.fitted_params_per_fold

end

@testset "Test internal evaluation of the stack in classification mode" begin
    X, y = make_blobs(;rng=rng, shuffle=false)
    resampling = StratifiedCV(;nfolds=3)
    measures = [accuracy, log_loss]
    constant = ConstantClassifier()
    knn = KNNClassifier()
    mystack = Stack(;metalearner=DecisionTreeClassifier(),
                    resampling=resampling,
                    constant=constant,
                    knn=knn,
                    measures=measures)

    mach = machine(mystack, X, y)
    fit!(mach, verbosity=0)
    internalreport = report(mach).cv_report
    # evaluate decisiontree and ridge out of stack and check results match
    std_evaluation = (
        constant = evaluate(constant, X, y,
                            measure=measures,
                            resampling=resampling,
                            verbosity=0),
        knn = evaluate(knn, X, y, measure=measures, resampling=resampling, verbosity=0)
        )

    test_internal_evaluation(internalreport, std_evaluation, (:knn, :constant))
    # Test fitted_params
    for i in 1:mystack.resampling.nfolds
        std_constant_fp = std_evaluation.constant.fitted_params_per_fold[i]
        intern_constant_fp = internalreport.constant.fitted_params_per_fold[i]
        @test std_constant_fp.target_distribution ≈ intern_constant_fp.target_distribution

        std_knn_fp = std_evaluation.knn.fitted_params_per_fold[i]
        intern_knn_fp = internalreport.knn.fitted_params_per_fold[i]
        @test std_knn_fp.tree.data == intern_knn_fp.tree.data
    end

end

@testset "Test Holdout CV" begin
    X, y = make_regression(100, 3; rng=rng)
    resampling = Holdout()
    constant = ConstantRegressor()
    ridge = FooBarRegressor()
    mystack = Stack(;metalearner=FooBarRegressor(),
                    resampling=resampling,
                    measures=[rmse],
                    ridge=ridge,
                    constant=constant)

    mach = machine(mystack, X, y)
    fit!(mach, verbosity=0)
    for modelname in (:ridge, :constant)
        model_perf = getproperty(report(mach).cv_report, modelname)
        @test length(model_perf.per_fold) == 1
        @test length(model_perf.train_test_rows) == 1
    end
end

@testset "Test cache is forwarded to submodels" begin
    X, y = make_regression(100, 3; rng=rng)
    constant = ConstantRegressor()
    ridge = FooBarRegressor()
    mystack = Stack(;metalearner=FooBarRegressor(),
                    cache=false,
                    ridge=ridge,
                    constant=constant)
    mach = machine(mystack, X, y)
    fit!(mach, verbosity = 0)
    # The data and resampled_data have not been populated
    for mach in fitted_params(mach).machines
        @test !isdefined(mach, :data)
        @test !isdefined(mach, :resampled_data)
    end
end

# a regression `Stack` which has `model` as one of the base models:
function _stack(model, resource)
    models = (constant=DeterministicConstantRegressor(),
                ridge_lambda=FooBarRegressor(;lambda=0.1),
                model=model)
    Stack(;
        metalearner=FooBarRegressor(;lambda=0.05),
        resampling=CV(;nfolds=3),
        acceleration=resource,
        models...
    )
end

# return a nested stack in which `model` appears at two levels, with
# both layers accelerated using `resource`:
_double_stack(model, resource) =
    _stack(_stack(model, resource), resource)

@testset "Test multithreaded version" begin
    X, y = make_regression(100, 5; rng=StableRNG(1234))

    stack = _double_stack(FooBarRegressor(;lambda=0.07), CPU1())

    mach = machine(stack, X, y)
    fit!(mach, verbosity=0)
    cpu_fp = fitted_params(mach)
    cpu_ypred = predict(mach)

    stack = _double_stack(FooBarRegressor(;lambda=0.07), CPUThreads())

    mach = machine(stack, X, y)
    fit!(mach, verbosity=0)
    thread_fp = fitted_params(mach)
    thread_ypred = predict(mach)

    @test cpu_ypred ≈ thread_ypred
    @test cpu_fp.metalearner ≈ thread_fp.metalearner
    @test cpu_fp.ridge_lambda ≈ thread_fp.ridge_lambda
end

end
true


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct DeploymentOperationsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets a deployments operation.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: operationId::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DeploymentOperation
"""
function deploymentOperationsGet(_api::DeploymentOperationsApi, resourceGroupName::String, deploymentName::String, operationId::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentOperationsGet", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentOperationsGet", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentOperationsGet", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentOperationsGet", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", DeploymentOperation, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/deployments/{deploymentName}/operations/{operationId}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "operationId", operationId)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all deployments operations for a deployment.
Param: resourceGroupName::String (required)
Param: deploymentName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: top::Int32
Return: DeploymentOperationsListResult
"""
function deploymentOperationsList(_api::DeploymentOperationsApi, resourceGroupName::String, deploymentName::String, api_version::String, subscriptionId::String; top=nothing, _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "deploymentOperationsList", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "deploymentOperationsList", :minLength, resourceGroupName, 1)

    Swagger.validate_param("deploymentName", "deploymentOperationsList", :maxLength, deploymentName, 64)
    Swagger.validate_param("deploymentName", "deploymentOperationsList", :minLength, deploymentName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", DeploymentOperationsListResult, "/subscriptions/{subscriptionId}/resourcegroups/{resourceGroupName}/deployments/{deploymentName}/operations", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "deploymentName", deploymentName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export deploymentOperationsGet, deploymentOperationsList


text:
## Hermite
@register0 Hermite AbstractCCOP0
export Hermite

"""
    Hermite

The [Hermite](https://en.wikipedia.org/wiki/Hermite_polynomials) polynomials have two versions the physicists (`Hermite`  or  `H`) and the probablalists (`ChebyshevHermite` or  `Hₑ`). They are  related through  `Hᵢ(x) =  2^(i/2) Hₑᵢ(√2 x)`.

The Hermite   polynomials have weight  function `w(x)=exp(-x^2/2)` and domain the real line.

## Examples
```jldoctest
julia> using Polynomials,  SpecialPolynomials

julia> x = variable(Polynomial{Rational{Int}})
Polynomials.Polynomial(x)

julia> [basis(Hermite, i)(x) for i in 0:5]
6-element Vector{Polynomial{Float64, :x}}:
 Polynomials.Polynomial(1.0)
 Polynomials.Polynomial(2.0*x)
 Polynomials.Polynomial(-2.0 + 4.0*x^2)
 Polynomials.Polynomial(-12.0*x + 8.0*x^3)
 Polynomials.Polynomial(12.0 - 48.0*x^2 + 16.0*x^4)
 Polynomials.Polynomial(120.0*x - 160.0*x^3 + 32.0*x^5)

julia> [basis(ChebyshevHermite, i)(x) for i in 0:5]
6-element Vector{Polynomial{Float64, :x}}:
 Polynomials.Polynomial(1.0)
 Polynomials.Polynomial(1.0*x)
 Polynomials.Polynomial(-1.0 + 1.0*x^2)
 Polynomials.Polynomial(-3.0*x + 1.0*x^3)
 Polynomials.Polynomial(3.0 - 6.0*x^2 + 1.0*x^4)
 Polynomials.Polynomial(15.0*x - 10.0*x^3 + 1.0*x^5)
```

!!! note
    The Hermite family needs help, as the computed values for `Bn`,and,`Cn` are  both 0.
"""
Hermite

basis_symbol(::Type{<:Hermite}) = "H"
Polynomials.domain(::Type{<:Hermite}) = Polynomials.Interval(-Inf, Inf)

abcde(::Type{<:Hermite}) = NamedTuple{(:a, :b, :c, :d, :e)}((1, 0, 0, -2, 0))

k0(P::Type{<:Hermite}) = one(eltype(P))
function k1k0(P::Type{<:Hermite}, k::Int)
    val = 2 * one(eltype(P))
    val
end

norm2(::Type{<:Hermite}, n) = sqrt(pi) * 2^n * gamma(n + 1)
weight_function(::Type{<:Hermite}) = x -> exp(-x^2)
generating_function(::Type{<:Hermite}) = (t, x) -> exp(2 * t * x - t^2)
function classical_hypergeometric(::Type{<:Hermite}, n, x)
    as = iseven(n) ? (-n ÷ 2, -(n - 1) / 2) : (-n / 2, -(n - 1) ÷ 2)
    bs = ()
    (2x)^n * pFq(as, bs, -inv(x)^2)
end

gauss_nodes_weights(p::Type{P}, n) where {P<:Hermite} = FastGaussQuadrature.gausshermite(n)

## Overrides
# Use override here, as we get  An,Bn,Cn =  1, 0,0, otherwise
An(P::Type{<:Hermite}, n::Int) = 2 * one(eltype(P))
Bn(P::Type{<:Hermite}, n::Int) = zero(eltype(P))
Cn(P::Type{<:Hermite}, n::Int) = 2 * n * one(eltype(P))

β̃n(P::Type{<:Hermite}, n::Int) = zero(eltype(P))
γ̃n(P::Type{<:Hermite}, n::Int) = zero(eltype(P))
b̂̃n(P::Type{<:Hermite}, ::Val{N}) where {N} = zero(eltype(P)) #where {M} = error("Don't call me")#zero(S)
ĉ̃n(P::Type{<:Hermite}, ::Val{N}) where {N} = zero(eltype(P)) #where {M} = error("Don't call me")#zero(S)

## https://arxiv.org/pdf/1901.01648.pdf. Connection formula (14)
##  x^n  = n! sum(H_{n-2j}/ (2^j(n-2j)!j!) j = 0:floor(n/2))

Base.convert(P::Type{<:Hermite}, q::Polynomial) = connection(P, q)
function Base.iterate(
    o::Connection{P,Q},
    state=nothing,
) where {P<:Hermite,Q<:Polynomials.StandardBasisPolynomial}
    k, n = o.k, o.n

    if state == nothing
        i = k
        j = 0
        i > n && return nothing
        val = __hermite_lambda(i, k)
    elseif state[1] + 2 > n # terminate
        return nothing
    else
        j, val1 = state
        #val1 *= (i+1)*(i+2)/4/(j+1/4)
        j += 1
        i = k + 2j
        val = __hermite_lambda(i, k)
    end

    return (i, val), (j, val)
end

function __hermite_lambda(n, k)
    val = gamma(1 + n) / 2^n
    val /= gamma(1 + k)
    val /= gamma(1 + (n - k) / 2)
    val
end

# 2x more performant
function Polynomials.integrate(p::P, C::Number=0) where {T,X,P<:Hermite{T,X}}
    # int H_n = 1/(n+1) H_{n+1}
    R = eltype(one(T) / 1)
    d = degree(p)
    qs = zeros(R, d + 2)

    for i in 0:d
        qs[i + 2] = p[i] / (2(i + 1))
    end

    q = ⟒(P){R,X}(qs)
    q = q - q(0) + R(C)

    return q
end

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

@register0 MonicHermite AbstractCCOP0
export MonicHermite
ϟ(::Type{<:MonicHermite}) = Hermite
ϟ(::Type{<:MonicHermite{T}}) where {T} = Hermite{T}
@register_monic(MonicHermite)

pqr_start(::Type{<:MonicHermite}) = 0
pqr_symmetry(::Type{<:MonicHermite}) = true
pqr_weight(P::Type{<:MonicHermite}, n, x, dπx) = (one(P) * 2) * exp(-x^2) / (dπx * dπx)

##
## --------------------------------------------------
##
@register0 ChebyshevHermite AbstractCCOP0
export ChebyshevHermite
"""
    ChebyshevHermite

Type for the Probabalist's  [Hermite](https://en.wikipedia.org/wiki/Hermite_polynomials) polynomials.

"""
ChebyshevHermite

basis_symbol(::Type{<:ChebyshevHermite}) = "Hₑ"
Polynomials.domain(::Type{<:ChebyshevHermite}) = Polynomials.Interval(-Inf, Inf)

# https://arxiv.org/pdf/1901.01648.pdf eqn 17
abcde(::Type{<:ChebyshevHermite}) = NamedTuple{(:a, :b, :c, :d, :e)}((0, 0, 1, -1, 0))

kn(P::Type{<:ChebyshevHermite}, n::Int)    = one(eltype(P))
k1k0(P::Type{<:ChebyshevHermite}, k::Int)  = one(eltype(P))
k1k_1(P::Type{<:ChebyshevHermite}, k::Int) = one(eltype(P))
ismonic(::Type{<:ChebyshevHermite})        = true

weight_function(P::Type{ChebyshevHermite}) = x -> exp(-(one(eltype(P)) * x^2) / 2)
generating_function(::Type{<:ChebyshevHermite}) = (t, x) -> exp(t * x - t^2 / 2)
function classical_hypergeometric(::Type{<:ChebyshevHermite}, n, x)
    2^(-n / 2) * classical_hypergeometric(Hermite, n, x / sqrt(2))
end

function gauss_nodes_weights(P::Type{<:ChebyshevHermite}, n) where {α,β}
    xs, ws = glaser_liu_rokhlin_gauss_nodes(basis(ChebyshevHermite, n))
    xs, ws
end

pqr_start(::Type{<:ChebyshevHermite}) = 0
pqr_symmetry(::Type{<:ChebyshevHermite}) = true
pqr_weight(P::Type{<:ChebyshevHermite}, n, x, dπx) =
    sqrt(2) * gamma(1 + n) * sqrt(pi) / (dπx)^2

## Overrides
#An(P::Type{<:ChebyshevHermite}, n::Int) = one(eltype(P))
#Bn(P::Type{<:ChebyshevHermite}, n::Int) = zero(eltype(P))
#Cn(P::Type{<:ChebyshevHermite}, n::Int) = -one(eltype(P))*n

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

##  Multiplication
## TODO: work  on  types
AbstractHermite{T,X} = Union{Hermite{T,X},ChebyshevHermite{T,X}}

# Default is slow. Instead
# directly use a linearization formula from
# https://arxiv.org/pdf/1901.01648.pdf
# and  for Hermite, Hn(x) =2^(n/2)Hₑ_n(sqrt(2) x)
# so Hm⋅Hn = ∑ C_{m,n,j} 2^j H_{m+n-2j}
function ⊗(
    ::Type{<:AbstractHermite},
    p::P,
    q::Q,
) where {T,X,S,Y,P<:AbstractHermite{T,X},Q<:AbstractHermite{S,Y}}
    R = eltype(one(promote_type(T, S)) / 1)
    N, M = degree(p), degree(q)
    N == -1 && return zero(⟒(P){R,Y})
    M == -1 && return zero(⟒(P){R,X})
    N == 0 && return q * p[0]
    M == 0 && return p * q[0]
    X != Y && throw(ArgumentError("Variable names must match"))

    as = zeros(R, 1 + N + M)
    for i in eachindex(p)
        pᵢ = p[i]
        iszero(pᵢ) && continue
        for j in eachindex(q)
            qⱼ = q[j]
            pᵢqⱼ = pᵢ * qⱼ
            iszero(qⱼ) && continue
            λᵤ = one(R)
            for k in 0:min(i, j)
                as[1 + i + j - 2k] += pᵢqⱼ * λᵤ
                λᵤ *= (i - k) * (j - k) / (k + 1) * hermite_α(P)
            end
        end
    end
    ⟒(P){R,X}(as)
end

hermite_α(P::Type{<:ChebyshevHermite}) = 1
hermite_α(P::Type{<:Hermite}) = 2

#⊗(p::Hermite, q::Hermite) = linearization_product(p, q)
#⊗(p::ChebyshevHermite, q::ChebyshevHermite) = linearization_product(p, q)

# function Base.iterate(o::Linearization{P,R}, state =  nothing) where {P <: AbstractHermite,R}

#     l, m, n = o.l, o.m, o.n
#     l  > m + n && return nothing

#     if state == nothing

#         #k =  0
#         j = 0
#         p = min(l, n)
#         q = l - p
#         val = linearization_α(P, R, j, l, p, q)

#     else
#         # we work with l + 2j as there is a parity consideration
#         j, p, q, val  = state
#         s = l + j # l +2j + p + q = l +2j + l = 2l + 2j, so s=l+j
#         if p == 0 ||  q == m || q  > s
#             j += 1
#             l + 2j > n+m && return nothing #  l + 2j  is too  big now
#             p = min(l + j, n)  # p <= s
#             q = l + 2j - p
#             q > s + 1 && return nothing
#             val = linearization_α(P, R, j, l, p, q)
#         else
#             p -= 1
#             q += 1

#             λ = q/(p+1)*(s-(q-1))/(s-p)
#             val *= λ
#         end

#     end

#     return (p,q,val), (j, p, q, val)
# end

# #  https://arxiv.org/pdf/1901.01648.pdf equation (74)
# function linearization_α(P::Type{<:AbstractHermite}, R, j, l, n, m)
#     s = l + j # pass in j to avoid s = divrem(l + m + n,  2) call
#     val = one(R)
#     val *= gamma(1+m)*gamma(1+n)/gamma(1+s-m)/gamma(1+s-n)/gamma(1+s-l)
#     val *= linearization_λ(P, l, m,n)
# end
# linearization_λ(::Type{<:Hermite}, l, m, n) = 2.0^((m+n-l)/2)
# linearization_λ(::Type{<:ChebyshevHermite}, l, m, n) = 1


text:
include("Map/Map.jl")
include("Map/ResidueRingPoly.jl")
include("Map/NfOrd.jl")
include("Map/NumberField.jl")
include("Map/GrpAb.jl")
include("Map/NfRel.jl")
include("Map/FqPolyRing.jl")