text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# tests for accurate updating of method tables

using Base: get_world_counter
tls_world_age() = ccall(:jl_get_tls_world_age, UInt, ())
@test typemax(UInt) > get_world_counter() == tls_world_age() > 0

# test simple method replacement
begin
    g265a() = f265a(0)
    f265a(x::Any) = 1
    @test g265a() == 1
    @test Base.return_types(g265a, ()) == Any[Int]
    @test Core.Compiler.return_type(g265a, ()) == Int

    f265a(x::Any) = 2.0
    @test g265a() == 2.0

    @test Base.return_types(g265a, ()) == Any[Float64]
    @test Core.Compiler.return_type(g265a, ()) == Float64
end

# test signature widening
begin
    f265b(x::Int) = 1
    let ty = Any[1, 2.0e0]
        global g265b(i::Int) = f265b(ty[i])
    end
    @test g265b(1) == 1
    @test Base.return_types(g265b, (Int,)) == Any[Int]
    @test Core.Compiler.return_type(g265b, (Int,)) == Int

    f265b(x::Any) = 2.0
    @test g265b(1) == 1
    @test g265b(2) == 2.0
    @test Base.return_types(g265b, (Int,)) == Any[Union{Int, Float64}]
    @test Core.Compiler.return_type(g265b, (Int,)) == Union{Int, Float64}
end

# test signature narrowing
begin
    g265c() = f265c(0)
    f265c(x::Any) = 1
    @test g265c() == 1
    @test Base.return_types(g265c, ()) == Any[Int]
    @test Core.Compiler.return_type(g265c, ()) == Int

    f265c(x::Int) = 2.0
    @test g265c() == 2.0

    @test Base.return_types(g265c, ()) == Any[Float64]
    @test Core.Compiler.return_type(g265c, ()) == Float64
end

# test constructor narrowing
mutable struct A265{T}
    field1::T
end
A265_() = A265(1)
@test (A265_()::A265{Int}).field1 === 1
A265(fld::Int) = A265(Float64(fld))
@test (A265_()::A265{Float64}).field1 === 1.0e0

# test constructor widening
mutable struct B265{T}
    field1::T
    # dummy arg is present to prevent (::Type{T}){T}(arg) from matching the test calls
    B265{T}(field1::Any, dummy::Nothing) where {T} = new(field1) # prevent generation of outer ctor
end
  # define some constructors
B265(x::Int, dummy::Nothing) = B265{Int}(x, dummy)
let ty = Any[1, 2.0e0, 3.0f0]
    global B265_(i::Int) = B265(ty[i], nothing)
end
  # test for correct answers
@test (B265_(1)::B265{Int}).field1 === 1
@test_throws MethodError B265_(2)
@test_throws MethodError B265_(3)
@test Base.return_types(B265_, (Int,)) == Any[B265{Int}]
@test Core.Compiler.return_type(B265_, (Int,)) == B265{Int}

  # add new constructors
B265(x::Float64, dummy::Nothing) = B265{Float64}(x, dummy)
B265(x::Any, dummy::Nothing) = B265{UInt8}(x, dummy)

  # make sure answers are updated
@test (B265_(1)::B265{Int}).field1 === 1
@test (B265_(2)::B265{Float64}).field1 === 2.0e0
@test (B265_(3)::B265{UInt8}).field1 === 0x03

@test Base.return_types(B265_, (Int,)) == Any[B265]
@test Core.Compiler.return_type(B265_, (Int,)) == B265

# test oldworld call / inference
function wfunc(c1,c2)
    while true
        (f, args) = take!(c1)
        put!(c2, f(args...))
    end
end
function put_n_take!(v...)
    put!(chnls[1], v)
    take!(chnls[2])
end

g265() = [f265(x) for x in 1:3.]
wc265 = get_world_counter()
f265(::Any) = 1.0
@test wc265 + 1 == get_world_counter()
chnls, tasks = Base.channeled_tasks(2, wfunc)
t265 = tasks[1]

wc265 = get_world_counter()
@test put_n_take!(get_world_counter, ()) == wc265
@test put_n_take!(tls_world_age, ()) == wc265
f265(::Int) = 1
@test put_n_take!(get_world_counter, ()) == wc265 + 1 == get_world_counter() == tls_world_age()
@test put_n_take!(tls_world_age, ()) == wc265

@test g265() == Int[1, 1, 1]
@test Core.Compiler.return_type(f265, (Any,)) == Union{Float64, Int}
@test Core.Compiler.return_type(f265, (Int,)) == Int
@test Core.Compiler.return_type(f265, (Float64,)) == Float64

@test put_n_take!(g265, ()) == Float64[1.0, 1.0, 1.0]
@test put_n_take!(Core.Compiler.return_type, (f265, (Any,))) == Float64
@test put_n_take!(Core.Compiler.return_type, (f265, (Int,))) == Float64
@test put_n_take!(Core.Compiler.return_type, (f265, (Float64,))) == Float64
@test put_n_take!(Core.Compiler.return_type, (f265, (Float64,))) == Float64

# test that reflection ignores worlds
@test Base.return_types(f265, (Any,)) == Any[Int, Float64]
@test put_n_take!(Base.return_types, (f265, (Any,))) == Any[Int, Float64]

# test for method errors
h265() = true
loc_h265 = "$(@__FILE__):$(@__LINE__() - 1)"
@test h265()
@test_throws TaskFailedException(t265) put_n_take!(h265, ())
@test_throws TaskFailedException(t265) fetch(t265)
@test istaskdone(t265)
let ex = t265.exception
    @test ex isa MethodError
    @test ex.f == h265
    @test ex.args == ()
    @test ex.world == wc265
    str = sprint(showerror, ex)
    wc = get_world_counter()
    cmps = """
        MethodError: no method matching h265()
        The applicable method may be too new: running in world age $wc265, while current world is $wc."""
    @test startswith(str, cmps)
    cmps = "\n  h265() at $loc_h265 (method too new to be called from this world context.)"
    @test occursin(cmps, str)
end

# test for generated function correctness
# and min/max world computation validity of cache_method
f_gen265(x) = 1
@generated g_gen265(x) = f_gen265(x)
@generated h_gen265(x) = :(f_gen265(x))
f_gen265(x::Int) = 2
f_gen265(x::Type{Int}) = 3
@generated g_gen265b(x) = f_gen265(x)
@test h_gen265(0) == 2
@test g_gen265(0) == 1
@test f_gen265(Int) == 3
@test g_gen265b(0) == 3

# Test that old, invalidated specializations don't get revived for
# intermediate worlds by later additions to the method table that
# would have capped those specializations if they were still valid
f26506(@nospecialize(x)) = 1
g26506(x) = f26506(x[1])
z = Any["ABC"]
f26506(x::Int) = 2
g26506(z) # Places an entry for f26506(::String) in mt.name.cache
f26506(x::String) = 3
let cache = typeof(f26506).name.mt.cache
    # The entry we created above should have been truncated
    @test cache.min_world == cache.max_world
end
c26506_1, c26506_2 = Condition(), Condition()
# Captures the world age
result26506 = Any[]
t = Task(()->begin
    wait(c26506_1)
    push!(result26506, g26506(z))
    notify(c26506_2)
end)
yield(t)
f26506(x::Float64) = 4
let cache = typeof(f26506).name.mt.cache
    # The entry we created above should have been truncated
    @test cache.min_world == cache.max_world
end
notify(c26506_1)
wait(c26506_2)
@test result26506[1] == 3

## Invalidation tests

function instance(f, types)
    m = which(f, types)
    inst = nothing
    tt = Tuple{typeof(f), types...}
    specs = m.specializations
    if isa(specs, Nothing)
    elseif isa(specs, Core.SimpleVector)
        for i = 1:length(specs)
            if isassigned(specs, i)
                mi = specs[i]::Core.MethodInstance
                if mi.specTypes === tt
                    inst = mi
                    break
                end
            end
        end
    else
        Base.visit(specs) do mi
            if mi.specTypes === tt
                inst = mi
            end
        end
    end
    return inst
end

function worlds(mi::Core.MethodInstance)
    w = Tuple{UInt,UInt}[]
    if isdefined(mi, :cache)
        ci = mi.cache
        push!(w, (ci.min_world, ci.max_world))
        while isdefined(ci, :next)
            ci = ci.next
            push!(w, (ci.min_world, ci.max_world))
        end
    end
    return w
end

# avoid adding this to Base
function equal(ci1::Core.CodeInfo, ci2::Core.CodeInfo)
    return ci1.code == ci2.code &&
           ci1.codelocs == ci2.codelocs &&
           ci1.ssavaluetypes == ci2.ssavaluetypes &&
           ci1.ssaflags == ci2.ssaflags &&
           ci1.method_for_inference_limit_heuristics == ci2.method_for_inference_limit_heuristics &&
           ci1.linetable == ci2.linetable &&
           ci1.slotnames == ci2.slotnames &&
           ci1.slotflags == ci2.slotflags &&
           ci1.slottypes == ci2.slottypes &&
           ci1.rettype == ci2.rettype
end
equal(p1::Pair, p2::Pair) = p1.second == p2.second && equal(p1.first, p2.first)

## Union-splitting based on state-of-the-world: check that each invalidation corresponds to new code
applyf35855(c) = f35855(c[1])
f35855(::Int) = 1
f35855(::Float64) = 2
applyf35855([1])
applyf35855([1.0])
applyf35855(Any[1])
wint   = worlds(instance(applyf35855, (Vector{Int},)))
wfloat = worlds(instance(applyf35855, (Vector{Float64},)))
wany2  = worlds(instance(applyf35855, (Vector{Any},)))
src2 = code_typed(applyf35855, (Vector{Any},))[1]
f35855(::String) = 3
applyf35855(Any[1])
@test worlds(instance(applyf35855, (Vector{Int},))) == wint
@test worlds(instance(applyf35855, (Vector{Float64},))) == wfloat
wany3 = worlds(instance(applyf35855, (Vector{Any},)))
src3 = code_typed(applyf35855, (Vector{Any},))[1]
@test !(wany3 == wany2) || equal(src3, src2) # code doesn't change unless you invalidate
f35855(::AbstractVector) = 4
applyf35855(Any[1])
wany4 = worlds(instance(applyf35855, (Vector{Any},)))
src4 = code_typed(applyf35855, (Vector{Any},))[1]
@test !(wany4 == wany3) || equal(src4, src3) # code doesn't change unless you invalidate
f35855(::Dict) = 5
applyf35855(Any[1])
wany5 = worlds(instance(applyf35855, (Vector{Any},)))
src5 = code_typed(applyf35855, (Vector{Any},))[1]
@test (wany5 == wany4) == equal(src5, src4)
f35855(::Set) = 6    # with current settings, this shouldn't invalidate
applyf35855(Any[1])
wany6 = worlds(instance(applyf35855, (Vector{Any},)))
src6 = code_typed(applyf35855, (Vector{Any},))[1]
@test wany6 == wany5
@test equal(src6, src5)

applyf35855_2(c) = f35855_2(c[1])
f35855_2(::Int) = 1
f35855_2(::Float64) = 2
applyf35855_2(Any[1])
wany3 = worlds(instance(applyf35855_2, (Vector{Any},)))
src3 = code_typed(applyf35855_2, (Vector{Any},))[1]
f35855_2(::AbstractVector) = 4
applyf35855_2(Any[1])
wany4 = worlds(instance(applyf35855_2, (Vector{Any},)))
src4 = code_typed(applyf35855_2, (Vector{Any},))[1]
@test !(wany4 == wany3) || equal(src4, src3) # code doesn't change unless you invalidate

## ambiguities do not trigger invalidation
using Printf
Printf.gen("%f")
mi = instance(+, (AbstractChar, UInt8))
w = worlds(mi)

abstract type FixedPoint35855{T <: Integer} <: Real end
struct Normed35855 <: FixedPoint35855{UInt8}
    i::UInt8
    Normed35855(i::Integer, _) = new(i % UInt8)
end
(::Type{X})(x::Real) where {T, X<:FixedPoint35855{T}} = X(round(T, typemax(T)*x), 0)
@test worlds(mi) == w

mi = instance(convert, (Type{Nothing}, String))
w = worlds(mi)
abstract type Colorant35855 end
Base.convert(::Type{C}, c) where {C<:Colorant35855} = false
@test worlds(mi) == w

# invoke_in_world
f_inworld(x) = "world one; x=$x"
g_inworld(x; y) = "world one; x=$x, y=$y"
wc_aiw1 = get_world_counter()
# redefine f_inworld, g_inworld, and check that we can invoke both versions
f_inworld(x) = "world two; x=$x"
g_inworld(x; y) = "world two; x=$x, y=$y"
wc_aiw2 = get_world_counter()
@test Base.invoke_in_world(wc_aiw1, f_inworld, 2) == "world one; x=2"
@test Base.invoke_in_world(wc_aiw2, f_inworld, 2) == "world two; x=2"
@test Base.invoke_in_world(wc_aiw1, g_inworld, 2, y=3) == "world one; x=2, y=3"
@test Base.invoke_in_world(wc_aiw2, g_inworld, 2, y=3) == "world two; x=2, y=3"


text:
using Flux
using Statistics

# Note: original Flux.logitbinarycrossentropy broadcasting calls fail on GPU, this is a batch version
batchlogitbinarycrossentropy(logY, T) = mean((1 .- T).*logY .- logσ.(logY))

# Note: original Flux.mse broadcasting calls fail on GPU, this is a batch version
function batchmse(Y, T)
    diff = (Y .- T)
    mean(diff .* diff)
end

ext: jl

max_stars_repo_name: dhonza/FluxGoodies

max_issues_repo_name: dhonza/FluxGoodies

max_forks_repo_name: dhonza/FluxGoodies


text:
@testset "FMM" begin
    λ0 = 1
    k0 = 2π/λ0
    kin = 3k0
    θ_i = π/4 #incident wave e^{i k_0 (1/sqrt{2},1/sqrt{2}) \cdot \mathbf{r}}
    N_squircle = 200
    N_star = 210
    P = 10

    M = 10
    shapes = [rounded_star(0.1λ0, 0.03λ0, 5, N_star);
                squircle(0.15λ0, N_squircle)]
    ids = rand(1:2, M^2)
    φs = 2π*rand(M^2) #random rotation angles
    dist = 2*maximum(shapes[i].R for i=1:2)
    try
        centers = randpoints(M^2, dist, 5λ0, 5λ0, [0.0 0.0; 0.01+dist 0.01],
                    failures = 1_000)
        sp = ScatteringProblem(shapes, ids, centers, φs)
        @test verify_min_distance(sp, [0.0 0.0; 0.01+dist 0.01])
    catch
        @warn("Could not find random points (1)")
    end
    try
        centers = randpoints(M^2, dist, 5λ0, 5λ0, failures = 1_000)
        sp = ScatteringProblem(shapes, ids, centers, φs)
        @test verify_min_distance(sp)
    catch
        @warn("Could not find random points (2)")
    end

    centers =  rect_grid(M, M, 0.8λ0, λ0)
    sp = ScatteringProblem(shapes, ids, centers, φs)
    fmm_options = FMMoptions(true, acc = 6, nx = div(M,2))

    groups,boxSize = divideSpace(centers, fmm_options; drawGroups = false)
    P2,Q = FMMtruncation(6, boxSize, k0)
    mFMM1 = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize; tri = false)
    mFMM2 = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize; tri = true)
    @test norm(mFMM1.Znear - mFMM2.Znear, Inf) == 0

    points = [range(-λ0*M, stop=λ0*M, length=200) zeros(200)]
    u1 = calc_near_field(k0, kin, P, sp, points, PlaneWave(θ_i); opt = fmm_options,
            verbose = false)
    u2 = calc_near_field(k0, kin, P, sp, points, PlaneWave(θ_i); opt = fmm_options,
            verbose = true, method = "density")
    @test norm(u1 - u2)/norm(u1) < 1e-6

    u3 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0);
            opt = fmm_options, verbose = false)
    u4 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0);
            opt = fmm_options, verbose = true, method = "recurrence")
    @test norm(u3 - u4)/norm(u3) < 1e-6

    #test transpose fmm
    mFMM, sMatrices, sLU, buf = ParticleScattering.prepare_fmm_reusal_φs(k0,
                                    kin, P, shapes, centers, ids, fmm_options)
    v1 = 5*rand(ComplexF64, M^2*(2P+1))
    v2 = rand(ComplexF64, M^2*(2P+1))
    o1 = Array{ComplexF64}(undef, M^2*(2P+1))
    o2 = Array{ComplexF64}(undef, M^2*(2P+1))
    ParticleScattering.FMM_mainMVP_pre!(o1, v1, sMatrices, φs, ids, P, mFMM,
        buf.pre_agg, buf.trans)
    res1 = transpose(v2)*o1
    ParticleScattering.FMM_mainMVP_transpose!(o2, v2, sMatrices, φs, ids, P,
        mFMM, buf.pre_agg, buf.trans)
    res2 = transpose(v1)*o2
    @test res1 ≈ res2
end


text:
println("Dp) Test of the ImpactExcitationAutoion module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
#
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(26.),
                        initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                        intermediateConfigs=[Configuration("1s 2s^2 2p"), Configuration("1s 2s 2p^2") ],
                        finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p") ],
                        process = ImpactExcAuto(),
                        processSettings=ImpactExcitationAutoion.Settings()  )

wb = perform(wa)


text:
using StatsBase: mean

mutable struct HKAgent <: AbstractAgent
    id::Int
    old_opinion::Float64
    new_opinion::Float64
    previous_opinion::Float64
end

"""
``` julia
hk(;
    numagents = 100,
    ϵ = 0.2
)
```
Same as in [HK (Hegselmann and Krause) opinion dynamics model](@ref).
"""
function hk(; numagents = 100, ϵ = 0.2)
    model = ABM(HKAgent, scheduler = Schedulers.fastest, properties = Dict(:ϵ => ϵ))
    for i in 1:numagents
        o = rand(model.rng)
        add_agent!(model, o, o, -1)
    end
    return model, hk_agent_step!, hk_model_step!, terminate
end

function boundfilter(agent, model)
    filter(
        j -> abs(agent.old_opinion - j) < model.ϵ,
        [a.old_opinion for a in allagents(model)],
    )
end

function hk_agent_step!(agent, model)
    agent.previous_opinion = agent.old_opinion
    agent.new_opinion = mean(boundfilter(agent, model))
end

function hk_model_step!(model)
    for a in allagents(model)
        a.old_opinion = a.new_opinion
    end
end

function terminate(model, s)
    if any(
        !isapprox(a.previous_opinion, a.new_opinion; rtol = 1e-12) for a in allagents(model)
    )
        return false
    else
        return true
    end
end


text:
# @testset "Verify log_pdf, grad_x_log_pdf and hess_x_log_pdf functions" begin
#     atol = 1e-8
#     Ne = 50
#     m = 10
#     Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8)]#; CstPhyHermiteBasis(8)]#; CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
#     for Nx = 1:3
#         for b in Blist
#             X = randn(Nx, Ne)
#             B = MultiBasis(b, Nx)
#             Nψ = 5
#             idx = vcat(zeros(Int64, 1, Nx), rand(0:2,Nψ-1, Nx))
#             coeff = randn(Nψ)
#             C = HermiteMapComponent(m, Nx, idx, coeff; b = string(typeof(b)))
#
# 			logC = log_pdf(C, X)
#             dxlogC  = grad_x_log_pdf(C, X)
#             d2xlogC = hess_x_log_pdf(C, X)
#
#
#             function evaluatef0(x)
#                 y = copy(x)
#                 y[end] = 0.0
#                 return C.I.f(y)
#             end
#
#             integrand(t,x) = C.I.g(ForwardDiff.gradient(y->C.I.f(y), vcat(x[1:end-1],t*x[end]))[end])
#
#             function Ct(x)
#                 lb = 0.0
#                 ub = 1.0
#                 prob = QuadratureProblem(integrand,lb,ub,x)
#                 out = evaluatef0(x) + x[end]*solve(prob,QuadGKJL())[1]
#                 return out
#             end
#
#             log_pdfCt(x) = log_pdf(Ct(x)) + log(C.I.g(ForwardDiff.gradient(z->C.I.f(z),x)[end]))
# 			@show "hello"
#             @inbounds for i=1:Ne
# 				@test norm(log_pdfCt(X[:,i]) - logC[i])<atol
#                 @test norm(ForwardDiff.gradient(log_pdfCt, X[:,i]) - dxlogC[i,:])<atol
#                 # Works but too costly, use FiniteDiff instead
# #                 @test norm(ForwardDiff.fhessian(log_pdfCt, X[:,i]) - d2xlogC[i,:,:])<atol
#                 @test norm(FiniteDiff.finite_difference_hessian(log_pdfCt, X[:,i]) - d2xlogC[i,:,:])<1000*atol
#             end
#             end
#
#         end
#     end
# end

@testset "Verify that initial map is identity" begin
	Blist = ["ProHermiteBasis"; "PhyHermiteBasis";
            "CstProHermiteBasis"; "CstPhyHermiteBasis";
            "CstLinProHermiteBasis"; "CstLinPhyHermiteBasis"]
		for b in Blist
			H = HermiteMapComponent(3, 2; α = 1e-6, b = b)
			x = randn(2,1)
			Hx = evaluate(H.I, x)
			@test abs(Hx[1] - x[2])<1e-10
		end
end

@testset "Verify getcoeff, setcoeff!, clearcoeff!" begin
	m = 5
	Nx = 10
	Nψ = 5
	idx = rand(0:m, Nψ, Nx)

	coeff = randn(Nψ)

	M = HermiteMapComponent(m, Nx, idx, deepcopy(coeff))
	@test norm(getcoeff(M) - coeff)<1e-12

	coeff2 = randn(Nψ)
	setcoeff!(M, deepcopy(coeff2))

	@test norm(getcoeff(M) - coeff2)<1e-12

	clearcoeff!(M)

	@test norm(getcoeff(M) - zeros(Nψ))<1e-12
end

# @testset "Verify loss function and its gradient" begin
#
#     Nx = 2
#     Ne = 8
#     m = 5
#
#     X     =  [0.267333   1.43021;
#               0.364979   0.607224;
#              -1.23693    0.249277;
#              -2.0526     0.915629;
#              -0.182465   0.415874;
#               0.412907   1.01672;
#               1.41332   -0.918205;
#               0.766647  -1.00445]';
#     B = MultiBasis(CstProHermiteBasis(3), Nx)
#
#     idx = [0 0; 0 1; 1 0; 2 1; 1 2]
#     Nψ = 5
#
#     coeff =   [0.6285037650645056;
#      -0.4744029092496623;
#       1.1405280011620331;
#      -0.7217760771894809;
#       0.11855056306742319]
#     f = ExpandedFunction(B, idx, coeff)
#     R = IntegratedFunction(f)
#     H = HermiteMapComponent(R; α = 0.0)
#     S = Storage(H.I.f, X);
#
#    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, H, X)), coeff, Optim.BFGS())
#    coeff_opt = Optim.minimizer(res)
#
#     @test norm(coeff_opt - [3.015753764546621;
#                           -2.929908252283099;
#                           -3.672401233483867;
#                            2.975554571687243;
#                            1.622308437415610])<1e-4
#
#     # Verify with L-2 penalty term
#
#     H = HermiteMapComponent(R; α = 0.1)
#     S = Storage(H.I.f, X);
#
#     res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, H, X)), coeff, Optim.BFGS())
#     coeff_opt = Optim.minimizer(res)
#
#     @test norm(coeff_opt - [ -0.11411931034615422;
#                              -0.21942146397348522;
#                              -0.17368042128948974;
#                               0.37348086659548607;
#                               0.02434745831060741])<1e-4
# end
#
# @testset "Verify log_pdf function" begin
#
#   Nx = 2
#   X  =  Matrix([0.267333   1.43021;
#           0.364979   0.607224;
#          -1.23693    0.249277;
#          -2.0526     0.915629;
#          -0.182465   0.415874;
#           0.412907   1.01672;
#           1.41332   -0.918205;
#           0.766647  -1.00445]');
#   B = MultiBasis(CstProHermiteBasis(6), Nx)
#
#   idx = [0 0; 0 1; 1 0; 0 2; 2 0; 1 1]
#
#   Nψ = 5
#
#   coeff = [ -0.9905841755746164;
#         0.6771992097558741;
#        -2.243695806805015;
#        -0.34591297359447354;
#        -1.420159186008486;
#        -0.5361337327704369]
#   f = ExpandedFunction(B, idx, coeff)
#   C = HermiteMapComponent(f)
#
#   @test norm(log_pdf(C, X) - [  -1.572509004118956
#                                 -2.870725221050853
#                                 -1.285696671132943
#                                 -1.088115266085997
#                                 -2.396567575345843
#                                 -2.238446803642176
#                                 -5.999143198546611
#                                 -7.082679248037675])<1e-8
# end

@testset "Verify evaluation of HermiteMapComponent" begin

	Ne = 500
	Blist = [CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]

	for Nx = 1:4
		for b in Blist
			X = randn(Nx, Ne)
			B = MultiBasis(b, Nx)
			Nψ = 5
			idx = vcat(zeros(Int64, 1, Nx), rand(0:3,Nψ-1, Nx))

			coeff = randn(Nψ)
			f = ExpandedFunction(B, idx, coeff)
			R = IntegratedFunction(f)

			C = HermiteMapComponent(R)

			# Test evaluate
			ψt = zeros(Ne)

			for i=1:Ne
				x = view(X,:,i)
				ψt[i] = R.f(vcat(x[1:end-1], 0.0)) + quadgk(t->R.g(ForwardDiff.gradient(y->R.f(y), vcat(x[1:end-1],t))[end]), 0, x[end])[1]
			end

			ψ = evaluate(C, X)

			@test norm(ψ - ψt)<1e-10
		end
	end
end

# @testset "Verify log_pdf, grad_x_log_pdf and hess_x_log_pdf functions" begin
#     atol = 1e-8
#     Ne = 50
#     m = 10
#     Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8)]#; CstPhyHermiteBasis(8)]#; CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
#     for Nx = 1:3
#         for b in Blist
#             X = randn(Nx, Ne)
#             B = MultiBasis(b, Nx)
#             Nψ = 5
#             idx = vcat(zeros(Int64, 1, Nx), rand(0:2,Nψ-1, Nx))
#             coeff = randn(Nψ)
#             C = HermiteMapComponent(m, Nx, idx, coeff; b = string(typeof(b)))
#
# 			logC = log_pdf(C, X)
#             dxlogC  = grad_x_log_pdf(C, X)
#             d2xlogC = hess_x_log_pdf(C, X)
#
#
#             function evaluatef0(x)
#                 y = copy(x)
#                 y[end] = 0.0
#                 return C.I.f(y)
#             end
#
#             integrand(t,x) = C.I.g(ForwardDiff.gradient(y->C.I.f(y), vcat(x[1:end-1],t*x[end]))[end])
#
#             function Ct(x)
#                 lb = 0.0
#                 ub = 1.0
#                 prob = QuadratureProblem(integrand,lb,ub,x)
#                 out = evaluatef0(x) + x[end]*solve(prob,QuadGKJL())[1]
#                 return out
#             end
#
#             log_pdfCt(x) = log_pdf(Ct(x)) + log(C.I.g(ForwardDiff.gradient(z->C.I.f(z),x)[end]))
#
#             @inbounds for i=1:Ne
# 				@test norm(log_pdfCt(X[:,i]) - logC[i])<atol
#                 @test norm(ForwardDiff.gradient(log_pdfCt, X[:,i]) - dxlogC[i,:])<atol
#                 # Works but too costly, use FiniteDiff instead
# #                 @test norm(ForwardDiff.fhessian(log_pdfCt, X[:,i]) - d2xlogC[i,:,:])<atol
#                 @test norm(FiniteDiff.finite_difference_hessian(log_pdfCt, X[:,i]) - d2xlogC[i,:,:])<1000*atol
#             end
#             end
#
#         end
#     end
# end

@testset "Test loss function and its gradient" begin

Nx = 2
Ne = 100

atol = 1e-8

X = Matrix([   0.096932076386797  -0.062113573952594
   1.504907056818036  -1.270761688825393
   0.122588515870727   0.138604054973047
   0.095040305032907  -0.097989669097228
   0.056360334325948  -0.098106117790437
   0.034757567319988  -0.152864122871058
   0.260831995265119  -0.564465060524812
   0.191690179163063   0.097123146831467
   0.131835643270215  -0.303532750891873
   0.121404994711150  -0.414223062306829
  -0.046027217149733   0.267918560101165
   0.222891154723706   0.054734106232808
   0.264238681579693  -0.561044970267127
   0.021345749974379   0.023020769850203
   0.395423632217738  -0.619284581169702
   0.231607199656522   0.319089487692097
   0.389505867820041   0.572610562725449
  -0.146845693861564  -0.007959126279965
   0.138489687840940   0.330481872227252
   1.617286792209003  -1.272765331253319
   0.013078993247482   0.006243503157725
   0.304266481947197  -0.508127858160603
  -0.134224572142272  -0.057784550290700
   0.058245459428341  -0.388163133932979
   0.143773871017238  -0.569911233902741
  -0.020852128497605   0.159888469685458
   0.042924979621615  -0.285726823344180
   0.566514904142701  -0.815010017214964
   0.022782788493659  -0.468191560991914
  -0.072005082647238  -0.139317736738881
  -0.088861617714439  -0.042323531263036
   0.146019714866274  -0.441835803838725
   0.062827433058116   0.255887684932700
   0.024887528661752  -0.036316627235643
   0.079288422911736   0.009625876863572
   0.124929172288805  -0.196711126223702
   0.886085483369468   1.012837297288688
  -0.062868708331948   0.065147939699519
   0.117047537973873  -0.301180703268764
   0.147117619773736   0.466610279392604
   0.690523534483231   0.722717906410571
  -0.110919157441807   0.000008657527041
   2.055034349742512   1.423365847218695
   0.091130503033388  -0.192933472923564
   0.075510323340039   0.075087088602955
  -0.181121769738736  -0.158594842114751
  -0.080866378396288  -0.081393537160189
   0.184561993178819   0.409693133491638
  -0.032018126656783   0.040851821749292
   0.702808256935611   0.978185567019995
  -0.029041940950512   0.097268795862442
   1.406479923500753  -1.166340130567743
  -0.109758856235453  -0.297180107753266
   0.236494055349260   0.548223372779986
  -0.108550945852467  -0.097740533128956
   0.908726493797006  -0.957752761353643
   0.178092997560103  -0.359297653971968
   1.103683368372646  -1.064482193617671
   1.907065203740822   1.349542162065178
   0.523722913693736   0.725695151794695
   0.261131020580618   0.576363120484564
   0.118044539009197   0.196304662470752
   0.289261334786348   0.399639383890177
   0.902906400981006  -0.957301599045371
  -0.054657884786803  -0.292760625531884
  -0.021735291600425   0.029650166664389
   0.065200888050752  -0.295894159582647
   1.486186253485923  -1.217814186300608
   0.889545420155124   0.939789761164950
  -0.174386606410644  -0.092037014472893
  -0.065037226616579   0.009771974040525
   0.074486430996400  -0.287910597788305
   0.174336742307535  -0.400464726305446
   0.096781997529499  -0.153178844921250
   0.796408810115516   0.881930856163525
   0.005874471676882   0.067750993468743
   0.156654113730103  -0.239182272065197
   0.333688106104332  -0.629954291766549
   0.086388606764696   0.305488995071947
   0.211268899950691  -0.299878322704640
   0.104223240397571   0.199354790284364
   0.336858958710283  -0.620166113933000
   0.145071750152222   0.250136305618056
   0.032242317206686  -0.233223578816564
   0.064558616046395  -0.007577632839606
  -0.055022872335109  -0.190212128003969
  -0.169436992810515  -0.206948615170099
   0.150088939478557   0.090560438773547
   0.256256842094403   0.598874371523965
   0.340882741111244  -0.516594535669757
   0.278186325120137  -0.547968005821968
   0.645979568807173  -0.827959899268083
   0.436535833804569   0.689957746461832
   0.268437499571141   0.341281325011944
   0.120485843238972  -0.301999984465851
   0.160365386980321  -0.202012022307596
   0.154560496835611   0.244912011403144
   0.117966782622546   0.342990099492354
   0.280465408470057  -0.526206432878627
   1.174002195932550   1.136790584581798]')

L = LinearTransform(X)
X = transform(L, X)

Nψ = 5

idx = [0 0; 0 1;  1 0; 1 1; 1 2]
coeff =  [  -0.08490347737829938;
              1.0107939341310863;
              0.8792816228540664;
             -0.2535482121318712;
              0.7462002599336889]

	Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]

	for b in Blist
		B = MultiBasis(b, 2)
		f = ExpandedFunction(B, idx, coeff)
		C = HermiteMapComponent(IntegratedFunction(f); α = 0.0)
		S = Storage(C.I.f, X)

		ψoff = evaluate_offdiagbasis(C.I.f, X)
		ψd0 = evaluate_diagbasis(C.I.f, vcat(X[1:end-1,:], zeros(1, Ne)))
		∂ψd = grad_xd_diagbasis(C.I.f, X)

		function loss_coeff(c)
		    R = (ψoff .* ψd0)*c +
		    X[end,:] .* quadgk(t-> C.I.g.(ψoff .*
		    grad_xd_diagbasis(C.I.f, vcat(X[1:end-1,:], t*X[end:end,:]))*c), 0.0, 1.0)[1]
		    δ = 1e-9
		    R .+= δ*X[end,:]
		    quad = 1/Ne*sum(x->0.5*x^2+ 0.5*log(2*π), R)
		    logterm = 1/Ne*(-1.0)*sum(log.(C.I.g.( (ψoff .* ∂ψd)*c) .+ δ))
		    return quad + logterm
		end

		Jt = loss_coeff(coeff)
		dJt = ForwardDiff.gradient(loss_coeff, coeff)

		dJ = zeros(Nψ)
		J = negative_log_likelihood!(0.0, dJ, getcoeff(C), S, C, X)

		@test norm(Jt - J)<atol
		@test norm(dJt - dJ)<atol
	end
end

# Code for optimization with Hessian
# X = randn(Nx, Ne) .* randn(Nx, Ne)
# S = Storage(H.I.f, X; hess = true);
#
# J = 0.0
# dJ = zeros(Nψ)
# d2J = zeros(Nψ, Nψ)
# hess_negative_log_likelihood!(J, dJ, d2J, coeff, S, H, X)
#
# res = Optim.optimize(Optim.only_fgh!(hess_negative_log_likelihood!($S, $H, $X)), $coeff, Optim.NewtonTrustRegion())


text:
"""
    StagedTimeStep{T,N}

T: the value type of the basis function.
N: the number of stages.
It corresponds to a time-space basis function where each time step has
intermediary stages given by the vertor c in a Butcher tableau (A,b,c)
"""
struct StagedTimeStep{SB, T, N}
	spatialBasis :: SB
	c  :: SVector{N,T}
	Δt :: T
	Nt :: Int
end

scalartype(sts :: StagedTimeStep{SB, T, N}) where {SB, T, N} = T;
temporalbasis(sts :: StagedTimeStep{SB, T, N}) where {SB, T, N} = timebasisdelta(sts.Δt, sts.Nt)


text:
module Smoothing
using SparseArrays
using LarSurf

# function eye(sz)
#     arr = spzeros(Int8, sz, sz)
#     for i=1:sz
#         arr[i,i] = 1
#     end
#     return arr
# end
#
# function setIzero!(arr)
#     for i=1:size(arr)[1]
#         arr[i,i] = 0
#     end
#     return arr
# end
#
# function getDiag(arr)
#     sz = size(arr)[1]
#     arro = zeros(eltype(arr), sz)
#     for i=1:sz
#         arro[i] = arr[i,i]
#     end
#     return arro
# end
#
# function make_full(arr)
#     arrout = Array{eltype(arr), 2}(undef, size(arr)...)
#     for i=1:size(arr)[1]
#         for j=1:size(arr)[2]
#             arrout[i,j] = arr[i,j]
#         end
#     end
#     return arrout
# end
# #neighboor vertexes
#
# """
# Smoothing based on EV
# """
# function smoothing_EV(V, EVch::SparseArrays.SparseMatrixCSC, k=0.35)
# #     Matrix(setIzero!(EVch' * EVch))
#
#     neighboorNumber = getDiag(EVch' * EVch)
#     neighboors = setIzero!(EVch' * EVch)
#
# #     targetV =  neighboors * V' ./ neighboorNumber
# #     diff = targetV - V'
# #     newV = (V' + k * diff)'
#
#     targetV = V * neighboors ./ neighboorNumber'
#     # println("targetV shape: ", size(targetV))
#     diff = targetV - V
#     # println("diff shape: ", size(diff))
#     newV = V + k * diff
#
#     return make_full(newV)
# end

function eye(sz)
    arr = spzeros(Int8, sz, sz)
    for i=1:sz
        arr[i,i] = 1
    end
    return arr
end

function setIzero!(arr)
    for i=1:size(arr)[1]
        arr[i,i] = 0
    end
    return arr
end

function getDiag(arr)
    sz = size(arr)[1]
    arro = zeros(eltype(arr), sz)
    for i=1:sz
        arro[i] = arr[i,i]
    end
    return arro
end

function make_full(arr)
    arrout = Array{eltype(arr), 2}(undef, size(arr)...)
    for i=1:size(arr)[1]
        for j=1:size(arr)[2]
            arrout[i,j] = arr[i,j]
        end
    end
    return arrout
end
#neighboor vertexes

function smoothing_EV(V, EVch::SparseMatrixCSC, k=0.35)
#     Matrix(setIzero!(EVch' * EVch))

	EVchPow = EVch' * EVch
    neighboorNumber = getDiag(EVchPow)
    neighboors = setIzero!(EVchPow)
	dropzeros!(neighboors)

#     targetV =  neighboors * V' ./ neighboorNumber
#     diff = targetV - V'
#     newV = (V' + k * diff)'

    targetV = (V * neighboors) ./ neighboorNumber'
    # println("targetV shape: ", size(targetV))
    diff = targetV - V
	# @info "sm V orig   $(V[:,1:5])"
	# @info "sm V target $(targetV[:,1:5])"
	# @info "sm V diff   $(diff[:,1:5])"
    # println("diff shape: ", size(diff))
    newV = V + k * diff

    return make_full(newV)
end

"""
Get EV array from FV of surface of quads. Filtration step does not work for non surface
objects. But implementation is prepared for triangles
"""
function get_EV_quads(FV::Array{Array{Int64,1},1}; return_unfiltered=false, data=nothing)
	@debug "getFV " FV
	# EV = []
	# for f in FV
	# 	push!(EV, [[f[1],f[2]],[f[3],f[4]],  [f[1],f[3]],[f[2],f[4]]])
	# 	# push!(EV, [[f[1],f[2]],[f[3],f[4]],  [f[1],f[3]],[f[2],f[4]]])
	# end
	if size(FV[1],1) == 4

		# couples = [[1,2], [3,4], [1,3], [2,4]]
		couples = [[1,2], [2,3], [3,4], [4,1]]
	else
		couples = [[1,2], [2,3], [3,1]]
	end
	szc = size(couples,1)

	ev_time = @elapsed EV = reshape([sort([f[couples[c][1]], f[couples[c][2]]]) for f in FV, c=1:szc],:)
	# doubleedges = Base.sort(mycat(EV))
	if return_unfiltered
		if data != nothing
			data["smoothing EV construction time [s]"] = ev_time
		end
		return EV
	end
    ev_unique_time = @elapsed EVnew = collect(Set(EV))
	if data != nothing
		data["smoothing EV construction time [s]"] = ev_time
		data["smoothing make EV unique time [s]"] = ev_unique_time
	end
	return EVnew
	# doubleedges = Base.sort(EV)
	# doubleedges = convert(LarSurf.Lar.Cells, doubleedges)
	# newEV = [doubleedges[k] for k=1:2:length(doubleedges)]
	# return newEV
end

function get_EV_quads2(FV::Array{Array{Int64,1},1}; return_unfiltered=false)
	# @info "getFV " FV
	# EV = []
	# for f in FV
	# 	push!(EV, [[f[1],f[2]],[f[3],f[4]],  [f[1],f[3]],[f[2],f[4]]])
	# 	# push!(EV, [[f[1],f[2]],[f[3],f[4]],  [f[1],f[3]],[f[2],f[4]]])
	# end
	if size(FV[1],1) == 4
		couples = [[1,2], [3,4], [1,3], [2,4]]
	else
		couples = [[1,2], [2,3], [3,1]]
	end
	EV = reshape([sort([f[couples[c][1]], f[couples[c][2]]]) for f in FV, c=1:size(couples,1)],:)
	# doubleedges = Base.sort(mycat(EV))
	if return_unfiltered
		return EV
	end
	doubleedges = Base.sort(EV)
	doubleedges = convert(LarSurf.Lar.Cells, doubleedges)
	newEV = [doubleedges[k] for k=1:2:length(doubleedges)]
	return newEV
end

"""
Smoothing with defined k. Works for quads or triangles.
"""
function smoothing_FV(V::Array, FV::Array{Array{Int64,1},1}, k=0.35, n=1)
	@info "smoothing V by FV, size(V) = $(size(V)), size(FV) = $(size(FV))"

	EV = LarSurf.Smoothing.get_EV_quads(FV)
	# kEV = LarSurf.characteristicMatrix(aEV, size(bigV)[2])
	# kEV = LarSurf.characteristicMatrix(aEV, size(V)[2])
	aEV = LarSurf.ll2array(EV)
	kEV = LarSurf.characteristicMatrix(aEV, size(V)[2])
	# kEV = Lar.characteristicMatrix(EV)
	# @info "computing new V" @__MODULE__
	# @info "computing new V" @__FILE__
	newV = V
	for i=1:n
		newV = smoothing_EV(newV, kEV, k)
	end
	return newV
end

function smoothing_FV_taubin(V::Array, FV::Array{Array{Int64,1},1}, k1=0.35, k2=-0.1, n=1, data=nothing)
	@info "smoothing V by FV, size(V) = $(size(V)), size(FV) = $(size(FV))"
	if size(FV[1],1) == 4
		@info "FV are quads"
	else
		@info "FV are triangles"
	end
	EV = get_EV_quads(FV, data=data)

	# LarSurf
	aEV = LarSurf.ll2array(EV)

	# kEV = LarSurf.characteristicMatrix(aEV, size(bigV)[2])
	kEV = LarSurf.characteristicMatrix(aEV, size(V)[2])
	# kEV = Lar.characteristicMatrix(EV)
	@info "computing new V"
	newV = V
	taubin_time = @elapsed for i=1:n
		@info "taubin iteration $(i)"
		newV = smoothing_EV(newV, kEV, k1)
		newV = smoothing_EV(newV, kEV, k2)
	end

	if data != nothing
		data["smoothing taubin iterations time [s]"] = taubin_time
	end
	return newV
end

# """
# Iterative smoothing
# """
# function smoothing_FV(V::Array, FV::Array{Array{Int64,1},1}, k::Real, n::Integer)
# 	Vtmp = V
# 	for i=1:n
# 		Vtmp = smoothing_FV(Vtmp, FV, k)
# 	end
# 	return Vtmp
#
# end
#
# """
# Iterative smoothing
# """
# function smoothing_FV_taubin(V::Array, FV::Array{Array{Int64,1},1}, k1::Real, k2::Real, n::Integer)
# 	Vtmp = V
# 	for i=1:n
# 		Vtmp = smoothing_FV(Vtmp, FV, k1)
# 		Vtmp = smoothing_FV(Vtmp, FV, k2)
# 	end
# 	return Vtmp
#
# end

end


text:
"""
    Oscillator{T<:Float64}

Custom type to represent a SDOF oscillator. The type has two fields:
  - `f_n` is the natural frequency of the oscillator
  - `ζ_n` is the damping ratio

# Examples
```julia-repl
    sdof = Oscillator( 1.0, 0.05 )
```
"""
struct Oscillator{T<:Float64}
    f_n::T
    ζ_n::T
end

"""
    Oscillator(f_n)

Default initializer setting the damping ratio to 5% of critical
  - `f_n` is the natural frequency of the oscillator (in Hz)

# Examples
```julia-repl
    f_n = 2.0
    sdof = Oscillator(f_n)
```
"""
Oscillator(f_n::T) where {T<:Float64} = Oscillator(f_n, 0.05)

"""
	period(sdof::Oscillator)

Natural period (s) of the sdof Oscillator.
"""
function period(sdof::Oscillator)
	return 1.0/sdof.f_n
end

@doc raw"""
    transfer(f::T, sdof::Oscillator) where {T<:Real}

Compute the modulus of the transfer function for a SDOF system.

The transfer function is defined as:
```math
|H(f,f_n,\zeta_n)| = \frac{1}{\sqrt{ \left(1 - \beta^2 \right)^2 + \left(2\zeta_n\beta\right)^2 }}
```
where ``\beta`` is the tuning ratio defined by ``f/f_n``.

# Examples
```julia-repl
    f = 2.0
    sdof = Oscillator(1.0, 0.05)
    Hf = transfer(f, sdof)
```

See also: [`squared_transfer`](@ref)
"""
function transfer(f, sdof::Oscillator)
    # tuning ratio
    β = f / sdof.f_n
    return 1.0 / sqrt( (1.0 - β^2)^2 + (2sdof.ζ_n*β)^2 )
end

"""
	squared_transfer(f, sdof::Oscillator)

Compute the square of the transfer function for a SDOF system, `sdof`, at frequency `f`.

# Examples
```julia-repl
	f = 2.0
	# create sdof with natural frequency f_n=1.0 and damping ζ=0.05
	sdof = Oscillator( 1.0, 0.05 )
	Hf2 = squared_transfer( f, sdof )
```

See also: [`transfer`](@ref)
"""
function squared_transfer(f, sdof::Oscillator)
    # tuning ratio
    β = f / sdof.f_n
    return 1.0 / ( (1.0 - β^2)^2 + (2sdof.ζ_n*β)^2 )
end

"""
    transfer(f::Vector{T}, sdof::Oscillator) where T<:Real

Computes the modulus of the transfer function of a SDOF for a vector of frequencies
  - `f::Vector` is the vector of frequencies
  - `sdof::Oscillator` is the oscillator instance

# Examples
```julia-repl
  f = collect(range(0.1, stop=10.0, step=0.01))
  sdof = Oscillator(1.0)
  Hf = transfer(f, sdof)
```
"""
function transfer(f::Vector{T}, sdof::Oscillator) where T<:Real
    # tuning ratio
    Hf = similar(f)
    for i in 1:length(f)
        @inbounds Hf[i] = transfer(f[i], sdof)
    end
    return Hf
end

"""
    transfer!(Hf::Vector{T}, f::Vector{T}, sdof::Oscillator) where T<:Real

Computes the modulus of the transfer function of a SDOF for a vector of frequencies in place
	- `Hf::Vector` is the pre-allocated vector into which the results are stored
    - `f::Vector` is the vector of frequencies
    - `sdof::Oscillator` is the oscillator instance

# Examples
```julia-repl
    f = collect(range(0.1, stop=10.0, step=0.01))
    sdof = Oscillator(1.0)
    Hf = similar(f)
    transfer!(Hf, f, sdof)
```
"""
function transfer!(Hf::Vector{T}, f::Vector{T}, sdof::Oscillator) where T<:Real
	for i in 1:length(f)
		@inbounds Hf[i] = transfer(f[i],sdof)
	end
	return nothing
end

"""
    squared_transfer!(Hf2::Vector{T}, f::Vector{T}, sdof::Oscillator) where T<:Real

Computes the square of the modulus of the transfer function of a SDOF for a vector of frequencies in place:
	- `Hf2::Vector` is the pre-allocated vector into which the results are stored
    - `f::Vector` is the vector of frequencies
    - `sdof::Oscillator` is the oscillator instance
Inputs derive from the `Real` type and so are differentiable.

# Examples
```julia-repl
    f = collect(range(0.1, stop=10.0, step=0.01))
    sdof = Oscillator(1.0)
    Hf2 = similar(f)
    squared_transfer!(Hf2, f, sdof)
```
"""
function squared_transfer!(Hf2::Vector{T}, f::Vector{T}, sdof::Oscillator) where T<:Real
	for i in 1:length(f)
		@inbounds Hf2[i] = squared_transfer(f[i],sdof)
	end
	return nothing
end

"""
    squared_transfer!(Hf2::Vector, f::Vector, sdof::Oscillator)

Computes the square of the modulus of the transfer function of a SDOF for a vector of frequencies in place:
	- `Hf2::Vector` is the pre-allocated vector into which the results are stored
    - `f::Vector` is the vector of frequencies
    - `sdof::Oscillator` is the oscillator instance

# Examples
```julia-repl
    f = collect(range(0.1, stop=10.0, step=0.01))
    sdof = Oscillator(1.0)
    Hf2 = similar(f)
    squared_transfer!(Hf2, f, sdof)
```
"""
function squared_transfer!(Hf2::Vector, f::Vector, sdof::Oscillator)
	for i in 1:length(f)
		Hf2[i] = squared_transfer(f[i],sdof)
	end
	return nothing
end


text:
# ------------------------------------
# set volume
# ------------------------------------
function set_volume(nodes, cellxmax, cellymax, volume)
    for j in 1:cellymax
        for i in 1:cellxmax
            vec_r1x = nodes[i+1,j+1,1] - nodes[i,j,1]
            vec_r1y = nodes[i+1,j+1,2] - nodes[i,j,2]
            vec_r2x = nodes[i,j+1,1] - nodes[i+1,j,1]
            vec_r2y = nodes[i,j+1,2] - nodes[i+1,j,2]

            volume[i,j] = abs(vec_r1x*vec_r2y - vec_r1y*vec_r2x) /2
        end
    end
    return volume
end

# ------------------------------------
# set cellcenter
# ------------------------------------
function set_cellcenter(cellcenter, nodes, cellxmax, cellymax)
    for l in 1:2
        for j in 1:cellymax
            for i in 1:cellxmax
                xl = 0.25*(nodes[i+1,j+1,l] + nodes[i+1,j,l] + nodes[i,j+1,l] + nodes[i,j,l])
                cellcenter[i,j,l] = xl
            end
        end
    end
    return cellcenter
end

# ------------------------------------
# set dx and dy for local time stepping
# ------------------------------------
function set_dx_lts(dx, dy, nodes, cellxmax, cellymax, icell)
    # define at cell boundaries
    # dx = zeros(cellxmax+1, cellymax)
    # dy = zeros(cellxmax, cellymax+1)

    for j in 1+icell:cellymax -icell
        for i in 1+icell:cellxmax+1 -icell
            x1 = nodes[i,j,1]
            y1 = nodes[i,j,2]
            x2 = nodes[i,j+1,1]
            y2 = nodes[i,j+1,2]

            if (x2 - x1) == 0.0
                a = 0.0
                b  = y1 - a*x1
            else
                a = (y2-y1) / (x2-x1)
                b = y1 - a*x1
            end

            ccx = 0.125 * (nodes[i,j,1] + nodes[i+1,j,1] + nodes[i,j+1,1] + nodes[i+1,j+1,1])
            ccy = 0.125 * (nodes[i,j,2] + nodes[i+1,j,2] + nodes[i,j+1,2] + nodes[i+1,j+1,2])

            dx1 = 0.0
            if a == 0.0
                dx1 = abs(ccx - b)
            else
                dx1 = abs(-a*ccx + ccy -b)/abs(a)
            end

            ccx = 0.125 * (nodes[i-1,j,1] + nodes[i,j,1] + nodes[i-1,j+1,1] + nodes[i,j+1,1])
            ccy = 0.125 * (nodes[i-1,j,2] + nodes[i,j,2] + nodes[i-1,j+1,2] + nodes[i,j+1,2])

            dx2 = 0.0
            if a == 0.0
                dx2 = abs(ccx - b)
            else
                dx2 = abs(-a*ccx + ccy -b)/abs(a)
            end
            dx[i,j] = 0.5 * (dx1+dx2)
        end
    end

    for j in 1+icell:cellymax+1 -icell
        for i in 1+icell:cellxmax -icell
            x1 = nodes[i,j,1]
            y1 = nodes[i,j,2]
            x2 = nodes[i+1,j,1]
            y2 = nodes[i+1,j,2]

            if (x2 - x1) == 0.0
                a = 0.0
                b  = y1 - a*x1
            else
                a  = (y2-y1) / (x2-x1)
                b  = y1 - a*x1
            end

            ccx = 0.125 * (nodes[i,j,1] + nodes[i+1,j,1] + nodes[i,j+1,1] + nodes[i+1,j+1,1])
            ccy = 0.125 * (nodes[i,j,2] + nodes[i+1,j,2] + nodes[i,j+1,2] + nodes[i+1,j+1,2])

            dy1 = 0.0
            if a == 0.0
                dy1 = abs(ccy - b)
            else
                dy1 = abs(-a*ccx + ccy -b)/abs(a)
            end

            ccx = 0.125 * (nodes[i,j,1] + nodes[i+1,j,1] + nodes[i,j-1,1] + nodes[i+1,j-1,1])
            ccy = 0.125 * (nodes[i,j,2] + nodes[i+1,j,2] + nodes[i,j-1,2] + nodes[i+1,j-1,2])

            dy2 = 0.0
            if a == 0.0
                dy2 = abs(ccy - b)
            else
                dy2 = abs(-a*ccx + ccy -b)/abs(a)
            end
            dy[i,j] = 0.5 * (dy1+dy2)
        end
    end

    return dx, dy
end

# ------------------------------------
# set dtau for local time stepping
# ------------------------------------
function set_lts(dtau, lambda_facex, lambda_facey, Qbase, cellxmax, cellymax, mu, dx, dy,
                vecAx, vecAy, volume, specific_heat_ratio, cfl, icell)
    g = specific_heat_ratio

    for j in 1+icell:cellymax -icell
        for i in 1+icell:cellxmax+1 -icell
            rho_av = 0.5 * (Qbase[i,j,1] + Qbase[i-1,j,1])
            u_av   = 0.5 * (Qbase[i,j,2] + Qbase[i-1,j,2])
            v_av   = 0.5 * (Qbase[i,j,3] + Qbase[i-1,j,3])
            mu_av  = 0.5 * (   mu[i,j]   +    mu[i-1,j]  )

            ap = (g * Qbase[  i,j,4] / Qbase[  i,j,1])^0.5
            am = (g * Qbase[i-1,j,4] / Qbase[i-1,j,1])^0.5
            a_av = 0.5 * (ap + am)

            U   = u_av*vecAx[i,j,1] + v_av*vecAx[i,j,2]
            lambda_facex[i,j] = abs(U) + a_av + 2*mu_av/(rho_av*dx[i,j])
        end
    end

    for j in 1+icell:cellymax+1 -icell
        for i in 1+icell:cellxmax -icell
            rho_av = 0.5 * (Qbase[i,j,1] + Qbase[i,j-1,1])
            u_av   = 0.5 * (Qbase[i,j,2] + Qbase[i,j-1,2])
            v_av   = 0.5 * (Qbase[i,j,3] + Qbase[i,j-1,3])
            mu_av  = 0.5 * (   mu[i,j]   +    mu[i,j-1]  )

            ap = (g * Qbase[  i,j,4] / Qbase[  i,j,1])^0.5
            am = (g * Qbase[i,j-1,4] / Qbase[i,j-1,1])^0.5
            a_av = 0.5 * (ap + am)

            V   = u_av*vecAy[i,j,1] + v_av*vecAy[i,j,2]
            lambda_facey[i,j] = abs(V) + a_av + 2*mu_av/(rho_av*dy[i,j])
        end
    end
    for j in 1+icell:cellymax-icell
        for i in 1+icell:cellxmax-icell
            a1 = lambda_facex[  i,  j]
            a2 = lambda_facex[i+1,  j]
            a3 = lambda_facey[  i,  j]
            a4 = lambda_facey[  i,j+1]
            lmax = maximum([a1,a2,a3,a4])

            dtau[i,j] = cfl * volume[i,j] / lmax
        end
    end
    return dtau
end

# ------------------------------------
# set viscosity by Sutherland's formula
# https://cattech-lab.com/science-tools/sutherland/
# ------------------------------------
function set_mu(mu, Qbase, cellxmax, cellymax, specific_heat_ratio, Rd)
    mu0 = 1.82e-5     # Reference Viscosity, Pa s
    T0  = 293.15      # Reference Temperature, K
    C   = 117         # Sutherland's constant, K
    for i in 1:cellxmax
        for j in 1:cellymax
            # T = p/(rho Rd )
            T = Qbase[i,j,4]/(Qbase[i,j,1]*Rd)
            mu[i,j] = mu0 * (T/T0)^1.5 * (T0+C)/(T+C)
        end
    end
    return mu
end

# ------------------------------------
# set thermal Conductivity by Sutherland's formula
# https://doi.org/10.11357/jsam1937.37.694
# ------------------------------------
function set_lambda(lambda, Qbase, cellxmax, cellymax, mu, specific_heat_ratio, Rd)
    lam0 = 22.3*10^(-3)  # Reference thermal Conductivity, W/(m K)
    T0   = 273.15        # Reference Temperature, K
    C    = 125           # Sutherland's constant, K

    for j in 1:cellymax
        for i in 1:cellxmax
            T = Qbase[i,j,4]/(Qbase[i,j,1]*Rd)
            lambda[i,j] = lam0*((T0+C)/(T+C))*(T/T0)^1.5
        end
    end

    return lambda
end

# ------------------------------------
# set Minf for AUSM+up (don't use)
# ------------------------------------
function set_Minf(bdcon, specific_heat_ratio, Rd)
    rho = 0
    u   = 0
    v   = 0
    p   = 0

    for i in 1:4
        if Int(bdcon[i][1]) == 0 || Int(bdcon[i][1]) == 5
            rho = bdcon[i][2]
            u   = bdcon[i][3]
            v   = bdcon[i][4]
            p   = bdcon[i][5]
        elseif Int(bdcon[i][1]) == 6
            rho = bdcon[i][2]
            u   = bdcon[i][3]
            v   = bdcon[i][4]
            T   = bdcon[i][8]
            p   = rho*Rd*T
        end
    end

    a = (specific_heat_ratio * p / rho)^0.5
    u = (u^2 + v^2)^0.5
    M = u/a

    return M
end

# ------------------------------------
# Conversion from primitive variables to conserved variables
# ------------------------------------
function base_to_conservative(Qbase, Qcon, cellxmax, cellymax, specific_heat_ratio)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """

    for j in 1:cellymax
        for i in 1:cellxmax
            Qcon[i,j,1] = Qbase[i,j,1]
            Qcon[i,j,2] = Qbase[i,j,1]*Qbase[i,j,2]
            Qcon[i,j,3] = Qbase[i,j,1]*Qbase[i,j,3]
            Qcon[i,j,4] = Qbase[i,j,4]/(specific_heat_ratio-1)+Qbase[i,j,1]*(Qbase[i,j,2]^2+Qbase[i,j,3]^2)/2
        end
    end
    return Qcon
end

# ------------------------------------
# Conversion from conserved variables to primitive variables
# ------------------------------------
function conservative_to_base(Qbase, Qcon, cellxmax, cellymax, specific_heat_ratio)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """

    for j in 1:cellymax
        for i in 1:cellxmax
            Qbase[i,j,1] = Qcon[i,j,1]
            Qbase[i,j,2] = Qcon[i,j,2]/Qcon[i,j,1]
            Qbase[i,j,3] = Qcon[i,j,3]/Qcon[i,j,1]
            Qbase[i,j,4] = (Qcon[i,j,4]-Qcon[i,j,1]*(Qbase[i,j,2]^2+Qbase[i,j,3]^2)/2)*(specific_heat_ratio-1)
        end
    end
    return Qbase
end

# ------------------------------------
# Conversion from Cartesian coordinate system to general coordinate system
# ------------------------------------
function setup_Qcon_hat(Qcon, Qcon_hat, cellxmax, cellymax, volume, nval)
    for l in 1:nval
        for j in 1:cellymax
            for i in 1:cellxmax
                Qcon_hat[i,j,l] = Qcon[i,j,l] * volume[i,j]
            end
        end
    end
    return Qcon_hat
end

# ------------------------------------
# Conversion from general coordinate system to Cartesian coordinate system
# ------------------------------------
function Qhat_to_Q(Qcon, Qcon_hat, cellxmax, cellymax, volume, nval)
    for l in 1:nval
        for j in 1:cellymax
            for i in 1:cellxmax
                Qcon[i,j,l] = Qcon_hat[i,j,l] / volume[i,j]
            end
        end
    end
    return Qcon
end

# ------------------------------------
# cal average convective value
# ------------------------------------
function cal_Qave(Qbase, Qbase_ave, cellxmax, cellymax, nval)
    for l in 1:nval
        for j in 1:cellymax
            for i in 1:cellxmax
                Qbase_ave[i,j,l] += Qbase[i,j,l]
            end
        end
    end
    return Qbase_ave
end


text:
"""
    F, C = maxflow(G, s, t)

G is the adjacency matrix of an undirected graph, with edge weights equal to capacity.
Computes the maximum s-t flow, F, and the minimum cut, C.

F[a,b] is flow from a to b, and F[b,a] = -F[a,b].
C is the set reachable from s using unsaturated edges.

Uses Dinic's algorithm.  It would probably be easy to make this work for directed graphs.
"""
function maxflow(G::SparseMatrixCSC{Tv, Ti}, s::Int64, t::Int64) where {Tv, Ti}

  n = max(G.n, G.m)

  # initialize the capacity matrix
  C = copy(G)
  backInd = backIndices(G)

  Q = zeros(Ti, n)
  inQ = zeros(Bool, n)
  #=
   prev[i][1] = the vertex preceding vertex i in the BFS
   prev[i][2] = the edge index of prev[i][1]; nbri(G, prev[i][1], prev[i][2]) = i
  =#
  prev = [(0,0) for i in 1:n]
  totalflow = 0

  improve = true
  while improve
    improve = false

    # initialize the queue at the current step
    left = 0
    right = 1
    Q[right] = s
    inQ[s] = true

    # do bfs
    while left < right
      left = left + 1
      u = Q[left]
      if (u == t)
        continue
      end

      for i in 1:deg(G, u)
        v = nbri(G, u, i)

        if inQ[v] == false && weighti(C, u, i) > 0
          right = right + 1
          Q[right] = v
          inQ[v] = true
          prev[v] = (u, i)
        end
      end
    end

    # do block flow
    if inQ[t] == true
      improve = true

      for i in 1:deg(G, t)
        # ignore a neighbor that isn't part of the BF-tree
        if inQ[nbri(G, t, i)] == false
          continue
        end

        # go through every neighbor of t and see the increase in flow
        prev[t] = (nbri(G, t, i), weighti(backInd, t, i))
        u = t

        currentflow = typemax(Tv)
        while u != s
          prevU = prev[u][1]
          indexPrevU = prev[u][2]
          currentflow = min(currentflow, weighti(C, prevU, indexPrevU))
          u = prevU
        end

        if currentflow != 0
          #=
            prevU is the vertex immmediately before u on the path from s to t
            indexU is such that nbri(G, u, indexU) = prevU
            indexPrevU is such that nbri(G, prevU, indexPrevU) = u
            capPrevUU and capUPrevU have similar meanings
          =#
          u = t
          while u != s
            prevU = prev[u][1]
            indexPrevU = prev[u][2]
            indexU = weighti(backInd, prevU, indexPrevU)

            capPrevUU = weighti(C, prevU, indexPrevU)
            capUPrevU = weighti(C, u, indexU)

            setValue(C, prevU, indexPrevU, capPrevUU - currentflow)
            setValue(C, u, indexU, capUPrevU + currentflow)

            u = prevU
          end
        end

        totalflow = totalflow + currentflow
      end

    end

    # reset the queue
    for i in 1:right
      u = Q[i]

      Q[i] = 0
      inQ[u] = false
      prev[u] = (0,0)
    end
    left = right = 0

  end

  # compute the min-cut
  left = 0
  right = 1
  Q[right] = s
  inQ[s] = true

  while left < right
    left = left + 1

    u = Q[left]
    for i in 1:deg(G, u)
      v = nbri(G, u, i)

      if inQ[v] == 0 && weighti(C, u, i) > 0
        right = right + 1
        Q[right] = v
        inQ[v] = 1
      end
    end
  end

  cut = Q[1:right]
  flow = C - G

  return flow, cut
end


text:
#=
145 is a curious number, as 1! + 4! + 5! = 1 + 24 + 120 = 145.

Find the sum of all numbers which are equal to the sum of the factorial of their digits.

Note: as 1! = 1 and 2! = 2 are not sums they are not included.
=#
function calc()
  total = 0
  for i in 3:9_999_999
    if sum(map((x)->factorial(x),digits(i))) == i
      total += i
    end
  end
  total
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
function ResidualOnFaces_v2_GPU!(F, Mesh_bc, Mesh_f2e, Mesh_dA_f, Mesh_n_x_f, Mesh_n_y_f, Mesh_vole_f, Mesh_vole, Mesh_e2f, Mesh_dA, Mesh_n_x, Mesh_n_y, Th, Te, Tneu, qx, qy, ae, be, ze, tau, mesh_nf, mesh_nf_el)

    idof = (blockIdx().x-1) * blockDim().x + threadIdx().x

    if idof <= mesh_nf
        # Identify BC faces
        bc = Mesh_bc[idof]
        Xi = 0.0 + (bc==2)*1.0   # Neumann
        # ti = Tneu[idof]        # Neumann

        # Element 1
        yes1  = Mesh_f2e[idof,1] > 0
        iel1  = (yes1==1) * Mesh_f2e[idof,1] + (yes1==0)*1 # if the neigbouring element is out then set index to 1
        # dAi1  = (yes1==1) * Mesh_dA_f[idof,1]
        # ni_x1 = (yes1==1) * Mesh_n_x_f[idof,1]
        # ni_y1 = (yes1==1) * Mesh_n_y_f[idof,1]
        # taui1 = (yes1==1) * tau #StabParam(tau,dAi1,Mesh_vole_f[idof,1],mesh_type)
        # >>>>>>>> Here they are computed on the fly
        Tel1   = (yes1==1) * be[iel1]/ae[iel1]
        qxel1  = (yes1==1) * -ze[iel1,1]/Mesh_vole[iel1]
        qyel1  = (yes1==1) * -ze[iel1,2]/Mesh_vole[iel1]
        for ifac=1:mesh_nf_el
            # Face
            nodei = Mesh_e2f[iel1,ifac]
            bcf   = Mesh_bc[nodei]
            dAi   = Mesh_dA[iel1,ifac]
            ni_x  = Mesh_n_x[iel1,ifac]
            ni_y  = Mesh_n_y[iel1,ifac]
            taui  = tau #StabParam(tau,dAi,Mesh_vole[iel1],mesh_type)      # Stabilisation parameter for the face
            # Assemble
            Tel1  = Tel1  + (yes1==1) * (bcf!=1) * dAi*taui*Th[Mesh_e2f[iel1,ifac]]/ae[iel1]
            qxel1 = qxel1 - (yes1==1) * (bcf!=1) * dAi*ni_x*Th[Mesh_e2f[iel1,ifac]]/Mesh_vole[iel1]
            qyel1 = qyel1 - (yes1==1) * (bcf!=1) * dAi*ni_y*Th[Mesh_e2f[iel1,ifac]]/Mesh_vole[iel1]
        end
        # F1 = (bc!=1) * (yes1==1) * dAi1* (ni_x1*qxel1 + dAi1*ni_y1*qyel1 + dAi1*taui1*Tel1 - dAi1*taui1*Th[idof] + dAi1*Xi*Tneu[idof])
        F1 = (bc!=1) * (yes1==1) * Mesh_dA_f[idof,1] * (Mesh_n_x_f[idof,1]*qxel1 + Mesh_n_y_f[idof,1]*qyel1 + tau*Tel1 - tau*Th[idof] + Xi*Tneu[idof])
        # Element 2
        yes2  = Mesh_f2e[idof,2] > 0
        iel2  = (yes2==1) * Mesh_f2e[idof,2] + (yes2==0)*1 # if the neigbouring element is out then set index to 1
        dAi2  = (yes2==1) * Mesh_dA_f[idof,2]
        # ni_x2 = (yes2==1) * Mesh_n_x_f[idof,2]
        # ni_y2 = (yes2==1) * Mesh_n_y_f[idof,2]
        # taui2 = (yes2==1) * tau #StabParam(tau,dAi2,Mesh_vole_f[idof,2],mesh_type)
        # >>>>>>>> Here they are computed on the fly
        Tel2   = (yes2==1) * be[iel2]/ae[iel2]
        qxel2  = (yes2==1) * -ze[iel2,1]/Mesh_vole[iel2]
        qyel2  = (yes2==1) * -ze[iel2,2]/Mesh_vole[iel2]
        for ifac=1:mesh_nf_el
            # Face
            nodei = Mesh_e2f[iel2,ifac]
            bcf   = Mesh_bc[nodei]
            dAi   = Mesh_dA[iel2,ifac]
            ni_x  = Mesh_n_x[iel2,ifac]
            ni_y  = Mesh_n_y[iel2,ifac]
            taui  = tau #StabParam(tau,dAi,Mesh_vole[iel2],mesh_type)      # Stabilisation parameter for the face
            # Assemble
            Tel2  = Tel2  + (yes2==1) * (bcf!=1) *  dAi*taui*Th[Mesh_e2f[iel2, ifac]]/ae[iel2]
            qxel2 = qxel2 - (yes2==1) * (bcf!=1) *  dAi*ni_x*Th[Mesh_e2f[iel2, ifac]]/Mesh_vole[iel2]
            qyel2 = qyel2 - (yes2==1) * (bcf!=1) *  dAi*ni_y*Th[Mesh_e2f[iel2, ifac]]/Mesh_vole[iel2]
        end
        # F2 = (bc!=1) * (yes2==1) * (dAi2*ni_x2*qxel2 + dAi2*ni_y2*qyel2 + dAi2*taui2*Tel2 - dAi2*taui2*Th[idof] + dAi2*Xi*Tneu[idof])
        F2 = (bc!=1) * (yes2==1) * Mesh_dA_f[idof,2] * (Mesh_n_x_f[idof,2]*qxel2 + Mesh_n_y_f[idof,2]*qyel2 + tau*Tel2 - tau*Th[idof] + Xi*Tneu[idof])

        # Contributions from the 2 elements
        F[idof] = F1 + F2
    end
    return
end

function Update_F_GPU!(F, Th_PT, F0, dTdtau, dmp, mesh_nf)

    idof = (blockIdx().x-1) * blockDim().x + threadIdx().x

    if idof <= mesh_nf
        F[idof]     = (1.0 - dmp)*F0[idof] + F[idof]
        Th_PT[idof] = Th_PT[idof] + dTdtau*F[idof]
        F0[idof]    = F[idof]
    end
    return
end

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

function ResidualOnFaces(mesh, Th, Te, qx, qy, tau)
    F = zeros(mesh.nf)
    @avx for idof=1:mesh.nf

        # Identify BC faces
        bc = mesh.bc[idof]

        # Element 1
        yes1  = mesh.f2e[idof,1]>0
        iel1  = (yes1==1) * mesh.f2e[idof,1]  + (yes1==0) * 1 # if the neigbouring element is out then set index to 1
        dAi1  = (yes1==1) * mesh.dA_f[idof,1]
        ni_x1 = (yes1==1) * mesh.n_x_f[idof,1]
        ni_y1 = (yes1==1) * mesh.n_y_f[idof,1]
        taui1 = (yes1==1) * StabParam(tau,dAi1,mesh.vole_f[idof,1],mesh.type)
        Tel1  = (yes1==1) * Te[iel1]
        qxel1 = (yes1==1) * qx[iel1]
        qyel1 = (yes1==1) * qy[iel1]
        F1    = (bc!=1) * (yes1==1) * (dAi1*ni_x1*qxel1 + dAi1*ni_y1*qyel1 + dAi1*taui1*Tel1 - dAi1*taui1*Th[idof])

        # Element 2
        yes2  = mesh.f2e[idof,2]>0
        iel2  = (yes2==1) * mesh.f2e[idof,2]   + (yes2==0) * 1 # if the neigbouring element is out then set index to 1
        dAi2  = (yes2==1) * mesh.dA_f[idof,2]
        ni_x2 = (yes2==1) * mesh.n_x_f[idof,2]
        ni_y2 = (yes2==1) * mesh.n_y_f[idof,2]
        taui2 = (yes2==1) * StabParam(tau,dAi2,mesh.vole_f[idof,2],mesh.type)
        Tel2  = (yes2==1) * Te[iel2]
        qxel2 = (yes2==1) * qx[iel2]
        qyel2 = (yes2==1) * qy[iel2]
        F2    = (bc!=1) * (yes2==1) * (dAi2*ni_x2*qxel2 + dAi2*ni_y2*qyel2 + dAi2*taui2*Tel2 - dAi2*taui2*Th[idof])

        # Contributions from the 2 elements
        F[idof] = F1 + F2
    end
    return F
end

function ComputeElementValuesFaces(mesh, uh, ae, be, ze, Tdir, tau)
    ue          = zeros(mesh.nel);
    qx          = zeros(mesh.nel);
    qy          = zeros(mesh.nel);

    for idof=1:mesh.nf #

        # Identify BC faces
        bc = mesh.bc[idof]

        # Element 1
        yes1  = mesh.f2e[idof,1]>0
        iel1  = (yes1==1) * mesh.f2e[idof,1]  + (yes1==0) * 1 # if the neigbouring element is out then set index to 1
        dAi1  = (yes1==1) * mesh.dA_f[idof,1]
        ni_x1 = (yes1==1) * mesh.n_x_f[idof,1]
        ni_y1 = (yes1==1) * mesh.n_y_f[idof,1]
        taui1 = (yes1==1) * StabParam(tau,dAi1,mesh.vole_f[idof,1],mesh.type)

        ue[iel1] += (yes1==1) * 1.0/mesh.nf_el*be[iel1]/ae[iel1]
        ue[iel1] += (bc!=1) * (yes1==1) * dAi1*taui1*uh[idof]/ae[iel1]

        qx[iel1] += (yes1==1) * 1.0/mesh.nf_el* -1.0/mesh.vole[iel1]*ze[iel1,1]
        qx[iel1] -= (bc!=1) * (yes1==1) * 1.0/mesh.vole[iel1]*dAi1*ni_x1*uh[idof]

        qy[iel1] += (yes1==1) * 1.0/mesh.nf_el* -1.0/mesh.vole[iel1]*ze[iel1,2]
        qy[iel1] -= (bc!=1) * (yes1==1) * 1.0/mesh.vole[iel1]*dAi1*ni_y1*uh[idof]

        # Element 2
        yes2  = mesh.f2e[idof,2]>0
        iel2  = (yes2==1) * mesh.f2e[idof,2]   + (yes2==0) * 1 # if the neigbouring element is out then set index to 1
        dAi2  = (yes2==1) * mesh.dA_f[idof,2]
        ni_x2 = (yes2==1) * mesh.n_x_f[idof,2]
        ni_y2 = (yes2==1) * mesh.n_y_f[idof,2]
        taui2 = (yes2==1) * StabParam(tau,dAi2,mesh.vole_f[idof,2],mesh.type)

        ue[iel2] += (yes2==1) * 1.0/mesh.nf_el*be[iel2]/ae[iel2]
        ue[iel2] += (bc!=1) * (yes2==1) * dAi2*taui2*uh[idof]/ae[iel2]

        qx[iel2] -= (yes2==1) * 1.0/mesh.nf_el* 1.0/mesh.vole[iel2]*ze[iel2,1]
        qx[iel2] -= (bc!=1) * (yes2==1) * 1.0/mesh.vole[iel2]*dAi2*ni_x2*uh[idof]

        qy[iel2] -= (yes2==1) * 1.0/mesh.nf_el* 1.0/mesh.vole[iel2]*ze[iel2,2]
        qy[iel2] -= (bc!=1) * (yes2==1) * 1.0/mesh.vole[iel2]*dAi2*ni_y2*uh[idof]

    end
    return ue, qx, qy
end


text:
nsamples, nchains, nchainsatone = 200001, 16, 1
Tmax = 2.5

addprocs(nchains)
@info "workers are $(workers())"
@everywhere any(pwd() .== LOAD_PATH) || push!(LOAD_PATH, pwd())
@everywhere using Distributed
@everywhere import MCMC_Driver
## run McMC
@time MCMC_Driver.main(opt, d, Tmax=Tmax, nsamples=nsamples, nchains=nchains, nchainsatone=nchainsatone)
rmprocs(workers())


text:
using Test
using GradientBoost.Util
using Random

@testset "utils" begin
    @testset "weighted_median" begin
        weights = [1, 1, 1, 1, 1]
        values = [1, 2, 3, 4, 5]
        @test weighted_median(weights, values) == 3

        weights = [1, 1, 1, 1, 1]
        values = [5, 4, 3, 2, 1]
        @test weighted_median(weights, values) == 3

        weights = [1, 1, 1, 1, 1]
        values = [4, 5, 2, 3, 1]
        @test weighted_median(weights, values) == 3

        weights = [1, 3, 1, 1, 1]
        values = [1, 2, 3, 4, 5]
        @test weighted_median(weights, values) == 2

        weights = [1, 1, 1, 1, 100]
        values = [1, 2, 3, 4, 5]
        @test weighted_median(weights, values) == 5

        weights = [100, 1, 1, 1, 1]
        values = [1, 2, 3, 4, 5]
        @test weighted_median(weights, values) == 1

        weights = [0.5, 0.5, 1, 1, 1]
        values = [1, 2, 3, 4, 5]
        @test weighted_median(weights, values) == 3
    end

    @testset "holdout returns proportional partitions" begin
        n = 10
        right_prop = 0.3
        (left, right) = holdout(n, right_prop)

        @test size(left, 1) == n - (n * right_prop)
        @test size(right, 1) == n * right_prop
        @test isempty(intersect(left, right))
        @test size(union(left, right), 1) == n
    end

end # module


text:
"""
    erig(b::Basis, i::Int, j::Int; alpha=1.0, scale=1.0)
    erig(b::Basis, i::Int, j::Int, n::Int, m::Int; alpha=1.0, scale=1.0)

Electron repulsion integral with Gaussian repulsion scale*exp(-α*r^2).

For erig4 `i` and `j` are indices for particale 1 and `n` and `m` for particale 2.
"""
function erig(b::AbstractBasis, i::Int, j::Int; Ve=x->exp(-x^2))
    return Ve(b[i]-b[j])
end

function erig(b::AbstractBasis, i::Int, j::Int, n::Int, m::Int; Ve=x->exp(-x^2))::Float64
    # TODO this function might not be type stable
    T = Float64   # This need to change for AD to work eg. Base.return_type()
    # Basis functions are orthogonal
    i != j && return zero(T)
    n != m && return zero(T)
    return Ve(b[i]-b[n])
end

function fock_matrix(b::AbstractBasis; Ve=x->exp(-x^2), Vn=x->-2exp(-0.5x^2))
    orbitals = initial_orbitals(b)
    return fock_matrix(b, orbitals; Ve=Ve, Vn=Vn)
end

function fock_matrix(b::AbstractBasis, orbitals::AbstractMatrix; Ve=x->exp(-x^2), Vn=x->-2exp(-0.5x^2))
    @argcheck size(orbitals) == (length(b), length(b))
    h₁ = one_electron_operator(b, Vn)
    J = coulomb_matrix(b, orbitals)
    #K = exchange_matrix(b, orbitals) # For more than 2 electrons
    return h₁ + J  # h₁ + 2*J - K
end

function fock_matrix!(f::AbstractMatrix, b::AbstractBasis, orbitals::AbstractMatrix; Ve=x->exp(-x^2), Vn=x->-2exp(-0.5x^2))
    @argcheck size(f) == size(orbitals) == (length(b), length(b))
    h₁ = one_electron_operator(b, Vn)
    J = coulomb_matrix(b, orbitals)
    #K = exchange_matrix(b, orbitals) # For more than 2 electrons
    return f .= h₁ + J  # h₁ + 2*J - K
end

function coulomb_matrix(b::AbstractBasis, orbitals::AbstractMatrix)
    C = zeros(size(orbitals))
    return coulomb_matrix!(C, b, orbitals)
end

function coulomb_matrix!(C::AbstractMatrix, b::AbstractBasis, orbitals::AbstractMatrix)
    @argcheck size(C) == size(orbitals) == (length(b), length(b))
    l = length(b)
    w = get_weight(b)

    # Two electrons in total
    ρ = orbitals[:,1] * orbitals[:,1]'
    for i in 1:l
        # C[i,j]=0 for i!=j
        # also erig(b,i,i,n,m) = 0 for n!=m
        C[i,i] = sum( n -> erig(b,i,i,n,n) * ρ[n,n] * w[n], 1:l)
        C[i,i] *= w[i]
    end
    return C
end

function exchange_matrix(b::AbstractBasis, orbitals::AbstractMatrix)
    @argcheck size(orbitals) == (length(b), length(b))
    l = length(b)
    K = zeros(l,l)
    w = get_weight(b)

    # Two electrons in total
    ρ = orbitals[:,1] * orbitals[:,1]'
    Threads.@threads for i in 1:l
        for j in 1:l
            for n in 1:l
                K[i,j] += sum( k-> erig(b, i,n,j,k) * ρ[n,k], 1:l)
            end
            K[i,j] *= w[i] * w[j]
        end
    end
    return K
end

function bracket(
        b::AbstractBasis,
        psi1::AbstractVector,
        op::AbstractMatrix,
        psi2::AbstractVector
    )
    w = get_weight(b)
    return (conj.(psi1).*w)' * op * psi2
end

function bracket(b::AbstractBasis, psi1::AbstractVector, psi2::AbstractVector)
    w = get_weight(b)
    return sum( conj.(psi1) .* w .* psi2 )
end

function kinetic_energy(b::AbstractBasis; mass=1.0)
    ∇ = derivative_matrix(b)
    g = metric_tensor(b)
    return (0.5/mass) * ∇' * g * ∇ # = 0.5 * (g∇)' * ∇
end

function one_electron_operator(b, Vn; mass=1.0)
    g = metric_tensor(b)
    V = g * diagm( Vn.(b) )
    return kinetic_energy(b; mass=mass) + V
end

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
game = basic_game(;window_visible=false)

ViZDoom.init(game)

actions = [[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]
episodes = 1
sleep_time = 1.0 / ViZDoom.DEFAULT_TICRATE

for i in 1:episodes
    println("Episode #$i")
    ViZDoom.new_episode(game)

    while !ViZDoom.is_episode_finished(game)
        state = ViZDoom.get_screen_buffer(game)
        print(size(state))
        r = ViZDoom.make_action(game, rand(actions))
        println("Reward $r")
        sleep(sleep_time)
    end
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module CyclotomicBinaryForms
using Nemo, Counts, NumberTheory

export ModuleCyclotomicBinaryForms
export is206864, F206864, I206864, L206864
export is206942, F206942, I206942, L206942
export is293654, F293654, I293654, L293654
export is296095, F296095, I296095, L296095
export V299214, L299214
export is299498, F299498, I299498, L299498
export is299733, L299733
export is299928, F299928, I299928, L299928
export is299929, F299929, I299929, L299929
export is299930, F299930, I299930, L299930
export is325143, F325143, I325143, L325143
export is325145, F325145, I325145, L325145

"""

E. Fouvry, C. Levesque, M. Waldschmidt,
[Representation of integers by cyclotomic binary forms](https://arxiv.org/pdf/1712.09019.pdf), arXiv:1712.09019 [math.NT], 2017.

* is206864, F206864, I206864, L206864, is206942, F206942, I206942, L206942, is293654, F293654, I293654, L293654, is296095, F296095, I296095, L296095, V299214, L299214, is299498, F299498, I299498, L299498, is299733, L299733, is299928, F299928, I299928, L299928, is299929, F299929, I299929, L299929, is299930, F299930, I299930, L299930, is325143, F325143, I325143, L325143, is325145, F325145, I325145, L325145
"""
const ModuleCyclotomicBinaryForms = ""

#A206942########################################################################
"""

Is ``n`` a numbers of the form ``𝚽_k(m)`` with ``k > 2`` and ``|m| > 1``
where ``𝚽_k(m)`` denotes the ``k``-th cyclotomic polynomial evaluated at ``m``.
"""
function is206942(n)
    if n < 3
        return false
    end
    R, x = PolynomialRing(ZZ, "x")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3))        # Fouvry & Levesque & Waldschmidt
    for k ∈ 3:K
        c = cyclotomic(k, x)
        for m ∈ 2:M
            n == subst(c, m) && return true
        end
    end
    return false
end

"""

Filter the integers which are A206942 and  ``<= n``.
"""
F206942(n) = filter(is206942, 1:n)

"""

Iterate over the first ``n`` integers which are A206942.
"""
I206942(n) = takeFirst(is206942, n)

"""

Return the list of the first ``n`` integers which are A206942.
"""
L206942(n) = collect(I206942(n))

#A206864########################################################################

"""

Is ``n`` a prime of the form ``𝚽_k(m)`` with ``k > 2`` and ``|m| > 1``
where ``𝚽_k(m)`` denotes the ``k``-th cyclotomic polynomial evaluated at ``m``.
"""
is206864(n) = isPrime(n) && is206942(n)

"""

Filter the integers which are A206864 and ``<= n``.
"""
F206864(n) = (j for j ∈ 1:n if is206864(j))

"""

Iterate over the first ``n`` integers which are A206864.
"""
I206864(n) = takeFirst(is206864, n)

"""

Return the list of the first ``n`` integers which are A206864.
"""
L206864(n) = collect(I206864(n))

#A299498########################################################################

"""

Is ``n`` primitively represented by a cyclotomic binary forms?
"""
function is299498(n)
    n == 0 && return false
    isPrimeTo(n, k) = gcd(fmpz(n), fmpz(k)) == fmpz(1)
    R, x = PolynomialRing(ZZ, "x")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3))  # Fouvry & Levesque & Waldschmidt
    N = QQ(n)

    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, x)
        for m ∈ 1:M, j ∈ m+1:M
            if isPrimeTo(m, j)
                N == m^e * subst(c, QQ(j, m)) && return true
            end
        end
    end
    return false
end

"""

Filter the integers which are A299498 and <= ``n``.
"""
F299498(n) = (j for j ∈ 1:n if is299498(j))

"""

Iterate over the first ``n`` integers which are A299498.
"""
I299498(n) = takeFirst(is299498, n)

"""

Return the list of the first ``n`` integers which are A299498.
"""
L299498(n) = collect(I299498(n))

"""

Count the number of cyclotomic binary forms which primitively represent ``n``.
"""
function countA299498(n)
    if n < 3
        return 0
    end
    isPrimeTo(n, k) = gcd(fmpz(n), fmpz(k)) == fmpz(1)
    R, x = PolynomialRing(ZZ, "x")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3))  # Fouvry & Levesque & Waldschmidt
    N = QQ(n)
    count = 0

    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, x)
        for m ∈ 1:M, j ∈ m+1:M
            if isPrimeTo(m, j)
                N == m^e * subst(c, QQ(j, m)) && (count += 1)
            end
        end
    end
    count
end
# [countA299498(n) for n ∈ 1:100] |> println

#A299928########################################################################

"""

Is ``n`` represented by a cyclotomic binary form f(x, y) where x and y are prime numbers and 0 < y < x?
"""
function is299928(n)
    n == 0 && return false
    R, z = PolynomialRing(ZZ, "z")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3))  # Fouvry & Levesque & Waldschmidt
    N = QQ(n)
    P(u) = (p for p ∈ u:M if isprime(fmpz(p)))
    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, z)
        for y ∈ P(2), x ∈ P(y + 1)
            N == y^e * subst(c, QQ(x, y)) && return true
        end
    end
    return false
end

"""

Filter the integers which are A299928 and <= ``n``.
"""
F299928(n) = (j for j ∈ 1:n if is299928(j))

"""

Iterate over the first ``n`` integers which are A299928.
"""
I299928(n) = takeFirst(is299928, n)

"""

Return the list of the first ``n`` integers which are A299928.
"""
L299928(n) = collect(I299928(n))

#A299929########################################################################

"""

Is ``n`` a prime represented by a cyclotomic binary form f(x, y) where x and y are prime numbers and 0 < y < x?
"""
is299929(n) = isPrime(n) && is299928(n)

"""

Filter the integers which are A299929 and <= ``n``.
"""
F299929(n) = (j for j ∈ 1:n if is299929(j))

"""

Iterate over the first ``n`` integers which are A299929.
"""
I299929(n) = takeFirst(is299929, n)

"""

Return the list of the first ``n`` integers which are A299929.
"""
L299929(n) = collect(I299929(n))

#A299930########################################################################

"""

Is ``n`` a prime represented by a cyclotomic binary form f(x, y) with x and y odd prime numbers and x > y.
"""
function is299930(n)
    !isprime(fmpz(n)) && return false
    R, z = PolynomialRing(ZZ, "z")
    K = floor(Int, 5.383 * log(n)^1.161)
    M = floor(Int, 2 * sqrt(n / 3))
    N = QQ(n)
    P(u) = (p for p ∈ u:M if isprime(fmpz(p)))
    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, z)
        for y ∈ P(3), x ∈ P(y + 2)
            N == y^e * subst(c, QQ(x, y)) && return true
        end
    end
    return false
end

"""

Filter the integers which are A299930 and <= ``n``.
"""
F299930(n) = (j for j ∈ 1:n if is299930(j))

"""

Iterate over the first ``n`` integers which are A299930.
"""
I299930(n) = takeFirst(is299930, n)

"""

Return the list of the first ``n`` integers which are A299930.
"""
L299930(n) = collect(I299930(n))

#A296095########################################################################

"""

Is ``n`` represented by a cyclotomic binary form?
"""
function is296095(n)
    n < 3 && return false
    R, z = PolynomialRing(ZZ, "z")
    N = QQ(n)
    # Bounds from Fouvry & Levesque & Waldschmidt
    logn = log(n)^1.161
    K = floor(Int, 5.383 * logn)
    M = floor(Int, 2 * (n / 3)^(1 / 2))
    k = 3
    while true
        c = cyclotomic(k, z)
        e = Int(euler_phi(fmpz(k)))
        if k == 7
            K = ceil(Int, 4.864 * logn)
            M = ceil(Int, 2 * (n / 11)^(1 / 4))
        end
        for y ∈ 2:M, x ∈ 1:y
            N == y^e * subst(c, QQ(x, y)) && return true
        end
        k += 1
        k > K && break
    end
    return false
end

"""

Filter the integers which are A296095 and <= ``n``.
"""
F296095(n) = (j for j ∈ 1:n if is296095(j))

"""

Iterate over the first ``n`` integers which are A296095.
"""
I296095(n) = takeFirst(is296095, n)

"""

Return the list of the first ``n`` integers which are A296095.
"""
L296095(n) = collect(I296095(n))

#A293654########################################################################

"""

Is ``n`` unrepresentable by a cyclotomic binary forms?
"""
function is293654(n)
    if n < 3
        return true
    end
    R, x = PolynomialRing(ZZ, "x")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3)) # Fouvry & Levesque & Waldschmidt
    N = QQ(n)
    count = 0

    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, x)
        for m ∈ 1:M, j ∈ 0:M
            if max(j, m) > 1
                N == m^e * subst(c, QQ(j, m)) && return false
            end
        end
    end
    return true
end

"""

Filter the integers which are A293654 and <= ``n``.
"""
F293654(n) = (j for j ∈ 1:n if is293654(j))

"""

Iterate over the first ``n`` integers which are A293654.
"""
I293654(n) = takeFirst(is293654, n)

"""

Return the list of the first ``n`` integers which are A293654.
"""
L293654(n) = collect(I293654(n))

#A325143########################################################################

"""

Is ``n`` a prime represented by a cyclotomic binary form?
"""
function is325143(n)
    (n < 3 || !isprime(fmpz(n))) && return false
    R, x = PolynomialRing(ZZ, "x")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3)) # Fouvry & Levesque & Waldschmidt
    N = QQ(n)

    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, x)
        for m ∈ 1:M, j ∈ 0:M
            if max(j, m) > 1
                N == m^e * subst(c, QQ(j, m)) && return true
            end
        end
    end
    return false
end

"""

Filter the integers which are A325143 and <= ``n``.
"""
F325143(n) = (j for j ∈ 1:n if is325143(j))

"""

Iterate over the first ``n`` integers which are A325143.
"""
I325143(n) = takeFirst(is325143, n)

"""

Return the list of the first ``n`` integers which are A325143.
"""
L325143(n) = collect(I325143(n))

#A325145########################################################################

"""

Is ``n`` a prime unrepresentable by a cyclotomic binary form?
"""
is325145(n) = isprime(fmpz(n)) && !is325143(n)

"""

Filter the integers which are A325145 and <= ``n``.
"""
F325145(n) = (j for j ∈ 1:n if is325145(j))

"""

Iterate over the first ``n`` integers which are A325145.
"""
I325145(n) = takeFirst(is325145, n)

"""

Return the list of the first ``n`` integers which are A325145.
"""
L325145(n) = collect(I325145(n))

#A299214########################################################################

"""

Return the number of representations of ``n`` by cyclotomic binary forms.
"""
function V299214(n)
    if n < 3
        return 0
    end
    R, x = PolynomialRing(ZZ, "x")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3))  # Fouvry & Levesque & Waldschmidt
    N = QQ(n)
    count = 0

    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, x)
        for m ∈ 1:M, j ∈ 0:M
            if max(j, m) > 1
                N == m^e * subst(c, QQ(j, m)) && (count += 1)
            end
        end
    end
    4 * count
end

"""

Return the initial list of V299214 of length len .
"""
L299214(len) = [V299214(i) for i ∈ 1:len]

################################################################################

"""

Is n a prime represented in more than one way by cyclotomic binary forms f(x,y) with x and y prime numbers and y < x?
"""
function is299733(n)
    if n < 3 || !isprime(fmpz(n))
        return false
    end
    R, x = PolynomialRing(ZZ, "x")
    K = floor(Int, 5.383 * log(n)^1.161) # Bounds from
    M = floor(Int, 2 * sqrt(n / 3)) # Fouvry & Levesque & Waldschmidt
    N = QQ(n)
    multi = 0

    for k ∈ 3:K
        e = Int(euler_phi(fmpz(k)))
        c = cyclotomic(k, x)
        for m ∈ 2:M
            if isprime(fmpz(m))
                for j ∈ m:M
                    if isprime(fmpz(j))
                        if N == m^e * subst(c, QQ(j, m))
                            multi += 1
                        end
                    end
                end
            end
        end
    end
    multi > 1
end

"""

Return the list of the first ``n`` integers which are A299733.
"""
L299733(n) = [i for i ∈ 1:n if is299733(i)]

# L299733(1000) |> println

#START-TEST-########################################################

using Test, SeqTests

function test()
    @testset "CyclotomicForms" begin
        if data_installed()

            L = [
                L206942,
                L206864,
                L299498,
                L299928,
                L299929,
                L299930,
                L299214,
                L296095,
                L293654
            ]
            SeqTest(L, 'L')
        end
    end
end

function demo()
    println("\nIntegers which are A206942 and <= 32:")
    for n ∈ 1:32
        is206942(n) && println(n)
    end

    println("\nIterate over the first 32 integers which are A206942:")
    for i ∈ I206942(32)
        print("$i, ")
    end

    println("\n\nFilter the integers which are A206942 and <= 60:")
    F206942(60) |> println

    println("\nReturn the list of the first 32 integers which are A206942:")
    L206942(32) |> println
end

"""

L299214(100) :: 1.822155 seconds (7.22 M allocations: 185.560 MiB, 36.76% gc time)
"""
function perf()
    @time L299214(100)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
##################################################
# Use this file to compute the diffusion coefficients
# Use the cluster model given in the "sources" folder

# Use the following command to get truncated quantities up to lmax = 10
# julia Compute.jl --lmax 10
##################################################

include("../sources/Main.jl")

a = 10.0 # Semi-major axis (in mpc)
j = 0.6 # Reduced angular momentum J/Jc (between 0 and 1)

println("----------------------------------------------")
println("Computing both DRRj and DRRjj at the same time")
@time djBoth, djjBoth = DRR_j_jj(a,j)
println("DRRj(a,j)  = ",djBoth)
println("DRRjj(a,j) = ",djjBoth)

println("----------------------------------------------")

println("Computing DRRjj only")
@time djj_only = DRR_jj(a,j)
println("DRRjj(a,j) = ",djj_only)
println("----------------------------------------------")


text:
# Module taken from:  https://github.com/alexyarosh/hyperbolic

using Plots
using Eirene
using Ripserer
using Statistics

# compute the persistent betti numbers of the Vietoris-Rips complex given by the distance matrix `matr`
# if mintime, maxtime and numofsteps are specified -- returns a `numofsteps x maxdim` array
# if either of the keyword arguments is not specified or is set to Inf, returns `maxdim` arrays for method=:eirene, or error for :ripser

# function bettis(matr, maxdim; mintime=-Inf, maxtime=Inf, numofsteps=Inf, method=:ripser)
#     if (method == :ripser) || (method == :Ripser)
#         if VERSION < v"0.7.0"
#             error("Ripser requires at least Julia v 0.7.0")
#         end
#         return bettis_ripser(matr, maxdim, mintime=mintime, maxtime=maxtime, numofsteps=numofsteps)
#     elseif (method == :eirene) || (method == :Eirene)
#         return bettis_eirene(matr, maxdim, mintime=mintime, maxtime=maxtime, numofsteps=numofsteps)
#     else
#         error("Method $(method) is not supported. Supported methods are: method=:eirene, method=:ripser")
#     end
# end

function bettis_ripser(matr, maxdim; mintime=-Inf, maxtime=Inf, numofsteps=Inf)
    if (mintime == -Inf) || (maxtime == -Inf) || (numofsteps == -Inf)
        error("To use Ripser, specify parameters mintime, maxtime, numofsteps")
    end
    r = ripser(matr,  dim_max = maxdim, threshold = maxtime)

    int_length = maxtime-mintime
    step_length= int_length/numofsteps
    betts = zeros(numofsteps, maxdim)

    for dim=1:maxdim
        ints = r[dim+1]
        for intl in ints
            st = Int(ceil((intl[1]-mintime)/step_length))
            if intl[2] == Inf
                fin = numofsteps
            else
                fin = Int(ceil((intl[2]-mintime)/step_length))
            end
            betts[st:fin, dim] = map(x->x+1, betts[st:fin, dim])
         end
    end
    return betts
end
#
# # Original function returns 2 different types of betti curves. If no default
# # value parameters is given, it returns vector of matrices. If num of steps is
# # given, then it return matrix maxdim x numsteps.
# function bettis_eirene(matr, maxdim; mintime=-Inf, maxtime=Inf, numofsteps=Inf)
#     c = eirene(matr, minrad = mintime, maxrad= maxtime, numrad= numofsteps, maxdim=maxdim)
#
#     int_length = maxtime-mintime
#     step_length= int_length/numofsteps
#
#     if (mintime == -Inf) || (maxtime == Inf) || (numofsteps == Inf)
#         # return [betticurve(c, dim=maxdim) for d=1:maxdim]
#         return hcat([betticurve(c, dim=d)[:,2] for d=1:maxdim]...)
#     end
#
#     betts = zeros(numofsteps, maxdim)
#     # For every dimension compute betti curve
#     for dim=1:maxdim
#         bet = betticurve(c, dim=dim)
#
#         #for every element in betti curve return betti value if index is positive
#         for i=1:size(bet,1)
#             b = bet[i,:]
#             ind = Int(ceil((b[1]-mintime)/step_length))
#             if ind > 0
#                 betts[ind,dim]=b[2]
#             else
#                 betts[1,dim]=b[2]
#             end
#         end
#     end
#     return betts
# end

# average betti numbers over arrs
# assuming arrs is an array of arrays, where each arrs[j] is the same size

function average_bettis(arrs; maxdim=-1)
    if size(arrs,2) > 1
        return arrs
    end
    md = maxdim
    if maxdim == -1
        md = size(arrs[1],2)
    end
    numofints = size(arrs[1],1)
    av_bet = zeros(numofints,md)

    for i=1:numofints
        for d=1:md
            av_bet[i,d] = mean([arrs[j][i,d] for j=1:length(arrs)])
        end
    end
    return av_bet
end

# compute standard deviation of betti numbers over arrays in arrs
# assuming arrs is an array of arrays, where each arrs[j] is the same size

function std_bettis(arrs; maxdim=-1)
    md = maxdim
    if maxdim == -1
        md = size(arrs[1],2)
    end
    numofints = size(arrs[1],1)
    std_bet = zeros(numofints,md)

    if size(arrs,2) > 1
        return std_bet
    end

    for i=1:numofints
        for d=1:md
            std_bet[i,d] = std([arrs[j][i,d] for j=1:length(arrs)])
        end
    end
    return std_bet
end

# plot average curves at values `xval`, with averages given by `means` and standard deviations given by `std`

function plot_averages(xvals, means, stds; ribbon=true, label="", linestyle=:solid, color=:auto)
    if ribbon
         return plot(xvals, means, ribbon=stds,fillalpha=.3, labels=label, linestyle=linestyle, color=color)
    else
        return plot(xvals, means, labels=label, linestyle=linestyle, c=color)
    end
end

function plot_averages!(xvals, means, stds; ribbon=true, label="", linestyle=:solid, color=:auto)
    if ribbon
         return plot!(xvals, means, ribbon=stds,fillalpha=.3, labels=label, linestyle=linestyle, color=color)
    else
        return plot!(xvals, means, labels=label, linestyle=linestyle, c=color)
    end
end

function load_bettis(filename)
    dict = load(filename)
    for (matr_name, matr) in dict
        return matr
    end
end

# plot average curves at values `xval`, given that the bettis numbers are saved in `file`

function plot_averages(xvals, file::String; dim=1, ribbon=true, label="", linestyle=:solid, color=:auto)
    matr = load_bettis(file)

    av = average_bettis(matr)[:,dim]
    if ribbon
         st = std_bettis(matr)[:,dim]
         return plot(xvals, av, ribbon=st,fillalpha=.3, labels=label, linestyle=linestyle, c=color)
    else
        return plot(xvals, av, labels=label, linestyle=linestyle, c=color)
    end
end

function plot_averages!(xvals, file::String; dim=1, ribbon=true, label="", linestyle=:solid, color=:auto)
    matr = load_bettis(file)

    av = average_bettis(matr)[:,dim]
    if ribbon
         st = std_bettis(matr)[:,dim]
         return plot!(xvals, av, ribbon=st,fillalpha=.3, labels=label, linestyle=linestyle, c=color)
    else
        return plot!(xvals, av, labels=label, linestyle=linestyle, c=color)
    end
end


text:
using MeasureSystems, Test

@test 1000molarmass(SI2019) == molarmass(EMU2019)
@test luminousefficacy(SI2019)/1e7 ≈ luminousefficacy(EMU2019)


text:
using ChainRulesCore
using Documenter
using DocumenterTools: Themes
using Markdown

DocMeta.setdocmeta!(
    ChainRulesCore,
    :DocTestSetup,
    quote
        using Random
        Random.seed!(0)  # frule doctest shows output

        using ChainRulesCore
        # These rules are all actually defined in ChainRules.jl, but we redefine them here to
        # avoid the dependency.
        @scalar_rule(sin(x), cos(x))  # frule and rrule doctest
        @scalar_rule(sincos(x), @setup((sinx, cosx) = Ω), cosx, -sinx)  # frule doctest
        @scalar_rule(hypot(x::Real, y::Real), (x / Ω, y / Ω))  # rrule doctest
    end
)

Themes.compile(joinpath(@__DIR__, "src/assets/chainrules.scss"))

makedocs(
    modules=[ChainRulesCore],
    format=Documenter.HTML(
        prettyurls=false,
        assets=["assets/chainrules.css"],
        mathengine=MathJax(),
    ),
    sitename="ChainRules",
    authors="Jarrett Revels and other contributors",
    pages=[
        "Introduction" => "index.md",
        "FAQ" => "FAQ.md",
        "Writing Good Rules" => "writing_good_rules.md",
        "Complex Numbers" => "complex.md",
        "Deriving Array Rules" => "arrays.md",
        "Debug Mode" => "debug_mode.md",
        "Gradient Accumulation" => "gradient_accumulation.md",
        "Usage in AD" => [
            "Overview" => "autodiff/overview.md",
            "Operator Overloading" => "autodiff/operator_overloading.md",
        ],
        "Design" => [
            "Changing the Primal" => "design/changing_the_primal.md",
            "Many Differential Types" => "design/many_differentials.md",
        ],
        "API" => "api.md",
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(
    repo = "github.com/JuliaDiff/ChainRulesCore.jl.git",
    push_preview=true,
)


text:
# TODO there might be Julia library for handling such configurations
# TODO maybe add possibility to add config values from Environment
# TODO add way to validate value and type of configs

config_values = Dict{Symbol, Any}()

supported_config_keys = [:scale, :transfer_mode, :prefer_png_to_svg]

function set_kitty_config!(key, value)

    if key ∉ supported_config_keys
        throw(DomainError(key, "$key is not a valid config key. Valid keys are $supported_config_keys"))
    end

    # TODO this could be handled more general
    if key == :transfer_mode
        if value ∉ (:direct, :temp_file)
            throw(DomainError(value, "value for key :transfer_mode must be either :direct or :temp_file"))
        end
    end
    if key == :prefer_png_to_svg
        if value ∉ (true, false)
            throw(DomainError(value, "value for key :prefer_png_to_svg must be either true or false"))
        end
    end

    config_values[key] = value
    return nothing
end

get_kitty_config(key::Symbol) = config_values[key]
get_kitty_config(key::Symbol, default) = get(config_values, key, default)

# insert some default
set_kitty_config!(:transfer_mode, :direct)
set_kitty_config!(:prefer_png_to_svg, true)


text:
const DEFAULT_GRADIENT_ϵ = 8sqrt(eps())

abstract type GradientParameters{T} end

compute_gradient!(g::AbstractVector, x::AbstractVector, grad::GradientParameters) = grad(g,x)

function compute_gradient(x, grad::GradientParameters)
    g = alloc_g(x)
    grad(g,x)
    return g
end

function check_gradient(g::AbstractVector)
    println("norm(Gradient):               ", norm(g))
    println("minimum(|Gradient|):          ", minimum(abs.(g)))
    println("maximum(|Gradient|):          ", maximum(abs.(g)))
    println()
end

function print_gradient(g::AbstractVector)
    display(g)
    println()
end

struct GradientParametersUser{T, ∇T} <: GradientParameters{T}
    ∇F!::∇T
end

function (grad::GradientParametersUser{T})(g::AbstractVector{T}, x::AbstractVector{T}) where {T}
    grad.∇F!(g, x)
end

struct GradientParametersAD{T, FT, ∇T <: ForwardDiff.GradientConfig} <: GradientParameters{T}
    F::FT
    ∇config::∇T

    function GradientParametersAD(F::FT, x::VT) where {T, FT, VT <: AbstractVector{T}}
        ∇config = ForwardDiff.GradientConfig(F, x)
        new{T, FT, typeof(∇config)}(F, ∇config)
    end
end

function GradientParametersAD{T}(F::FT, nx::Int) where {T, FT}
    GradientParametersAD(F, zeros(T, nx))
end

function (grad::GradientParametersAD{T})(g::AbstractVector{T}, x::AbstractVector{T}) where {T}
    ForwardDiff.gradient!(g, grad.F, x, grad.∇config)
end

function compute_gradient_ad!(g::AbstractVector{T}, x::AbstractVector{T}, F) where {T}
    grad = GradientParametersAD{T}(F, length(x))
    grad(g,x)
end

struct GradientParametersFD{T, FT} <: GradientParameters{T}
    ϵ::T
    F::FT
    e::Vector{T}
    tx::Vector{T}
end

function GradientParametersFD{T}(F::FT, nx::Int; ϵ=DEFAULT_GRADIENT_ϵ) where {T, FT}
    e  = zeros(T, nx)
    tx = zeros(T, nx)
    GradientParametersFD{T,FT}(ϵ, F, e, tx)
end

function (grad::GradientParametersFD{T})(g::AbstractVector{T}, x::AbstractVector{T}) where {T}
    local ϵⱼ::T

    for j in eachindex(x,g)
        ϵⱼ = grad.ϵ * x[j] + grad.ϵ
        fill!(grad.e, 0)
        grad.e[j] = 1
        grad.tx .= x .- ϵⱼ .* grad.e
        f1 = grad.F(grad.tx)
        grad.tx .= x .+ ϵⱼ .* grad.e
        f2 = grad.F(grad.tx)
        g[j] = (f2 - f1)/(2ϵⱼ)
    end
end

function compute_gradient_fd!(g::AbstractVector{T}, x::AbstractVector{T}, F::FT; kwargs...) where {T, FT}
    grad = GradientParametersFD{T}(F, length(x); kwargs...)
    grad(g,x)
end

function GradientParameters{T}(ForG, nx::Int; mode = :autodiff, diff_type = :forward, kwargs...) where {T}
    if mode == :autodiff
        if diff_type == :forward
            Gparams = GradientParametersAD{T}(ForG, nx)
        else
            Gparams = GradientParametersFD{T}(ForG, nx; kwargs...)
        end
    else
        Gparams = GradientParametersUser{T, typeof(ForG)}(ForG)
    end
    return Gparams
end

GradientParameters{T}(∇F!, F, nx::Int; kwargs...) where {T} = GradientParameters{T}(∇F!, nx; mode = :user, kwargs...)

GradientParameters{T}(∇F!::Nothing, F, nx::Int; kwargs...) where {T} = GradientParameters{T}(F, nx;  mode = :autodiff, kwargs...)

GradientParameters(∇F!, F, x::AbstractVector{T}; kwargs...) where {T} = GradientParameters{T}(∇F!, F, length(x); kwargs...)

GradientParameters(F, x::AbstractVector{T}; kwargs...) where {T} = GradientParameters{T}(F, length(x); kwargs...)

function compute_gradient!(g::AbstractVector{T}, x::AbstractVector{T}, ForG; kwargs...) where {T}
    grad = GradientParameters{T}(ForG, length(x); kwargs...)
    grad(g,x)
end


text:
n = 7 # rand(1:7)
@test matrixdepot("sampling", n) == matrixdepot("sampling", Float64, n)

A = matrixdepot("sampling", n)
v = sort(eigvals(A))
@test v ≈ [0: n - 1;] atol=1e-10
println("'sampling' passed test...")


text:
# module SQLite_Octo.SQL

type SQLiteAdapter <: QueryBuilders.AbstractAdapter
    db::Union{Void, SQLite.DB}
end