text:
# FileID
getHashValue(id::FileID) = clang_FileID_getHashValue(id)

dispose(x::FileID) = clang_FileID_dispose(x)

# SourceLocation
SourceLocation() = SourceLocation(clang_SourceLocation_createInvalid())

isFileID(x::SourceLocation) = clang_SourceLocation_isFileID(x)
isMacroID(x::SourceLocation) = clang_SourceLocation_isMacroID(x)
isValid(x::SourceLocation) = clang_SourceLocation_isValid(x)
isInvalid(x::SourceLocation) = clang_SourceLocation_isInvalid(x)
getHashValue(x::SourceLocation) = clang_SourceLocation_getHashValue(x)

function isPairOfFileLocations(b::SourceLocation, e::SourceLocation)
    return clang_SourceLocation_isPairOfFileLocations(b, e)
end

function getLocWithOffset(x::SourceLocation, offset::Integer)
    return SourceLocation(clang_SourceLocation_getLocWithOffset(x, offset))
end

function printToString(x::SourceLocation, src_mgr::SourceManager)
    @check_ptrs src_mgr
    str = clang_SourceLocation_printToString(x, src_mgr)
    s = unsafe_string(str)
    clang_SourceLocation_disposeString(str)
    return s
end

# SourceRange
getBeginLoc(x::SourceRange) = x.begin_loc
getEndLoc(x::SourceRange) = x.end_loc


text:
struct CPUDevice <: Device end

export CPUDevice

Computation.name(device::CPUDevice) = :CPU
Computation.alloc(device::CPUDevice, array::AbstractArray) = Array(array)
Computation.arraytype(device::CPUDevice) = Array
Computation.alloc(device::CPUDevice, T::Type, indxs...) = Array{T}(undef, indxs)

ext: jl


text:
filePath = @__DIR__
(@isdefined(SphericalHarmonics_Mod)) ? nothing : include(string(filePath,"\\SphericalHarmonics.jl"))
using .SphericalHarmonics_Mod
using Plots

savepath = string(filePath,"\\anim.gif")

function main(savepath)

    pyplot(size=(1080,1080))

    anim = @animate for l in 0:20, m = 0:l

        x,y,z,Y = Calc_SphericalHarmonics(l,m,500)
        if l==0
            surface(x,y,z, legend = nothing, show = false, axis=([], false))
        else
            surface(x,y,z, legend = nothing, fill_z = abs.(Y), show = false, axis=([], false))
        end
        title!(string("l = ", l, ", m = ", m))

    end

    gif(anim,savepath,fps=2)

    return

end

main(savepath)

max_stars_repo_name: philipbergander/SphericalHarmonics

max_issues_repo_name: philipbergander/SphericalHarmonics

max_forks_repo_name: philipbergander/SphericalHarmonics


text:
export FourierTransform

struct FourierTransform{N, S}<:AbstractTransform
    modes::NTuple{N, S} # N == ndims(x)
end

Base.ndims(::FourierTransform{N}) where {N} = N

function transform(ft::FourierTransform, 𝐱::AbstractArray)
    return fft(Zygote.hook(real, 𝐱), 1:ndims(ft)) # [size(x)..., in_chs, batch]
end

function low_pass(ft::FourierTransform, 𝐱_fft::AbstractArray)
    return view(𝐱_fft, map(d->1:d, ft.modes)..., :, :) # [ft.modes..., in_chs, batch]
end

const truncate_modes = low_pass

function inverse(ft::FourierTransform, 𝐱_fft::AbstractArray)
    return real(ifft(𝐱_fft, 1:ndims(ft))) # [size(x_fft)..., out_chs, batch]
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using Test

function _timed_include(file)
    println("Testing ", file)
    start = time()
    include(file)
    run_time = round(time() - start, digits = 1)
    return println("    Took $(run_time) seconds")
end

@testset "$(dir)" for dir in [".", "Variable", "Constraint", "Objective"]
    @testset "$(file)" for file in readdir(joinpath(@__DIR__, dir))
        if !endswith(file, ".jl")
            continue
        elseif (dir == ".") && (file == "Bridges.jl" || file == "utilities.jl")
            continue
        end
        _timed_include(joinpath(@__DIR__, dir, file))
    end
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
abstract type GeometricIntegratorAlgorithm <: DiffEqBase.DEAlgorithm end
struct GIEuler <: GeometricIntegratorAlgorithm end
struct GIMidpoint <: GeometricIntegratorAlgorithm end
struct GIHeun2 <: GeometricIntegratorAlgorithm end
struct GIHeun3 <: GeometricIntegratorAlgorithm end
struct GIRalston2 <: GeometricIntegratorAlgorithm end
struct GIRalston3 <: GeometricIntegratorAlgorithm end
struct GIRunge <: GeometricIntegratorAlgorithm end
struct GIKutta <: GeometricIntegratorAlgorithm end
struct GIRK4 <: GeometricIntegratorAlgorithm end
struct GIRK416 <: GeometricIntegratorAlgorithm end
struct GIRK438 <: GeometricIntegratorAlgorithm end
struct GISSPRK3 <: GeometricIntegratorAlgorithm end
struct GICrankNicolson <: GeometricIntegratorAlgorithm end
struct GIKraaijevangerSpijker <: GeometricIntegratorAlgorithm end
struct GIQinZhang <: GeometricIntegratorAlgorithm end
struct GICrouzeix <: GeometricIntegratorAlgorithm end
struct GIImplicitEuler <: GeometricIntegratorAlgorithm end
struct GIImplicitMidpoint <: GeometricIntegratorAlgorithm end
struct GISRK3 <: GeometricIntegratorAlgorithm end
struct GIGLRK <: GeometricIntegratorAlgorithm
    s::Int
end
struct GIRadauIA <: GeometricIntegratorAlgorithm
    s::Int
end
struct GIRadauIIA <: GeometricIntegratorAlgorithm
    s::Int
end
struct GILobattoIIIA <: GeometricIntegratorAlgorithm
    s::Int
end
struct GILobattoIIIB <: GeometricIntegratorAlgorithm
    s::Int
end
struct GILobattoIIIC <: GeometricIntegratorAlgorithm
    s::Int
end
struct GILobattoIIIC̄ <: GeometricIntegratorAlgorithm
    s::Int
end
struct GILobattoIIID <: GeometricIntegratorAlgorithm
    s::Int
end
struct GILobattoIIIE <: GeometricIntegratorAlgorithm
    s::Int
end
struct GILobattoIIIF <: GeometricIntegratorAlgorithm
    s::Int
end
struct GISymplecticEulerA <: GeometricIntegratorAlgorithm end
struct GISymplecticEulerB <: GeometricIntegratorAlgorithm end
struct GILobattoIIIAIIIB2 <: GeometricIntegratorAlgorithm end
struct GILobattoIIIBIIIA2 <: GeometricIntegratorAlgorithm end


text:
using FractionalDiffEq
using Plots

function Qi!(du, u, p, t)
    a, b, c, d, r = 35, 8/3, 80, -1, 1
    du[1] = -a*u[1]+a*u[2]+r*u[2]*u[3]
    du[2] = c*u[1]+d*u[2]-u[1]*u[3]
    du[3] = -b*u[3]+u[1]*u[2]
end

alpha = [0.98, 0.98, 0.98]
h=0.01
T=50
x0=[0.1, 0.2, 0.3]
prob=FODESystem(Qi!, alpha, x0)
result = solve(prob, h, T, GLWithMemory())
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Qi System")


text:
function Base.real(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    n=ncoefficients(f)
    cfs=f.coefficients

    ret=Array{real(cfstype(f))}(undef, iseven(n) ? n+1 : n)
    ret[1]=real(cfs[1])

    for k=2:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]=imag(cfs[k])
        ret[k+1]=real(cfs[k])
    end
    for k=3:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]+=real(cfs[k])
        ret[k-1]-=imag(cfs[k])
    end

    Fun(Fourier(domain(f)),ret)
end

function Base.imag(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    n=ncoefficients(f)
    cfs=f.coefficients

    ret=Array{real(cfstype(f))}(undef, iseven(n) ? n+1 : n)
    ret[1]=imag(cfs[1])

    for k=2:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]=-real(cfs[k])
        ret[k+1]=imag(cfs[k])
    end
    for k=3:2:n
        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
        ret[k]+=imag(cfs[k])
        ret[k-1]+=real(cfs[k])
    end

    Fun(Fourier(domain(f)),ret)
end


text:
import ChainRulesCore: rrule, no_rrule
import ChainRulesCore: rrule_via_ad, RuleConfig, NoForwardsMode, HasReverseMode
import Umlaut: make_name, Input, to_expr, BcastCtx

###############################################################################
#                              Primitives                                     #
###############################################################################

struct ChainRulesCtx end

function isprimitive(::ChainRulesCtx, f, args...)
    F = Core.Typeof(f)
    Args = Core.Typeof.(args)
    Core.Compiler.return_type(rrule, Tuple{YotaRuleConfig, F, Args...}) !== Nothing && return true
    if is_kwfunc(F)
        Args_kwrrule = Tuple{Any, typeof(Core.kwfunc(f)), YotaRuleConfig, Args...,}
        Core.Compiler.return_type(Core.kwfunc(rrule), Args_kwrrule) !== Nothing && return true
    end
    return false
end

###############################################################################
#                              RuleConfig                                     #
###############################################################################

"""
    YotaRuleConfig()

ChainRules.RuleConfig passed to all `rrule`s in Yota.
Extends RuleConfig{Union{NoForwardsMode,HasReverseMode}}.
"""
struct YotaRuleConfig <: RuleConfig{Union{NoForwardsMode,HasReverseMode}} end

###############################################################################
#                              rrule_via_ad                                   #
###############################################################################

"""
    bcast_rrule(::YotaRuleConfig, ::typeof(broadcasted), f, args...; kw...)

Similar to rrule(config, broadcasted, f, args...), but works on for ChainRule-primitive
functions. For a more flexible handling of broadcasting use rrule(...) directly.
"""
function bcast_rrule(::YotaRuleConfig, ::typeof(broadcasted), f::F, args...; kw...) where F
    ys, pbs = unzip(rrule.(YOTA_RULE_CONFIG, f, args...; kw...))
    function pullback(Δ)
        if Δ isa NoTangent || Δ isa ZeroTangent
            return (NoTangent(), [Δ for _=1:length(pbs) + 1]...,)
        end
        Δ = unthunk(Δ)
        dxs = map((pb, Δ) -> pb(Δ), pbs, Δ) |> unzip
        dxs = [all(dx .== NoTangent()) ? NoTangent() : dx for dx in dxs]
        return (NoTangent(), dxs...,)
    end
    return ys, pullback
end

function to_rrule_expr(tape::Tape)
    # TODO (maybe): add YotaRuleConfig() as the first argument for consistency
    fn_name = gensym("rrule_$(tape[V(1)].val)")
    header = Expr(:call, fn_name)
    push!(header.args, Expr(:(::), :config, YotaRuleConfig))
    for v in inputs(tape)
        op = tape[v]
        push!(header.args, Expr(:(::), make_name(op), op.typ))
    end
    body = Expr(:block)
    # generate transformed forward pass
    seed_id = tape.meta[:seed].id
    for op in tape.ops[1:seed_id - 1]
        op isa Input && continue
        ex = to_expr(op)
        if ex isa Vector
            push!(body.args, ex...)
        else
            push!(body.args, ex)
        end
    end
    # generate pullback
    pb_name = gensym("pullback_$(tape[V(1)].val)")
    pb_ex = :(function $pb_name(dy) end)
    pb_body = pb_ex.args[2]
    empty!(pb_body.args)  # clean from useless linenumber nodes
    push!(pb_body.args, Expr(:(=), make_name(tape.meta[:seed].id), :dy))
    for op in tape.ops[seed_id + 1:length(tape) - 2]
        op isa Input && continue
        ex = to_expr(op)
        if ex isa Vector
            push!(pb_body.args, ex...)
        else
            push!(pb_body.args, ex)
        end
    end
    push!(body.args, pb_ex)
    # generate return
    result_name = make_name(tape[tape.result].args[1].id)
    push!(body.args, Expr(:tuple, result_name, pb_name))
    fn_ex = Expr(:function, header, body)
    return fn_ex
end

"""
    make_rrule(tape::Tape)
    make_rrule(f, args...)

Generate a function equivalent to (but not extending) ChainRulesCore.rrule(),
i.e. returning the primal value and the pullback.

Examples:
=========

    foo(x) = 2x + 1
    rr = make_rrule(foo, 2.0)
    val, pb = rr(foo, 3.0)
    pb(1.0)

"""
make_rrule(tape::Tape) = Base.eval(@__MODULE__, to_rrule_expr(tape))

function make_rrule(f, args...)
    return make_rrule(gradtape(f, args...; seed=:auto, ctx=GradCtx()))
end

function make_rrule(::typeof(broadcasted), f, args...)
    if isprimitive(GradCtx(), f, map(first, args)...)
        return bcast_rrule # (YOTA_RULE_CONFIG, broadcasted, f, args...)
    end
    ctx = BcastGradCtx(GradCtx())
    _, tape = trace(f, args...; ctx=ctx, fargtypes=(f, map(eltype, args)))
    tape = Tape(tape; ctx=ctx.inner)
    gradtape!(tape, seed=:auto)
    # insert imaginary broadcasted to the list of inputs
    insert!(tape, 1, Umlaut.Input(broadcasted))
    # insert ZeroTangent to the result to account for the additional argument
    grad_tuple_op = tape[V(tape.result.id - 2)]
    @assert grad_tuple_op isa Call && grad_tuple_op.fn == tuple
    grad_tuple_op.args = [ZeroTangent(), grad_tuple_op.args...]
    for id=grad_tuple_op.id:grad_tuple_op.id + 2
        Umlaut.exec!(tape, tape[V(id)])
    end
    return make_rrule(tape)
end

const GENERATED_RRULE_CACHE = Dict()

"""
    rrule_via_ad(::YotaRuleConfig, f, args...)

Generate `rrule` using Yota.
"""
function ChainRulesCore.rrule_via_ad(::YotaRuleConfig, f, args...)
    res = rrule(f, args...)
    !isnothing(res) && return res
    sig = map(typeof, (f, args...))
    if haskey(GENERATED_RRULE_CACHE, sig)
        rr = GENERATED_RRULE_CACHE[sig]
        # return Base.invokelatest(rr, f, args...)
        val, pb = Base.invokelatest(rr, YOTA_RULE_CONFIG, f, args...)
        return val, dy -> Base.invokelatest(pb, dy)
    else
        rr = make_rrule(f, args...)
        GENERATED_RRULE_CACHE[sig] = rr
        # return Base.invokelatest(rr, f, args...)
        val, pb = Base.invokelatest(rr, YOTA_RULE_CONFIG, f, args...)
        return val, dy -> Base.invokelatest(pb, dy)
    end
end


text:
using MultidimensionalTools
using Test

function aresamearrays(A::AbstractVector{T}, B::AbstractVector{R}) where {T, R}
	A === B && return true
	length(A) ≠ length(B) && return false

	for (a, b) in zip(sort(A), sort(B))
		if a ≠ b
			return false
		end
	end

	return true
end
aresamearrays(A::AbstractArray{T, N}, B::AbstractArray{R, M}) where {T, R, N, M} =
	aresamearrays(reshape(A, :), reshape(B, :))

# TESTS

@time @testset "MultidimensionalTools.jl Types" begin
	# TYPES
	@test (0, 0, 0) isa AbstractIndex{3}
	@test CartesianIndex(1, 2, 3, 4, 5) isa AbstractIndex{5}
	@test [(0, 0), (1, 1), (2, 2)] isa AbstractIndices{2}
	@test CartesianIndex.([(3, 6, 9, 12), (4, 4, 1, 0)]) isa AbstractIndices{4}
	@test CartesianIndices([1, 2, 3, 4]) isa AbstractIndices{1}
	@test CartesianIndices(rand(Int8, 3, 3, 3, 3, 3, 3, 3)) isa AbstractIndices{7}
	@test (0, 0, 0) isa AbstractIndexOrIndices{3}
	@test CartesianIndex(1, 2, 3, 4, 5) isa AbstractIndexOrIndices{5}
	@test [(0, 0), (1, 1), (2, 2)] isa AbstractIndexOrIndices{2}
	@test CartesianIndex.([(3, 6, 9), (4, 4, 1)]) isa AbstractIndexOrIndices{3}
	@test CartesianIndices([1, 2, 3, 4]) isa AbstractIndexOrIndices{1}
	@test CartesianIndices(rand(Int8, 3, 3, 3, 3, 3, 3, 3)) isa AbstractIndexOrIndices{7}
end

@time @testset "MultidimensionalTools.jl Main" begin
	# MAIN
    @test n_adjacencies(2) == 8
    @test n_adjacencies(3) == 26
    @test n_adjacencies(17) == 129140162
	@test n_adjacencies(cat([1 1; 1 1], [1 1; 1 1], [1 1; 1 1], dims = 3)) == n_adjacencies(3)
    @test aresamearrays(get_directions(3), [(-1, -1, -1), (0, -1, -1), (1, -1, -1), (-1, 0, -1), (0, 0, -1), (1, 0, -1), (-1, 1, -1), (0, 1, -1), (1, 1, -1), (-1, -1, 0), (0, -1, 0), (1, -1, 0), (-1, 0, 0), (1, 0, 0), (-1, 1, 0), (0, 1, 0), (1, 1, 0), (-1, -1, 1), (0, -1, 1), (1, -1, 1), (-1, 0, 1), (0, 0, 1), (1, 0, 1), (-1, 1, 1), (0, 1, 1), (1, 1, 1)])
	A = Int[-69 -63 -5; 119 67 1; -101 7 -88]
	@test tryindex(A, (1, 1), (2, 2), (3, 3), (4, 4), (5, 5)) == (-69, 67, -88, nothing, nothing)
	@test tryindex(A, (1, 1)) == (-69, )
	# @test tryindex(A, CartesianIndex((1, 1))) == (-69, )
	T = [(3, 1), (4, 5), (6, 2), (1, 1)]
	@test extrema_indices(T) == ((1, 6), (1, 5))
	@test extrema_indices(T...) == ((1, 6), (1, 5))
	@test append_n_times(A, 5, 1, dims = 1) == [-69 -63 -5; 119 67 1; -101 7 -88; 1 1 1; 1 1 1; 1 1 1; 1 1 1; 1 1 1]
	@test append_n_times(A, 5, 1, dims = 2) == [-69 -63 -5 1 1 1 1 1; 119 67 1 1 1 1 1 1; -101 7 -88 1 1 1 1 1]
	@test append_n_times(A, 2, 1, dims = 3) == cat([-69 -63 -5; 119 67 1; -101 7 -88], [1 1 1; 1 1 1; 1 1 1], [1 1 1; 1 1 1; 1 1 1], dims = 3)
	@test append_n_times_backwards(A, 5, 1, dims = 1) == [1 1 1; 1 1 1; 1 1 1; 1 1 1; 1 1 1; -69 -63 -5; 119 67 1; -101 7 -88]
	@test append_n_times_backwards(A, 5, 1, dims = 2) == [1 1 1 1 1 -69 -63 -5; 1 1 1 1 1 119 67 1; 1 1 1 1 1 -101 7 -88]
	@test promote_to_nD(A, 3, 0) == cat([0 0 0; 0 0 0; 0 0 0], [-69 -63 -5; 119 67 1; -101 7 -88], [0 0 0; 0 0 0; 0 0 0], dims = 3)
	@test promote_to_3D(A, 0) == cat([0 0 0; 0 0 0; 0 0 0], [-69 -63 -5; 119 67 1; -101 7 -88], [0 0 0; 0 0 0; 0 0 0], dims = 3)
	# @test expand_as_required(A, 0, CartesianIndex((4, 4))) == [-69 -63 -5 0; 119 67 1 0; -101 7 -88 0; 0 0 0 0]
	@test expand_as_required(A, 0, [(1, 1), (2, 2), (3, 3), (4, 4), (5, 5)]) == [-69 -63 -5 0 0; 119 67 1 0 0; -101 7 -88 0 0; 0 0 0 0 0; 0 0 0 0 0]
	# @test expand_as_required(A, 0, [(1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (-2, -2)]) == [-69 -63 -5 0 0; 119 67 1 0 0; -101 7 -88 0 0; 0 0 0 0 0; 0 0 0 0 0]
	@test aresamearrays(adjacencies(A, (2, 2)), [-69, 119, -101, -63, 7, -5, 1, -88])
	@test aresamearrays(adjacencies(A, (3, 3)), [67, 7, 1])
	@test aresamearrays(expanded_adjacencies(A, 0, (2, 2)), [-69, 119, -101, -63, 7, -5, 1, -88])
	@test aresamearrays(expanded_adjacencies(A, 0, (3, 3)), [67, 7, 0, 1, 0, 0, 0, 0])
	@test n_adjacent_to(A, (2, 2), 7) == 1
	@test n_adjacent_to(A, (3, 3), 7) == 1
	@test n_adjacent_to(promote_to_3D(A, 7), (3, 3, 1), 7) == 4
	@test expanded_n_adjacent_to(A, 7, (3, 3), 7) == 6
	@test expanded_n_adjacent_to(A, 7, (6, 6), 7) == 8 # (the filled element is surrounding the given index)
	@test expanded_n_adjacent_to(promote_to_3D(A, 7), 7, (3, 3, 1), 7) == 24
	@test aresamearrays(global_adjacencies_indices(A, (2, 2), 7), [(1, 1), (2, 1), (3, 1), (1, 2), (1, 3), (2, 3), (3, 3)])
	@test aresamearrays(global_adjacencies_indices(promote_to_3D(A, 7), (3, 3, 1), 7), [(2, 2, 2), (2, 3, 2), (3, 3, 2)])
	@test aresamearrays(global_adjacencies(A, (3, 3), 67), [-69, 7, 1])
	@test global_n_adjacent_to(A, (3, 3), -69, 67) == 1
    @test aresamearrays(reduce_dims(zeros(Int, 3, 3, 3)), fill(zeros(Int, 3), 3, 3))
    @test aresamearrays(reduce_dims(zeros(Int, 20, 3, 10, 100)), fill(zeros(Int, 100), 20, 3, 10))
    @test aresamearrays(reduce_dims(zeros(Int, 3, 3)), fill(zeros(Int, 3), 3))
end # end testset


text:
module HawkesProcesses

    using Distributions
    using StatsBase

    include("sample_parents.jl")
    include("event_difference.jl")
    include("parent_probabilites.jl")
    include("all_constant_hawkes.jl")
    include("hierarchical_all_constant.jl")
    include("simulate.jl")
    include("intensity.jl")
    include("likelihood.jl")
    include("time_change.jl")
    include("compensator.jl")
    include("utils.jl")
    include("event_probability.jl")
    include("simulate_forward.jl")

    export fit, hierarchical_fit,
    intensity, compensator,
    likelihood,
    simulate, simulate_forward

end # module


text:
module EMLcolumn_cost
	using DanaTypes
	using DotPlusInheritance
	using Reexport
end


text:
astar(args...; do_action = :stroke) = star(Point(-100, -100), 30, 5, 0.5, 0, do_action)
acirc(args...; do_action = :stroke) = circle(Point(100, 100), 30, do_action)

@testset "morphing star2circle and back" begin
    video = Video(500, 500)

    back = Background(1:20, (args...) -> ground_color("white", "black", args[3]))
    Object(1:10, (args...) -> circle(Point(-100, 0), val(back), :fill))
    star_obj = Object(1:10, astar)
    act!(star_obj, Action(linear(), morph_to(acirc)))
    circle_obj = Object(11:20, acirc)
    act!(circle_obj, Action(:same, morph_to(astar)))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/star2circle5.png" load("images/0000000005.png")
    @test_reference "refs/star2circle15.png" load("images/0000000015.png")
    for i in 1:20
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "morphing star2circle and back with fill" begin
    video = Video(500, 500)
    back = Background(1:20, (args...) -> ground_color("white", "black", args[3]))
    Object(1:10, (args...) -> circle(Point(-100, 0), val(back), :fill))
    star_obj = Object(1:10, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    circle_obj = Object(11:20, acirc)
    act!(circle_obj, Action(morph_to(astar; do_action = :fill)))

    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/star2circle_fill5.png" load("images/0000000005.png")
    @test_reference "refs/star2circle_fill15.png" load("images/0000000015.png")
    for i in 1:20
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "test default kwargs" begin
    video = Video(500, 500)
    back = Background(1:10, (args...) -> ground_color("white", "black", args[3]))
    Object(1:10, (args...) -> circle(Point(-100, 0), val(back), :fill))
    star_obj = Object(1:10, astar)
    act!(star_obj, Action(morph_to(acirc)))

    pathname = render(video)

    path, ext = splitext(pathname)
    @test ext == ".gif"
    @test isfile(pathname)
    rm(pathname)
end

@testset "transposing a matrix" begin
    function latex_ground(args...)
        translate(-200, -100)
        background("white")
        sethue("black")
        fontsize(50)
    end

    function matrix(args...; do_transpose = false, do_action = :stroke)
        fontsize(50)
        str =
            L"$\begin{equation}\left[\begin{array}{ccc}\alpha & \beta & \gamma \\x^{2} & \sqrt{y} & \lambda \\1 & 2 & y \\\end{array}\right]\end{equation}$"
        if do_transpose
            str =
                L"$\begin{equation}\left[\begin{array}{ccc}\alpha & x^{2} & 1 \\\beta & \sqrt{y} & 2 \\\gamma & \lambda & y \\\end{array}\right]\end{equation}$"
        end
        latex(str, O, do_action)
    end

    video = Video(600, 400)
    Background(1:62, latex_ground)
    m = Object(1:60, matrix)
    act!(
        m,
        Action(
            31:60,
            morph_to(
                (args...; do_action = :fill) ->
                    matrix(; do_transpose = true, do_action = do_action);
                do_action = :fill,
            ),
        ),
    )
    Object(61:62, (args...) -> matrix(; do_transpose = true))
    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/matrix_transpose1.png" load("images/0000000001.png")
    @test_reference "refs/matrix_transpose50.png" load("images/0000000050.png")
    @test_reference "refs/matrix_transpose55.png" load("images/0000000055.png")
    @test_reference "refs/matrix_transpose62.png" load("images/0000000062.png")
    for i in 1:62
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Matrix sequence" begin
    function sequence(args...; second = false, do_action = :stroke)
        fontsize(50)

        str =
            L"$\begin{equation}\left[\begin{array}{cc}2 & -1\\-1 & 2 \\\end{array}\right]\end{equation}$"
        if second
            str =
                L"$\begin{equation}\left[\begin{array}{ccc} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \\\end{array}\right]\end{equation}$"
        end
        if !second
            latex(str, Point(-150, -60), do_action)
        else
            latex(str, Point(-200, -120), do_action)
        end
    end

    video = Video(600, 400)
    Background(1:120, ground)
    seq = Object(1:60, sequence)
    act!(
        seq,
        Action(
            31:60,
            morph_to(
                (args...; do_action = :stroke) ->
                    sequence(; second = true, do_action = do_action);
                do_action = :fill,
            ),
        ),
    )
    seq2 = Object(
        61:120,
        (args...; do_action = :stroke) ->
            sequence(; second = true, do_action = do_action),
    )
    act!(
        seq2,
        Action(
            31:60,
            morph_to(
                (args...; do_action = :fill) -> sequence(; do_action = do_action);
                do_action = :fill,
            ),
        ),
    )
    render(video; tempdirectory = "images", pathname = "")

    @test_reference "refs/matrix_sequence01.png" load("images/0000000001.png")
    @test_reference "refs/matrix_sequence50.png" load("images/0000000050.png")
    @test_reference "refs/matrix_sequence80.png" load("images/0000000080.png")
    @test_reference "refs/matrix_sequence99.png" load("images/0000000099.png")
    for i in 1:120
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Morphing mutates the object" begin
    function ground(args...)
        background("black")
        sethue("white")
    end

    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    abox(args...; do_action = :stroke) = rect(-50, -50, 100, 100, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(0, 0), 50, do_action)

    video = Video(500, 500)
    back = Background(1:200, ground)
    star_obj = Object(1:200, abox)
    act!(star_obj, Action(10:20, morph_to(acirc)))
    act!(star_obj, Action(30:40, anim_translate(Point(100, -100))))

    # this is also a bug with the morph_to function
    # the star is formed at the incorrect position
    # the the origin is shifted and neve restored back
    # you can see this in the result of anim_translate actions
    # TODO: fix this at some point!
    act!(star_obj, Action(40:60, morph_to(astar)))
    act!(star_obj, Action(70:90, anim_translate(Point(100, -100))))
    act!(star_obj, Action(100:120, morph_to(abox)))
    act!(star_obj, Action(130:150, anim_translate(Point(-100, 50))))
    act!(star_obj, Action(160:180, morph_to(acirc)))
    render(video; tempdirectory = "images", pathname = "")

    for i in [
        1,
        10,
        20,
        21,
        30,
        31,
        40,
        41,
        60,
        61,
        70,
        91,
        101,
        121,
        126,
        131,
        150,
        151,
        161,
        180,
        190,
    ]
        @test_reference "refs/morph_mutate$i.png" load("images/$(lpad(i, 10, "0")).png")
    end

    for i in 1:200
        rm("images/$(lpad(i, 10, "0")).png")
    end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api;
#
# import typedb.client.api.database.DatabaseManager;
#
# import javax.annotation.CheckReturnValue;
#
# public interface TypeDBClient extends AutoCloseable {
#
#     @CheckReturnValue
#     boolean isOpen();
#
#     @CheckReturnValue
#     DatabaseManager databases();
#
#     @CheckReturnValue
#     TypeDBSession session(String database, TypeDBSession.Type type);
#
#     @CheckReturnValue
#     TypeDBSession session(String database, TypeDBSession.Type type, TypeDBOptions options);
#
#     @CheckReturnValue
#     boolean isCluster();
#
#     @CheckReturnValue
#     TypeDBClient.Cluster asCluster();
#
#     void close();
#
#     interface Cluster extends TypeDBClient {
#
#         @Override
#         @CheckReturnValue
#         DatabaseManager.Cluster databases();
#     }
# }


text:
#= -------------------------------------------------------------------------
# @file problema1_sin_disipador_parametric.jl
#
# @date 11/25/17 00:36:32
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
---------------------------------------------------------------------------=#
using DifferentialEquations
using PyPlot

step_power(t) = 3.871
pulse(t) = (175 * sin(2*π*100*t) + 175) / 2

family = @ode_def heat_simulation begin
   dT_c = (-T_c / C_c)*( 1 / (R_jc) + 1 / (R_ca))  + (T_j / (R_jc * C_c)) + (T_a / (R_ca * C_c))
   dT_j = (T_c/ (C_j * R_jc)) - T_j / (C_j * R_jc) + (step_power(t) / C_j)
end R_jc=>1.0 C_c=>6.8 R_ca=>30.0 C_j=>0.5 T_a=>30.0

T₀ = [0.0 0.0]
tspan = (0.0, 1500.0)
problems = [ODEProblem(heat_simulation(T_a=parameter), T₀, tspan) for parameter in 10:30]
solutions = solve.(problems, Vern7(), abstol=1/10^14, reltol=1/10^14)
for (num, solution) in enumerate(solutions)
   plot(solution.t, solution[2, :], label=latexstring("T_a="*"$(num + 10)"))
   legend(bbox_to_anchor=(1.09,1), loc="upper right", ncol=1)
   xlabel(L"t")
   ylabel(L"T")
end


text:
@testset "Boolean Decision Diagrams" begin
    l1 = pliteral(1)
    @show l1
    @show ¬l1
    @show apply(~,l1)
    @show nliteral(1)
    @show l1 ∧ (¬l1)
    @test value(l1 ∧ (¬l1)) == false
    @show l1 ∨ (¬l1)
    @test value(l1 ∨ (¬l1)) == true
    @info "Fig. 7 in Bryant's paper"
    l2, l3 = pliteral(2), pliteral(3)
    @show ¬(l1 ∧ l3) ∨ (l2 ∧ l3)
    println()
end # testset


text:
using BuildExecutable
import Playground

# Get our current working path
deps_dir = dirname(@__FILE__)

# Setup the build directory
build_dir = joinpath(deps_dir, "usr", "build")
mkpath(build_dir)

build_script = "script.jl"

# Actually build the playground executable
info("Trying to build playground executable in $build_dir ...")
build_executable(
    "playground",
    build_script,
    build_dir,
    "generic"; force=true
)

config_file = joinpath(build_dir, "config.yml")

info("Writing default config to $config_file.")
if ispath(config_file)
    backup_file = joinpath(build_dir, ".config.yml_$(Dates.today()).bak")
    info("Backing up existing config file to $backup_file")
    Playground.copy(config_file, backup_file)
end

open(config_file, "w+") do fstream
    write(fstream, Playground.DEFAULT_CONFIG)
end

@unix_only begin
    Playground.copy(joinpath(deps_dir, "usr", "bin", "INSTALL.sh"), joinpath(build_dir, "INSTALL.sh"))
end

Playground.copy(joinpath(deps_dir, "..", "LICENSE"), joinpath(build_dir, "LICENSE"))
Playground.copy(joinpath(deps_dir, "..", "README.md"), joinpath(build_dir, "README.md"))

# Only install the config and executable if the
# PLAYGROUND_INSTALL env variable has been set.
# This is just cause there isn't a `Pkg.install` or
# `Pkg.build("Pkg", install=true)`
install = haskey(ENV, "PLAYGROUND_INSTALL")

# Store our install paths
install_dir = Playground.config_path()
config_installed = joinpath(install_dir, "config.yml")
playground_installed = joinpath(install_dir, "bin", "playground")
playground_compiled = joinpath(build_dir, "playground")

if install
    # Set up the user level playground directory
    info("Setting up user playground directory...")
    mkpath(install_dir)
    mkpath(joinpath(install_dir, "bin"))

    info("Linking playground config to $config_installed.")

    if ispath(config_installed)
        info("~/.playground/config.yml already exists. Skipping.")
        info("Please see $config_file if you have any problems with your existing config.yml file.")
    else
        Playground.mklink(config_file, config_installed)
    end

    info("Linking playground executable to $playground_installed")

    if ispath(playground_installed)
        backup_file = joinpath(install_dir, "bin", ".playground_$(Dates.today()).bak")
        info("Backing up existing playground executable to $backup_file")
        Playground.copy(playground_installed, backup_file)
    end

    Playground.mklink(playground_compiled, playground_installed)

    info(
        "Adding $(joinpath(install_dir, "bin")) to your PATH " *
        "variable will make `playground` and any julia versions installed via" *
        "Playground.jl available on your search path."
    )
else
    warn(
        "Compiled playground executable $playground_compiled " *
        "not installed to $playground_installed"
    )
end


text:
using Weave
using Test

#Test if running document with and without cache works
isdir("documents/cache") && rm("documents/cache", recursive = true)
weave("documents/chunk_options.noweb", cache=:all)
result =  read("documents/chunk_options.md", String)
rm("documents/chunk_options.md")
weave("documents/chunk_options.noweb", cache=:all)
cached_result =  read("documents/chunk_options.md", String)
@test result == cached_result

# cache = :user
isdir("documents/cache") && rm("documents/cache", recursive = true)
out = "documents/chunk_cache.md"
Weave.weave("documents/chunk_cache.noweb", cache=:user);
result =  read(out, String)
rm(out)
Weave.weave("documents/chunk_cache.noweb", cache=:user);
cached_result =  read(out, String)
@test result == cached_result

# cache = :all
isdir("documents/cache") && rm("documents/cache", recursive = true)
out = "documents/chunk_cache.md"
Weave.weave("documents/chunk_cache.noweb", cache=:all);
result =  read(out, String)
rm(out)
Weave.weave("documents/chunk_cache.noweb", cache=:all);
cached_result =  read(out, String)
@test result == cached_result


text:
# Metida
# Copyright © 2019-2020 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>

__precompile__()

module MetidaFreq

#using StaticArrays
#using AxisArrays
using Tables, CategoricalArrays, Distributions, Roots, StatsBase

import HypothesisTests

import MetidaBase: AbstractData, AbstractIdData, DataSet, Proportion, PrettyTables
import HypothesisTests: ChisqTest, MultinomialLRTest, FisherExactTest
import Base: ht_keyindex, size, show

export contab, diffci, confint

export metaprop, metapropfixed, metaproprandom

export ChisqTest, MultinomialLRTest, FisherExactTest

include("contab.jl")
include("confint.jl")
include("metaprop.jl")
include("hypothesistest.jl")

end


text:
x1 = [3, 9, 18, 30, 91]
y1 = [1.01, 1.204, 1.54, 1.81, 2.12]
c1 = Curve(x1, y1)

x2 = [5, 12, 18, 30, 125, 291]
y2 = [1.01, 1.204, 1.54, 1.81, 2.12, 7.436]
clog = Curve(x2, y2, logx=true, logy=true)
clogy = Curve(x1, y1, logy=true)

c0 = Curve([3], [5.5]) # curve with single point

# equality
c2 = Curve(c1)
c3 = Curve(c1, logx=true)
@test c1 == c2
@test c1 ≈ c2
@test !(c1 == c3)
@test !(c1 ≈ c3)
@test !(c0 == c1)
@test c0 == Curve(3, 5.5)
@test Curve(Tenor("3D"), 7.8) == Curve([3], [7.8])

@test Curve(reverse(x1), reverse(y1)) == c1

# first / last
@test firstpoint(c1, dims=1) == 3
@test firstpoint(c1, dims=2) == 1.01
@test lastpoint(c1, dims=1) == 91
@test lastpoint(c1, dims=2) == 2.12
@test_throws ErrorException firstpoint(c1, dims=3)
@test_throws ErrorException lastpoint(c1, dims=0)
@test firstpoint(c0) == lastpoint(c0) == 3

@test first(c1, 3) == Curve([3, 9, 18], [1.01, 1.204, 1.54])
@test last(c1, 4) == Curve([9, 18, 30, 91], [1.204, 1.54, 1.81, 2.12])
@test first(c1, 10) == c1
@test last(c1, 100) == c1
@test_throws ErrorException first(c1, 0)
@test_throws ErrorException last(c1, -1)
@test first(c0, 10) == last(c0, 10) == c0

@test filter((x) -> x > 10, c1) == Curve([18, 30, 91], [1.54, 1.81, 2.12])
@test filter((x) -> x < 2, c1, axis=:y) == Curve([3, 9, 18, 30], [1.01, 1.204, 1.54, 1.81])
@test filter((x) -> x > 2, c0) == c0
@test_throws ErrorException filter((x) -> x > 10, c0)

# Interpolation
@test interpolate(5.5, c1) ≈ (1.204-1.01)/(9-3)*(5.5-3)+1.01
c_int = interpolate([3, 9, 14, 31, 33], c1)
@test c_int isa Curve && length(c_int) == 5

@test interpolate(25, clog) ≈ exp((log(1.81)-log(1.54))/(log(30)-log(18))*(log(25)-log(18))+log(1.54))
c_log_int = interpolate([3, 9, 14, 31, 33], clog)
@test c_log_int isa Curve && length(c_log_int) == 5
@test interpolate(5.5, clogy) ≈ exp((log(1.204)-log(1.01))/(9-3)*(5.5-3)+log(1.01))
@test interpolate(c1, clog) isa Curve

@test interpolate(10, c0) == 5.5

# Operations with Scalars
@test (c1 + 2).y == y1 .+ 2 && (c1 + 2).x == x1
@test (1/c1).y == 1 ./ y1 && (1/c1).x == x1
c_add = +(c1, 2, logy=true)
@test c_add.y == y1 .+ 2 && c_add.logy == true
# test if interpolation object gets correctly updated
@test interpolate(5.5, 2c1) ≈ 2((1.204-1.01)/(9-3)*(5.5-3)+1.01)
@test exp(c1).y == exp.(c1.y)
@test sin(c1)/cos(c1) ≈ tan(c1)

# Merges
c1d = Curve([1, 3, 3, 7, 9], [2, 4, 4, 8, 10])
@test c1d.x == [1, 3, 7, 9]
@test c1d.y ==  [2, 4, 8, 10]
@test concat(c1, c1d).x == [1, 3, 7, 9, 18, 30, 91]
@test concat(Curve(5.5, 42.1), c1) == Curve([3, 5.5, 9, 18, 30, 91], [1.01, 42.1, 1.204, 1.54, 1.81, 2.12])

# Operations on Curves
c_sum = c1+c1d
@test c_sum.x == sort(unique(vcat(c1.x, c1d.x)))

sumc1 = c1+c1
@test sumc1.y == (2*c1).y # tests correctness of result
@test sumc1 == 2c1 # tests comparison operator in addition
@test log(clogy/c1d).logx == false

@test ismissing(missing*c1)
@test ismissing(c2/missing)

# apply
@test apply((t,r) -> exp(-r/100*t), c1, logy=true).logy == true
@test apply(x -> 2x, c1, axis=:y) ≈ 2c1
res = apply(x -> 3x, c1, axis=:x)
@test res.x == 3c1.x && res.y == c1.y
res = apply((x, y) -> 2x+y, c1, axis=:xy)
@test res.x == c1.x && res.y == c1.y .+ 2.0 .* c1.x
@test_throws ErrorException apply(x -> 2x, c1, axis=:z)

# test non-standard interpolators
c_const = Curve(x1, y1, method=ItpConstant())
@test interpolate(28, c_const) == 1.81

# test extrapolations
@test interpolate(1, c1) == 1.01 # constant extrapolation
c_42 = Curve(x1, y1, extrapolation=42)
@test interpolate(1042, c_42) == 42
c_line = Curve(x1, y1, extrapolation=EtpLine())
@test interpolate(102, c_line) > 0 # just smoke test

# test pretty printing
io = IOBuffer(append=true)
print(io, Curve([t"2D", t"1W", t"3M", t"6M", t"12M"], [0.8, 0.9, 1.1, 1.15, 1.2], logy=true))
@test read(io, String) == "x = [2, 7, 90, 180, 365], y = [0.8, 0.9, 1.1, 1.15, 1.2], logx = false, logy = true"


text:
@testset "abstract_gp.jl" begin
    k = SqExponentialKernel()
    f = GP(k)
    @test_throws ErrorException mean(f)
    @test_throws ErrorException var(f)
    @test_throws ErrorException cov(f)
    @test_throws ErrorException mean_and_var(f)
    @test_throws ErrorException mean_and_cov(f)
end


text:
mutable struct List elt :: Any; next; end; # structure de liste
# Une liste vide va être : List(Vnothing, Vnothing)

"""function element from pjulia-stdlib:

element(l, n) returns the n-th element of list l."""
function element(x :: List, n :: Int64) :: Any
    a = x
    atbound = false
    for i = 0:n
        if a.next.next == nothing && !atbound
            println("Warning : index outside bound");
        end;
        if a.next.next == nothing
            atbound = true;
        else a = a.next end;
    end;
    return a.elt;
end;

"""function len from pjulia-stdlib:

len(l) returns the length of list l
"""
function len(x :: List) :: Int64
    if x.next == nothing
        return 0;
    else return (1 + len(x.next));
    end;
end;


text:
using MCMCBenchmarkSuite, Distributed
setprocs(4)

ProjDir = @__DIR__
cd(ProjDir)

isdir("tmp") && rm("tmp", recursive=true)
mkdir("tmp")
!isdir("results") && mkdir("results")
path = pathof(MCMCBenchmarkSuite)
@everywhere begin
  using MCMCBenchmarkSuite
  #Model and configuration patterns for each sampler are located in a
  #seperate model file.
  include(joinpath($path, "../Linear_Regression/Linear_Regression_Models.jl"))
end

#run this on primary processor to create tmp folder
include(joinpath(path, "../Linear_Regression/Linear_Regression_Models.jl"))

@everywhere Turing.turnprogress(false)
#set seeds on each processor
seeds = (939388, 39884, 28484, 495858, 544443)
for (i,seed) in enumerate(seeds)
    @fetch @spawnat i Random.seed!(seed)
end
#create a sampler object or a tuple of sampler objects

samplers=(
  CmdStanNUTS(CmdStanConfig, ProjDir),
  AHMCNUTS(AHMCregression, AHMCconfig),
  DHMCNUTS(sampleDHMC)
)

stanSampler = CmdStanNUTS(CmdStanConfig, ProjDir)
#Initialize model files for each instance of stan
initStan(stanSampler)

#Number of data points
Nd = [10, 100, 1000]

#Number of coefficients
Nc = [2, 3]

#Number of simulations
Nreps = 50

options = (Nsamples=2000, Nadapt=1000, delta=.8, Nd=Nd, Nc=Nc)
#perform the benchmark
results = pbenchmark(samplers, simulateRegression, Nreps; options...)

#save results
save(results, ProjDir)

# Make plots
include("primary_plots.jl")


text:
# user functions

"""
    identitypose6fg()

Initialize a and exiting or new factor graph with a first pose and prior as specified by default
keywords or user.
"""
function identitypose6fg(;fg=nothing,
      N::Int=100,
      initCov::Array{Float64,2}=0.001*eye(6),
      initpose::SE3=SE3(0) )
  #
  fg = fg != nothing ? fg : initfg()
  N = 100

  println("Adding PriorPose3 and :x0 to graph...")
  # pts = 0.001*randn(6,N)
  initPosePrior = PriorPose3( MvNormal(veeEuler(initpose), initCov)  )
  pts = getSample(initPosePrior,N)[1]
  v0 = addNode!(fg, :x0,  pts,  N=N, labels=["POSE"])
  f1  = addFactor!(fg,[v0], initPosePrior)
  return fg
end

# function solveandvisualize(fg::FactorGraph,
#       vc::DrakeVisualizer.Visualizer;
#       refresh::Number=1.0,
#       densitymeshes::Vector{Symbol}=Symbol[],
#       drawlandms::Bool=true,
#       N::Int=100,
#       drawtype::Symbol=:max,
#       gt::Dict{Symbol, Tuple{Symbol,Vector{Float64}}}=Dict{Symbol, Tuple{Symbol,Vector{Float64}}}() )
#   #
#   deletemeshes!(vc)
#
#   # draw while solving[1]==true
#   solving = [true]
#   @async begin
#     while solving[1]
#       println(".")
#       visualizeallposes!(vc, fg, drawlandms=drawlandms, drawtype=drawtype, gt=gt)
#       i = 1
#       deletemeshes!(vc)
#       for dm in densitymeshes
#         i+=1
#         visualizeDensityMesh!(vc, fg, dm)
#       end
#       sleep(1)
#     end
#   end
#
#   # solve
#   tree = wipeBuildNewTree!(fg)
#   @time inferOverTree!(fg, tree, N=N)
#   solving[1]=false;
#   nothing
# end

function hasval(d::Dict, va)
  for (k,v) in d
    if v == va
      return true
    end
  end
  return false
end

function projectrbe(fgl::FactorGraph, from::Symbol, to::Vector{Float64})
  x1 = getKDEMax(getVertKDE(fgl,from))
  x2 = to

  wTc = SE3(x1[1:3],Euler(x1[4:6]...))

  cX2 = (matrix(wTc) \ [x2[1:3];1.0])[1:3]

  # Dx = x2[1:3]-x1[1:3]
  Dx = cX2
  range = norm(Dx)
  bearing = atan2(Dx[2],Dx[1])
  elev = -atan2(Dx[3],norm(Dx[1:2]))
  return range, bearing, elev
end

function projectrbe(fgl::FactorGraph, from::Symbol, to::Symbol)
  # x1 = getKDEMax(getVertKDE(fg,from))
  x2 = getKDEMax(getVertKDE(fgl,to))
  return projectrbe(fgl, from, x2[1:3])

  # wTc = SE3(x1[1:3],Euler(x1[4:6]...))
  #
  # cX2 = (matrix(wTc) \ [x2;1.0])[1:3]
  #
  # # Dx = x2[1:3]-x1[1:3]
  # Dx = cX2
  # range = norm(Dx)
  # bearing = atan2(Dx[2],Dx[1])
  # elev = -atan2(Dx[3],norm(Dx[1:2]))
  # return range, bearing, elev
end

#


text:
#=
Given an array of integers, find the first missing positive integer in linear time and constant space. In other words, find the lowest positive integer that does not exist in the array. The array can contain duplicates and negative numbers as well.

For example, the input [3, 4, -1, 1] should give 2. The input [1, 2, 0] should give 3.
=#

function lowestmissing(X::Array{Int64})
    out = 1
    while out in X
        out = out +1
    end
    return out
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using ACEbonds, ACE, ACEbase, Test, StaticArrays, LinearAlgebra, JuLIP
using ACEbase.Testing
using ACEbonds: CylindricalBondEnvelope
using ACE: discrete_jacobi, Rn1pBasis, scal1pbasis, Scal1pBasis,
           evaluate, evaluate_d, Trig1pBasis, @λ,
           Product1pBasis, Categorical1pBasis, SimpleSparseBasis,
           SymmetricBasis, Invariant, PIBasis

## basics

@info("Test housholderreflection")
for _ = 1:15
   rr0 = randn(SVector{3, Float64})
   H = ACEbonds.housholderreflection(rr0)
   print_tf(@test( H * rr0 ≈ norm(rr0) * [0,0,1] ))
   print_tf(@test( H' * [0,0,1] ≈ rr0/norm(rr0) ))
   print_tf(@test H' * H ≈ I)
end
println()

##

@info("Test housholderreflection near the poles")
for _ = 1:15
   rr0 = SVector{3}([0, 0, 1]) + 1e-15 * randn(SVector{3, Float64})
   H = ACEbonds.housholderreflection(rr0)
   print_tf(@test( H * rr0 ≈ norm(rr0) * [0,0,1] ))
   print_tf(@test( H' * [0,0,1] ≈ rr0/norm(rr0) ))
   print_tf(@test H' * H ≈ I)
end
println()

println_slim(@test( ACEbonds.housholderreflection(SVector{3, Float64}([0, 0, 1])) == I ))

##

@info("test transformation from euclidean to cylindrical environments")

# TODO: HACK ACE States to achieve this without type piracy
#       introduce an `ace_isapprox` which defaults to `isapprox`
#       and then overload that one for Base types
Base.isapprox(s1::Symbol, s2::Symbol) = (s1 == s2)

r0cut = 4.0
rcut = 4.0
zcut = 2.0

for ntest = 1:30
   rr0, Zi, Zj, Rs, Zs, Xs = ACEbonds.rand_env(r0cut, rcut, zcut)
   Xenv = ACEbonds.eucl2cyl(rr0, Zi, Zj, Rs, Zs)
   print_tf(@test( all(Xenv[2:end] .≈ Xs) ))
end
println()

##

r0cut = 4.0
rcut = 4.0
zcut = 2.0

maxdeg = 8
maxL = 5

Jz = discrete_jacobi(maxdeg; pcut = 2, pin = 2, xcut = r0cut+zcut, xin = -r0cut-zcut)
Jr = discrete_jacobi(maxdeg; pcut = 2, pin = 2, xcut = rcut, xin = -rcut)
Jr0 = discrete_jacobi(maxdeg; pcut = 2, pin = 2, xcut = r0cut, xin = -r0cut)

Cbe = Categorical1pBasis([:bond, :env], :be, :be, "Cbe")
Zm = Scal1pBasis(:z, nothing, :m, Jz, "Zm")
Rn = Scal1pBasis(:r, nothing, :n, Jr, "Rn")
El = Trig1pBasis(maxL; varsym = :θ, lsym = :l, label = "El")
Pk = Scal1pBasis(:rij, nothing, :k, Jr0, "Jk")

B1p = Product1pBasis((Cbe, Pk, Rn, El, Zm))

##
B1p = Product1pBasis((Cbe, Pk, Rn, El, Zm))
Bsel = ACEbonds.SparseBondBasis(; maxorder = 3,
                                  default_maxdeg = maxdeg,
                                  weight = Dict{Symbol, Float64}(:m => 1.0, :n => 1.0, :k => 1.0, :l => 1.0),
                                 )
ACE.init1pspec!(B1p, Bsel)
length(B1p)

##

# basis = PIBasis(B1p, Bsel; isreal=true)
basis = SymmetricBasis(ACE.Invariant(), B1p, ACE.NoSym(), Bsel; isreal=true)
@show length(basis)

pot = ACE.LinearACEModel(basis)
θ = randn(ACE.nparams(pot)) ./ (1:ACE.nparams(pot)).^2
ACE.set_params!(pot, θ)

##

rr0, Zi, Zj, Rs, Zs, Xs = ACEbonds.rand_env(r0cut, rcut, zcut)
Xenv = ACEbonds.eucl2cyl(rr0, Zi, Zj, Rs, Zs)

evaluate(B1p, Xenv)
b = evaluate(basis, ACEConfig(Xenv))
v = evaluate(pot, ACEConfig(Xenv))
println_slim(@test (dot(ACE.val.(b), θ) ≈ v.val))

##

@info("Test invariance of the new basis")
for ntest = 1:30
   local rr0, Rs, Zs, Xs, Xenv, Zi, Zj
   rr0, Zi, Zj, Rs, Zs, Xs = ACEbonds.rand_env(r0cut, rcut, zcut)
   Xenv = ACEbonds.eucl2cyl(rr0, Zi, Zj, Rs, Zs)
   B1 = evaluate(basis, ACEConfig(Xenv))
   Q = ACE.Random.rand_rot()
   rr0_Q = Q * rr0
   Rs_Q = Ref(Q) .* Rs
   Xenv_Q = ACEbonds.eucl2cyl(rr0_Q, Zi, Zj, Rs_Q, Zs)
   B2 = evaluate(basis, ACEConfig(Xenv_Q))
   print_tf(@test( B1 ≈ B2 && !all(Xenv_Q .≈ Xenv) ))
end
println()

##

@info("Check derivatives of basis w.r.t. cylindrical coordinates")
B = evaluate(basis, ACEConfig(Xenv))
dB = evaluate_d(basis, ACEConfig(Xenv))
TDX = ACE.dstate_type(Xenv[1])

for ntest = 1:30
   U = [ randn(TDX) for _ = 1:length(Xenv) ]
   V = randn(length(B)) ./ (1:length(B))

   F = t -> dot(V, ACE.val.(evaluate(basis, ACEConfig(Xenv + t * U))))
   dF = t -> ( dB = evaluate_d(basis, ACEConfig(Xenv + t * U));
               ACE.contract(sum(V[i] * dB[i, :] for i = 1:length(V)), U) )
   F(0.0)
   dF(0.0)

   print_tf(@test all(ACEbase.Testing.fdtest(F, dF, 0.0; verbose=false)) )
end
println()

##

@info("Check derivatives of basis w.r.t. euclidean coordinates")

for ntest = 1:30
   local rr0, Rs, Zs, Xs, Xenv, Zi, Zj
   rr0, Zi, Zj, Rs, Zs, Xs = ACEbonds.rand_env(r0cut, rcut, zcut)
   Us = randn(SVector{3, Float64}, length(Rs))
   uu0 = randn(SVector{3, Float64})
   V = randn(length(B)) ./ (1:length(B))

   julip2ace = t -> ACEConfig(ACEbonds.eucl2cyl(rr0 + t * uu0, Zi, Zj, Rs + t * Us, Zs))
   F = t -> dot(V, ACE.val.(evaluate(basis, julip2ace(t))))

   dF = t -> begin
         dB = evaluate_d(basis, julip2ace(t))
         dB0, dBenv = ACEbonds.rrule_eucl2cyl(rr0, Zi, Zj, Rs, Zs, dB)
         ACE.contract( sum(V[i] * dBenv[i, :] for i = 1:length(V)), Us) +
                  dot( sum(V[i] * dB0[i] for i = 1:length(V)), uu0 )
      end

   print_tf(@test all(ACEbase.Testing.fdtest(F, dF, 0.0; verbose=false)) )
end
println()

##

@info("Check derivatives of a potential w.r.t. cylindrical coordinates")
v = evaluate(pot, ACEConfig(Xenv))
dv = ACE.grad_config(pot, ACEConfig(Xenv))
TDX = ACE.dstate_type(Xenv[1])

for ntest = 1:30
   U = [ randn(TDX) for _ = 1:length(Xenv) ]

   F = t -> evaluate(pot, ACEConfig(Xenv + t * U)).val
   dF = t -> ( dv = ACE.grad_config(pot, ACEConfig(Xenv + t * U));
               ACE.contract(dv, U) )
   F(0.0)
   dF(0.0)

   print_tf(@test all(ACEbase.Testing.fdtest(F, dF, 0.0; verbose=false)) )
end
println()

##

@info("Check derivatives of potential w.r.t. euclidean coordinates")

for ntest = 1:30
   local rr0, Rs, Zs, Xs, Xenv, Zi, Zj
   rr0, Zi, Zj, Rs, Zs, Xs = ACEbonds.rand_env(r0cut, rcut, zcut)
   Us = randn(SVector{3, Float64}, length(Rs))
   uu0 = randn(SVector{3, Float64})

   julip2ace = t -> ACEConfig(ACEbonds.eucl2cyl(rr0 + t * uu0, Zi, Zj, Rs + t * Us, Zs))
   F = t -> evaluate(pot, julip2ace(t)).val

   dF = t -> begin
         dv_cyl = ACE.grad_config(pot, julip2ace(t))
         dv0, dvenv = ACEbonds.rrule_eucl2cyl(rr0, Zi, Zj, Rs, Zs, dv_cyl)
         ACE.contract(dvenv, Us) + dot(dv0, uu0 )
      end

   print_tf(@test all(ACEbase.Testing.fdtest(F, dF, 0.0; verbose=false)) )
end
println()


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Hector Perez, Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

"""
    @process f(arg...) [cycles]

Create a process from a function `f(arg...)`.

Wrap a function and its arguments in a [`Prc`](@ref) and start it with
[`process!`](@ref).

# Arguments
- `f`: a function,
- `arg...`: arguments, the first argument must be an AbstractClock,
- `cycles::Int`: the number of cycles, `f` should run.

# Returns
- an `Int` process id.
"""
macro process(expr, args...)
    @assert expr isa Expr && expr.head == :call "1st term is not a function call!"
    f = expr.args[1] #extract function passed
    c = expr.args[2] #first function arg must be an AbstractClock
    fargs = expr.args[3:end] #extract other function args
    p = :(Prc($f, $(fargs...))) #create Prc struct
    if isempty(args)
        esc(:(process!($c,$p))) #execute process!
    else
        cycles = args[1] #extract cycle
        esc(:(process!($c,$p,$cycles))) #execute process!
    end
end

"""
    @event f(arg...) T t [n]

Schedule a function `f(arg...)` as an event to a clock.

# Arguments
- `f`: function to be executed at event time,
- `arg...`: its arguments, the first argument must be a clock,
- `T`: a [`Timing`](@ref) (at, after, every),
- `t`: a `Number` or a `Distribution`,
- `n::Int`: repetitions for a repeat event.
"""
macro event(expr, args...)
    @assert expr isa Expr && expr.head == :call "1st term is not a function call!"
    f = expr.args[1] #extract function passed
    c = expr.args[2] #first function arg must be an AbstractClock
    fargs = expr.args[2:end] #extract other function args
    ex = :(fun($f, $(fargs...))) #create Action
    if length(args) == 3
        esc(:(event!($c, $ex, $(args[1:2]...), n = $(args[3])))) #execute event!
    else
        esc(:(event!($c, $ex, $(args...)))) #execute event!
    end
end

"""
```
@event f(farg...) c(carg...)
@event f(farg...) ca
```
Schedule a function `f(farg...)` as a conditional event to a clock.

# Arguments
- `f`: function to be executed at event time,
- `farg...`: its arguments, the first argument must be a clock,
- `c`: function to be evaluated at the clock's sample rate, if it
    returns true, the event is triggered,
- `carg...`: arguments to c,
- `ca`: an anyonymous function of the form `()->...` to be evaluated
    at the clock's sample rate, if it returns true, the event
    is triggered.
"""
macro event(expr, cond)
    @assert expr isa Expr && expr.head == :call "1st term is not a function call!"
    f1 = expr.args[1]
    c = expr.args[2]
    f1args = expr.args[2:end]
    ex1 = :(fun($f1, $(f1args...)))
    if cond isa Expr
        if cond.head == :call
            f2 = cond.args[1]
            f2args = cond.args[2:end]
            ex2 = :(fun($f2, $(f2args...)))
        elseif cond.head == :->
            ex2 = :(()->$(cond.args...))
        end
        esc(:(event!($c, $ex1, $ex2)))
    else    # assume it should be a simple timed event
        esc(:(event!($c, $ex1, $cond)))
    end
end

"""
    @periodic f(arg...) [Δt]

Register a function `f(arg...)` for periodic execution at the
clock`s sample rate.

# Arguments
- `f`: function to be executed periodically,
- `arg...`: its arguments, the first argument must be a clock,
- `Δt`: for setting the clock's sample rate.
"""
macro periodic(expr, args...)
    @assert expr isa Expr && expr.head == :call "1st term is not a function call!"
    f = expr.args[1]
    c = expr.args[2]
    fargs = expr.args[2:end]
    ex = :(fun($f, $(fargs...)))
    if length(args) == 3
        esc(:(periodic!($c, $ex, $(args[1:2]...), $(args[3])))) #execute event!
    else
        esc(:(periodic!($c, $ex, $(args...)))) #execute event!
    end
end

"""
```
@delay clk Δt
@delay clk until t
```
Delay a process on a clock `clk` for a time interval `Δt` or until
a time `t`.
"""
macro delay(clk, delay...)
    esc(:(delay!($clk, $(delay...))))
end

"""
    @wait clk f(arg...)

Conditionally wait on a clock `clk` until `f(arg...)` returns true.
"""
macro wait(clk, expr)
    expr.head != :call && error("2nd term is not a function call.")
    f = expr.args[1]
    fargs = expr.args[2:end]
    ex = :(fun($f, $(fargs...)))
    esc(:(wait!($clk, $ex)))
end

"""
    @run! clk t

Run a clock `clk` for a duration `t`.
"""
macro run!(clk, duration)
    esc(:(run!($clk, $duration)))
end


text:
using Distributions
using Plots
using Optim

# Calculate dose-response
# flux reaction:
function reaction(conc,Km,Vmax)
    10 .^Vmax.*conc./(10 .^Km .+ conc)
end

Km = 0
Vmax = 1
conc = collect(0:.01:2)

flux = reaction(conc, Km, Vmax)
plot(conc,flux)

# simulate data with 7% std
conc_sim = [.1,.5,1,2]
flux_sim = reaction(conc_sim,  Km, Vmax)
flux_sim = randn(length(flux_sim))*0.07.*flux_sim+flux_sim
std_sim = flux_sim * .07

scatter!(conc_sim,flux_sim,yerror=std_sim,grid=false,linealpha = 0)

# negative log-likelihood for simulated data
function L(conc,flux,std,p)
    X = reaction(conc,p[1],p[2])
    LL = -sum(-length(X)*log(2π)/2 .- length(X).*log.(abs.(std)) .- (sum((flux .- X).^2) ./ (2std.^2)))
end

# maximum likelihood fit
init = [1.0,1.0]        # start parameters Km,Vmax
optimum = optimize(p -> L(conc_sim,flux_sim,std_sim,p),init,Newton())

# Profile Likelihood for one fixed parameter
function L_1(conc,flux,std,pfix,pn,p)
    if pn <2
        X = reaction(conc,pfix,p)
    else
        X = reaction(conc,p,pfix)
    end
    LL = -sum(-length(X)*log(2π)/2 .- length(X).*log.(abs.(std)) .- (sum((flux .- X).^2) ./ (2std.^2)))
end
#init = [2.0]
#pfix = 0.0
#pn = 1
#opt = optimize(p -> L_1(conc_sim,flux_sim,std_sim,pfix,pn,p),init,NewtonTrustRegion())

# Profile likelihood for pfix and pfit
function prof(optimum,whichPar)
    chi2 = Array{Float64}(undef,1)
    pfit = Array{Float64}(undef,1)
    pfix_collect = Array{Float64}(undef,1)
    step = .01
    if whichPar<2
        fix = 2
    else
        fix = 1
    end
    for sig = -1:2:1
        dev = 0
        n=0
        pfix = optimum.minimizer[fix]+sig*step*n
        while (dev<4 && pfix.>-4 && pfix.<4)
            println(dev)
            println(pfix)
            pfix = optimum.minimizer[fix]+sig*step*n
            init = [optimum.minimizer[whichPar]]
            profopt = optimize(p -> L_1(conc_sim,flux_sim,std_sim,pfix,fix,p),init,Newton())
            dev = (profopt.minimum - optimum.minimum)
            push!(chi2,dev)
            push!(pfit,profopt.minimizer[1])
            push!(pfix_collect,pfix)
            n+=1
        end
    end
    deleteat!(chi2,1)
    deleteat!(pfit,1)
    deleteat!(pfix_collect,1)
    return chi2,pfit,pfix_collect
end

chi2km,Kmfit,Vmaxfix = prof(optimum,1)
chi2vmax,Vmaxfit,Kmfix = prof(optimum,2)

# thresholds
cl = [.68,.95]
thresholds = quantile.(Chisq(1),cl)

# Plot profile likelihoods
p1 = scatter(Vmaxfix,chi2km,ylabel="deviance",xlabel="Vmax", legend= :none)
scatter!([Vmaxfix[1] Vmaxfix[1]],[chi2km[1] chi2km[1]], color = :red)
plot!(Vmaxfix,ones(length(Vmaxfix)).*thresholds[1],linestyle=:dash,label="0.68")
plot!(Vmaxfix,ones(length(Vmaxfix)).*thresholds[2],linestyle=:dash,label="0.95")

p2 = scatter(Kmfix,chi2vmax,xlabel="Km")
scatter!([Kmfix[1] Kmfix[1]],[chi2vmax[1] chi2vmax[1]], color = :red, legend= :none)
plot!(Kmfix,ones(length(Kmfix)).*thresholds[1],linestyle=:dash,label="0.68")
plot!(Kmfix,ones(length(Kmfix)).*thresholds[2],linestyle=:dash,label="0.95")
plot(p1,p2)

# plot parameter paths
p3 = plot(Kmfit,Vmaxfix, title = "Km fitted", ylabel="Vmax", xlabel = "Km")
scatter!([Kmfit[1] Kmfit[1]],[Vmaxfix[1] Vmaxfix[1]], color = :red)
p4 = plot(Kmfix,Vmaxfit, title = "Vmax fitted", xlabel="Km")
scatter!([Kmfix[1] Kmfix[1]],[Vmaxfit[1] Vmaxfit[1]], color = :red)
plot(p3,p4, legend = :none)

# simulate data for lower concentrations
conc_sim = [.01,.05,.08,.1]
flux_sim = reaction(conc_sim,  Km, Vmax)
flux_sim = rand(length(flux_sim))*0.07.*flux_sim+flux_sim
std_sim = flux_sim * .07

plot(conc,reaction(conc, Km, Vmax))
scatter!(conc_sim,flux_sim,yerror=std_sim,grid=false,linealpha = 0)

# Max-likelihood fit
init = [1.0,1.0]        # start parameters Km,Vmax
optimum = optimize(p -> L(conc_sim,flux_sim,std_sim,p),init,Newton())

# Compute profiles
chi2kmlow,Kmfitlow,Vmaxfixlow = prof(optimum,1)
chi2vmaxlow,Vmaxfitlow,Kmfixlow = prof(optimum,2)

# Plot profile likelihoods
p1 = scatter(Vmaxfixlow,chi2kmlow, ylabel="deviance",xlabel="Km")
scatter!([Vmaxfixlow[1] Vmaxfixlow[1]],[chi2kmlow[1] chi2kmlow[1]], color = :red, legend= :none)
plot!(Vmaxfixlow,ones(length(Vmaxfixlow)).*thresholds[1],linestyle=:dash,label="0.68")
plot!(Vmaxfixlow,ones(length(Vmaxfixlow)).*thresholds[2],linestyle=:dash,label="0.95")

p2 = scatter(Vmaxfitlow,chi2vmaxlow,xlabel="Vmax")
scatter!([Vmaxfitlow[1] Vmaxfitlow[1]],[chi2vmaxlow[1] chi2vmaxlow[1]], color = :red, legend= :none)
plot!(Vmaxfitlow,ones(length(Vmaxfitlow)).*thresholds[1],linestyle=:dash,label="0.68")
plot!(Vmaxfitlow,ones(length(Vmaxfitlow)).*thresholds[2],linestyle=:dash,label="0.95")
plot(p1,p2)

# same procedure for large concentrations
conc_sim = [5, 6, 8, 10]
flux_sim = reaction(conc_sim,  Km, Vmax)
flux_sim = rand(length(flux_sim))*0.07.*flux_sim+flux_sim
std_sim = flux_sim * .07

plot(conc,reaction(conc, Km, Vmax))
scatter!(conc_sim,flux_sim,yerror=std_sim,grid=false,linealpha = 0)

# Max-likelihood fit
init = [1.0,1.0]        # start parameters Km,Vmax
optimum = optimize(p -> L(conc_sim,flux_sim,std_sim,p),init,Newton())

# Compute profiles
chi2kmhigh,Kmfithigh,Vmaxfixhigh = prof(optimum,1)
chi2vmaxhigh,Vmaxfithigh,Kmfixhigh = prof(optimum,2)

# Plot profile likelihoods
p1 = scatter(Kmfixhigh,chi2vmaxhigh,ylabel="deviance",xlabel="Vmax")
scatter!([Kmfixhigh[1] Kmfithigh[1]],[chi2vmaxhigh[1] chi2vmaxhigh[1]], color = :red, legend= :none)
plot!(Kmfixhigh,ones(length(Kmfixhigh)).*thresholds[1],linestyle=:dash,label="0.68")
plot!(Kmfixhigh,ones(length(Kmfixhigh)).*thresholds[2],linestyle=:dash,label="0.95")

p2 = scatter(Kmfixhigh,chi2vmaxhigh,xlabel="Km")
scatter!([Kmfixhigh[1] Kmfixhigh[1]],[chi2vmaxhigh[1] chi2vmaxhigh[1]], color = :red, legend= :none)
plot!(Kmfixhigh,ones(length(Kmfixhigh)).*thresholds[1],linestyle=:dash,label="0.68")
plot!(Kmfixhigh,ones(length(Kmfixhigh)).*thresholds[2],linestyle=:dash,label="0.95")
plot(p1,p2)

# All profiles for Km for different concentrations
p3 = scatter(Kmfix,chi2vmax,ylabel="deviance",xlabel="Km", label = "wide range")
scatter!(Kmfixlow,chi2vmaxlow, label = "low conc")
scatter!(Kmfixhigh,chi2vmaxhigh, label = "high conc")
scatter!(Kmfixhigh,ones(length(Kmfixhigh)).*thresholds[1],linestyle=:dash,label="0.68")
scatter!(Kmfixhigh,ones(length(Kmfixhigh)).*thresholds[2],linestyle=:dash,label="0.95")

# All profiles for Vmax for different concentrations
p4 = scatter(Vmaxfix,chi2km,ylabel="deviance",xlabel="vmax", label = "wide range")
scatter!(Vmaxfixlow,chi2kmlow, label = "low conc")
scatter!(Vmaxfixhigh,chi2kmhigh, label = "high conc")
scatter!(Vmaxfix,ones(length(Vmaxfix)).*thresholds[1],linestyle=:dash,label="0.68")
scatter!(Vmaxfix,ones(length(Vmaxfix)).*thresholds[2],linestyle=:dash,label="0.95")
plot(p3,p4)


text:
@testset "initialization strategies" begin
    n = 100
    d = 5
    dummy_data, labels = generate_mvn_with_outliers(d, n, 42, true, false)
    pools = fill(:U, size(dummy_data, 2))
    model = SVDD.VanillaSVDD(dummy_data)

    @testset "RuleOfThumbSilverman" begin
        # see https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.stats.gaussian_kde.html
        expected = (n * (d + 2) / 4.0)^(-1.0 / (d + 4))
        @test expected == SVDD.calculate_gamma(model, SVDD.RuleOfThumbSilverman())
    end

    @testset "RuleOfThumbScott" begin
        # see https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.stats.gaussian_kde.html
        expected = n^(-1.0 / (d + 4))
        @test expected == SVDD.calculate_gamma(model, SVDD.RuleOfThumbScott())
    end

    @testset "TaxErrorEstimate" begin
        # see Tax and Duin 2004, Support Vector Data Description
        outlier_percentage = .1
        expected = 1 / (n*outlier_percentage)
        @test expected == SVDD.calculate_C(model, SVDD.TaxErrorEstimate(outlier_percentage))
    end

    @testset "BoundedTaxErrorEstimate-1" begin
        n = 10
        dummy_data_small, _ = generate_mvn_with_outliers(d, n, 42, true, false)
        lower_bound = 0.03
        upper_bound = 0.97
        outlier_percentage = 0.05
        @test upper_bound ≈ SVDD.calculate_C(SVDD.VanillaSVDD(dummy_data_small), SVDD.BoundedTaxErrorEstimate(outlier_percentage, lower_bound, upper_bound))
    end

    @testset "BoundedTaxErrorEstimate-2" begin
        n = 2000
        dummy_data_small, _ = generate_mvn_with_outliers(d, n, 42, true, false)
        lower_bound = 0.03
        upper_bound = 0.97
        outlier_percentage = 0.05
        @test lower_bound ≈ SVDD.calculate_C(SVDD.VanillaSVDD(dummy_data_small), SVDD.BoundedTaxErrorEstimate(outlier_percentage, lower_bound, upper_bound))
    end

    @testset "SimpleCombinedStrategy" begin
        gamma_expected = SVDD.calculate_gamma(model, SVDD.RuleOfThumbScott())
        C_expected = 0.5
        init_strategy = SVDD.SimpleCombinedStrategy(SVDD.RuleOfThumbScott(), SVDD.FixedCStrategy(C_expected))
        C, kernel = SVDD.get_parameters(model, init_strategy)
        @test C == C_expected
        @test kernel == MLKernels.GaussianKernel(gamma_expected)
    end

    @testset "BinarySearchCStrategy" begin
        model = SVDD.VanillaSVDD(dummy_data)

        target_outlier_percentage = 0.1
        C_init = 0.5
        max_iter = 10
        eps = 0.01

        strategy = SVDD.BinarySearchCStrategy(target_outlier_percentage, C_init, max_iter, eps, TEST_SOLVER)
        gamma = 5.0
        C_actual = SVDD.calculate_C(model, strategy, GaussianKernel(gamma))

        init_strategy = SVDD.FixedParameterInitialization(GaussianKernel(gamma), C_actual)

        SVDD.initialize!(model, init_strategy)
        SVDD.fit!(model, TEST_SOLVER)
        actual_outlier_percentage = countmap(SVDD.classify.(SVDD.predict(model, model.data)))[:outlier] / size(model.data, 2)
        @test abs.(actual_outlier_percentage - target_outlier_percentage) < 0.04
    end

    @testset "GammaFirstCombinedStrategy" begin
        model = SVDD.SVDDneg(dummy_data, pools)
        target_outlier_percentage = 0.15

        gamma_strategy = SVDD.RuleOfThumbSilverman()
        C_strategy = SVDD.BinarySearchCStrategy(target_outlier_percentage, 0.5, 10, 0.01, TEST_SOLVER)
        init_strategy = SVDD.GammaFirstCombinedStrategy(gamma_strategy, C_strategy)

        SVDD.initialize!(model, init_strategy)
        SVDD.fit!(model, TEST_SOLVER)
        actual_outlier_percentage = countmap(SVDD.classify.(SVDD.predict(model, model.data)))[:outlier] / size(model.data, 2)
        @test abs.(actual_outlier_percentage - target_outlier_percentage) < 0.04
    end
end


text:
function launch_ml!(
    obs_file :: String,
    cp_file  :: String,
    symmetry :: String,
    l        :: Lattice,
    r        :: Reduced_lattice,
    m        :: Mesh,
    a        :: Action,
    p_σ      :: NTuple{2, Float64},
    p_Ωs     :: NTuple{5, Float64},
    p_νs     :: NTuple{5, Float64},
    p_Ωt     :: NTuple{5, Float64},
    p_νt     :: NTuple{5, Float64},
    p_χ      :: NTuple{5, Float64},
    lins     :: NTuple{5, Float64},
    bounds   :: NTuple{5, Float64},
    loops    :: Int64,
    Σ_corr   :: Bool,
    Λi       :: Float64,
    Λf       :: Float64,
    dΛi      :: Float64,
    bmin     :: Float64,
    bmax     :: Float64,
    min_eval :: Int64,
    max_eval :: Int64,
    Σ_tol    :: NTuple{2, Float64},
    Γ_tol    :: NTuple{2, Float64},
    χ_tol    :: NTuple{2, Float64},
    ODE_tol  :: NTuple{2, Float64},
    t        :: DateTime,
    t0       :: DateTime,
    cps      :: Vector{Float64},
    wt       :: Float64,
    ct       :: Float64
    ;
    A        :: Float64 = 0.0,
    S        :: Float64 = 0.5
    )        :: Nothing

    # init ODE solver buffers
    da      = get_action_empty(symmetry, r, m, S = S)
    a_stage = get_action_empty(symmetry, r, m, S = S)
    a_inter = get_action_empty(symmetry, r, m, S = S)
    a_err   = get_action_empty(symmetry, r, m, S = S)
    init_action!(l, r, a_inter)
    set_repulsion!(A, a_inter)

    # init left (right part by symmetry) and central part, full loop contribution and self energy corrections
    da_l    = get_action_empty(symmetry, r, m, S = S)
    da_c    = get_action_empty(symmetry, r, m, S = S)
    da_temp = get_action_empty(symmetry, r, m, S = S)
    da_Σ    = get_action_empty(symmetry, r, m, S = S)

    # init buffers for evaluation of rhs
    num_comps = length(a.Γ)
    num_sites = length(r.sites)
    tbuffs    = NTuple{3, Matrix{Float64}}[(zeros(Float64, num_comps, num_sites), zeros(Float64, num_comps, num_sites), zeros(Float64, num_comps, num_sites)) for i in 1 : Threads.nthreads()]
    temps     = Array{Float64, 3}[zeros(Float64, num_sites, num_comps, 4) for i in 1 : Threads.nthreads()]
    corrs     = zeros(Float64, 2, 3, m.num_Ω)

    # init cutoff, step size, monotonicity and correlations
    Λ        = Λi
    dΛ       = dΛi
    monotone = true
    χ        = get_χ_empty(symmetry, r, m)

    # set up required checkpoints
    push!(cps, Λi)
    push!(cps, Λf)
    cps = sort(unique(cps), rev = true)

    for i in eachindex(cps)
        if cps[i] < Λ
            dΛ = min(dΛ, 0.85 * (Λ - cps[i]))
            break
        end
    end

    # compute renormalization group flow
    while Λ > Λf
        println()
        println("ODE step at cutoff Λ / |J| = $(Λ) ...")
        flush(stdout)

        # set eval for integration
        eval = min(max(ceil(Int64, min_eval / Λ), min_eval), max_eval)

        # prepare da and a_err
        replace_with!(da, a)
        replace_with!(a_err, a)

        # compute k1 and parse to da and a_err
        compute_dΣ!(Λ, r, m, a, a_stage, Σ_tol)
        compute_dΓ_ml!(Λ, r, m, loops, a, a_stage, da_l, da_c, da_temp, da_Σ, tbuffs, temps, corrs, eval, Γ_tol)
        if Σ_corr compute_dΣ_corr!(Λ, r, m, a, a_stage, da_Σ, Σ_tol) end
        mult_with_add_to!(a_stage, -2.0 * dΛ / 9.0, da)
        mult_with_add_to!(a_stage, -7.0 * dΛ / 24.0, a_err)

        # compute k2 and parse to da and a_err
        replace_with!(a_inter, a)
        mult_with_add_to!(a_stage, -0.5 * dΛ, a_inter)
        compute_dΣ!(Λ - 0.5 * dΛ, r, m, a_inter, a_stage, Σ_tol)
        compute_dΓ_ml!(Λ - 0.5 * dΛ, r, m, loops, a_inter, a_stage, da_l, da_c, da_temp, da_Σ, tbuffs, temps, corrs, eval, Γ_tol)
        if Σ_corr compute_dΣ_corr!(Λ - 0.5 * dΛ, r, m, a_inter, a_stage, da_Σ, Σ_tol) end
        mult_with_add_to!(a_stage, -1.0 * dΛ / 3.0, da)
        mult_with_add_to!(a_stage, -1.0 * dΛ / 4.0, a_err)

        # compute k3 and parse to da and a_err
        replace_with!(a_inter, a)
        mult_with_add_to!(a_stage, -0.75 * dΛ, a_inter)
        compute_dΣ!(Λ - 0.75 * dΛ, r, m, a_inter, a_stage, Σ_tol)
        compute_dΓ_ml!(Λ - 0.75 * dΛ, r, m, loops, a_inter, a_stage, da_l, da_c, da_temp, da_Σ, tbuffs, temps, corrs, eval, Γ_tol)
        if Σ_corr compute_dΣ_corr!(Λ - 0.75 * dΛ, r, m, a_inter, a_stage, da_Σ, Σ_tol) end
        mult_with_add_to!(a_stage, -4.0 * dΛ / 9.0, da)
        mult_with_add_to!(a_stage, -1.0 * dΛ / 3.0, a_err)

        # compute k4 and parse to a_err
        replace_with!(a_inter, da)
        compute_dΣ!(Λ - dΛ, r, m, a_inter, a_stage, Σ_tol)
        compute_dΓ_ml!(Λ - dΛ, r, m, loops, a_inter, a_stage, da_l, da_c, da_temp, da_Σ, tbuffs, temps, corrs, eval, Γ_tol)
        if Σ_corr compute_dΣ_corr!(Λ - dΛ, r, m, a_inter, a_stage, da_Σ, Σ_tol) end
        mult_with_add_to!(a_stage, -1.0 * dΛ / 8.0, a_err)

        # estimate integration error
        subtract_from!(a_inter, a_err)
        Δ     = get_abs_max(a_err)
        scale = ODE_tol[1] + max(get_abs_max(a_inter), get_abs_max(a)) * ODE_tol[2]
        err   = Δ / scale

        println("   Relative integration error err = $(err).")
        println("   Current vertex maximum Γmax = $(get_abs_max(a_inter)).")
        println("   Performing sanity checks and measurements ...")

        if err <= 1.0 || dΛ <= bmin
            # update cutoff
            Λ -= dΛ

            # update step size
            dΛ = max(bmin, min(bmax * Λ, 0.85 * (1.0 / err)^(1.0 / 3.0) * dΛ))

            # check if we pass by required checkpoint and adjust step size accordingly
            cps_lower = filter(x -> x < 0.0, cps .- Λ)
            Λ_cp      = 0.0

            if length(cps_lower) >= 1
                Λ_cp = Λ + first(cps_lower)
            end

            dΛ = min(dΛ, Λ - Λ_cp)

            # terminate if vertex diverges
            if get_abs_max(a_inter) > max(min(50.0 / Λ, 1000), 10.0)
                println("   Vertex has diverged, terminating solver ...")
                t, monotone = measure(symmetry, obs_file, cp_file, Λ, dΛ, χ, χ_tol, t, t0, r, m, a_inter, wt, 0.0)
                break
            end

            # do measurements and checkpointing
            mk_cp = false

            for i in eachindex(cps)
                if Λ ≈ cps[i]
                    mk_cp = true
                    break
                end
            end

            if mk_cp
                t, monotone = measure(symmetry, obs_file, cp_file, Λ, dΛ, χ, χ_tol, t, t0, r, m, a_inter, wt, 0.0)
            else
                t, monotone = measure(symmetry, obs_file, cp_file, Λ, dΛ, χ, χ_tol, t, t0, r, m, a_inter, wt, ct)
            end

            # terminate if correlations show non-monotonicity
            if monotone == false
                println("   Flowing correlations show non-monotonicity, terminating solver ...")
                t, monotone = measure(symmetry, obs_file, cp_file, Λ, dΛ, χ, χ_tol, t, t0, r, m, a_inter, wt, 0.0)
                break
            end

            # update frequency mesh
            m = resample_from_to(Λ, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, m, a_inter, a, χ)

            if Λ > Λf
                println("Done. Proceeding to next ODE step.")
            end
        else
            # update step size
            dΛ = max(bmin, min(bmax * Λ, 0.85 * (1.0 / err)^(1.0 / 3.0) * dΛ))

            # check if we pass by required checkpoint and adjust step size accordingly
            cps_lower = filter(x -> x < 0.0, cps .- Λ)
            Λ_cp      = 0.0

            if length(cps_lower) >= 1
                Λ_cp = Λ + first(cps_lower)
            end

            dΛ = min(dΛ, Λ - Λ_cp)

            println("Done. Repeating ODE step with smaller dΛ.")
        end
    end

    # open files
    obs = h5open(obs_file, "cw")
    cp  = h5open(cp_file, "cw")

    # set finished flags
    obs["finished"] = true
    cp["finished"]  = true

    # close files
    close(obs)
    close(cp)

    println("Renormalization group flow finished.")
    flush(stdout)

    return nothing
end


text:
module ChainRules

using Base.Broadcast: materialize, materialize!, broadcasted, Broadcasted, broadcastable
using ChainRulesCore
using Compat
using IrrationalConstants: logtwo, logten
using LinearAlgebra
using LinearAlgebra.BLAS
using Random
using RealDot: realdot
using SparseArrays
using Statistics

# Basically everything this package does is overloading these, so we make an exception
# to the normal rule of only overload via `ChainRulesCore.rrule`.
import ChainRulesCore: rrule, frule

# numbers that we know commute under multiplication
const CommutativeMulNumber = Union{Real,Complex}

include("rulesets/Core/core.jl")

include("rulesets/Base/utils.jl")
include("rulesets/Base/nondiff.jl")
include("rulesets/Base/base.jl")
include("rulesets/Base/fastmath_able.jl")
include("rulesets/Base/evalpoly.jl")
include("rulesets/Base/array.jl")
include("rulesets/Base/arraymath.jl")
include("rulesets/Base/indexing.jl")
include("rulesets/Base/sort.jl")
include("rulesets/Base/mapreduce.jl")

include("rulesets/Statistics/statistics.jl")

include("rulesets/LinearAlgebra/utils.jl")
include("rulesets/LinearAlgebra/blas.jl")
include("rulesets/LinearAlgebra/lapack.jl")
include("rulesets/LinearAlgebra/dense.jl")
include("rulesets/LinearAlgebra/norm.jl")
include("rulesets/LinearAlgebra/matfun.jl")
include("rulesets/LinearAlgebra/structured.jl")
include("rulesets/LinearAlgebra/symmetric.jl")
include("rulesets/LinearAlgebra/factorization.jl")
include("rulesets/LinearAlgebra/uniformscaling.jl")

include("rulesets/SparseArrays/sparsematrix.jl")

include("rulesets/Random/random.jl")

end # module


text:
module GmshDiscreteModelsTests

using Gridap
using GridapGmsh
using Gridap.Geometry
using Gridap.Helpers
using Test

function check_interpolation(model)
  order = 2
  reffe = ReferenceFE(lagrangian,Float64,order)
  V = FESpace(model,reffe)
  if num_dims(model) == 2
    u = x -> x[1]^2 + 2*x[2]^2
  elseif num_dims(model) == 3
    u = x -> x[1]^2 + 2*x[2]^2 + 3*x[3]^2
  else
    @unreachable
  end
  uh = interpolate(u,V)
  e = u - uh
  Ω = Triangulation(model)
  dΩ = Measure(Ω,2*order)
  @test sum(∫( e*e )dΩ) < 1.0e-9
end

mshfile = joinpath(@__DIR__,"t1.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)
writevtk(model,"model")
check_interpolation(model)

mshfile = joinpath(@__DIR__,"twoTetraeder.msh")
model = GmshDiscreteModel(mshfile; renumber=true)
test_discrete_model(model)
check_interpolation(model)

mshfile = joinpath(@__DIR__,"..","demo","demo.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)
check_interpolation(model)

mshfile = joinpath(@__DIR__,"square.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)
check_interpolation(model)

mshfile = joinpath(@__DIR__,"cube.msh")
model = GmshDiscreteModel(mshfile)
test_discrete_model(model)
check_interpolation(model)

end # module


text:
#=
--- Day 4: Passport Processing ---

You arrive at the airport only to realize that you grabbed your North Pole Credentials instead of your passport. While these documents are extremely similar, North Pole Credentials aren't issued by a country and therefore aren't actually valid documentation for travel in most of the world.

It seems like you're not the only one having problems, though; a very long line has formed for the automatic passport scanners, and the delay could upset your travel itinerary.

Due to some questionable network security, you realize you might be able to solve both of these problems at the same time.

The automatic passport scanners are slow because they're having trouble detecting which passports have all required fields. The expected fields are as follows:

    byr (Birth Year)
    iyr (Issue Year)
    eyr (Expiration Year)
    hgt (Height)
    hcl (Hair Color)
    ecl (Eye Color)
    pid (Passport ID)
    cid (Country ID)

Passport data is validated in batch files (your puzzle input). Each passport is represented as a sequence of key:value pairs separated by spaces or newlines. Passports are separated by blank lines.

Here is an example batch file containing four passports:

ecl:gry pid:860033327 eyr:2020 hcl:#fffffd
byr:1937 iyr:2017 cid:147 hgt:183cm

iyr:2013 ecl:amb cid:350 eyr:2023 pid:028048884
hcl:#cfa07d byr:1929

hcl:#ae17e1 iyr:2013
eyr:2024
ecl:brn pid:760753108 byr:1931
hgt:179cm

hcl:#cfa07d eyr:2025 pid:166559648
iyr:2011 ecl:brn hgt:59in

The first passport is valid - all eight fields are present. The second passport is invalid - it is missing hgt (the Height field).

The third passport is interesting; the only missing field is cid, so it looks like data from North Pole Credentials, not a passport at all! Surely, nobody would mind if you made the system temporarily ignore missing cid fields. Treat this "passport" as valid.

The fourth passport is missing two fields, cid and byr. Missing cid is fine, but missing any other field is not, so this passport is invalid.

According to the above rules, your improved system would report 2 valid passports.

Count the number of valid passports - those that have all required fields. Treat cid as optional. In your batch file, how many passports are valid?
=#

function read_input()::Array{String}

    passports::Array{String} = []
    current_passport::Array{String} = []
    open("./input", "r") do f
        while ! eof(f)
            s::String = readline(f)
            if s == ""
                pass_info = join(current_passport, " ")
                push!(passports,pass_info)
                current_passport = []
            else
                push!(current_passport,s)
            end
        end
    end

    pass_info = join(current_passport, " ")
    push!(passports,pass_info)

    return passports
end

function count_valid_passports(passports::Array{String})::Int32

    required_fields = ["byr:","iyr:","eyr:","hgt:","hcl:","ecl:","pid:"]
    valid_count = 0

    for passport in passports
        is_valid = true
        for field in required_fields
            if !(occursin(field,passport))
                is_valid=false
                break
            end
        end
        if is_valid
            valid_count+=1
        end
    end

    return valid_count
end

passports = read_input()
println("Total passports: ", length(passports))
println("Valid passports found: ", count_valid_passports(passports))

#=
--- Part Two ---

The line is moving more quickly now, but you overhear airport security talking about how passports with invalid data are getting through. Better add some data validation, quick!

You can continue to ignore the cid field, but each other field has strict rules about what values are valid for automatic validation:

    byr (Birth Year) - four digits; at least 1920 and at most 2002.
    iyr (Issue Year) - four digits; at least 2010 and at most 2020.
    eyr (Expiration Year) - four digits; at least 2020 and at most 2030.
    hgt (Height) - a number followed by either cm or in:
        If cm, the number must be at least 150 and at most 193.
        If in, the number must be at least 59 and at most 76.
    hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f.
    ecl (Eye Color) - exactly one of: amb blu brn gry grn hzl oth.
    pid (Passport ID) - a nine-digit number, including leading zeroes.
    cid (Country ID) - ignored, missing or not.

Your job is to count the passports where all required fields are both present and valid according to the above rules. Here are some example values:

byr valid:   2002
byr invalid: 2003

hgt valid:   60in
hgt valid:   190cm
hgt invalid: 190in
hgt invalid: 190

hcl valid:   #123abc
hcl invalid: #123abz
hcl invalid: 123abc

ecl valid:   brn
ecl invalid: wat

pid valid:   000000001
pid invalid: 0123456789

Here are some invalid passports:

eyr:1972 cid:100
hcl:#18171d ecl:amb hgt:170 pid:186cm iyr:2018 byr:1926

iyr:2019
hcl:#602927 eyr:1967 hgt:170cm
ecl:grn pid:012533040 byr:1946

hcl:dab227 iyr:2012
ecl:brn hgt:182cm pid:021572410 eyr:2020 byr:1992 cid:277

hgt:59cm ecl:zzz
eyr:2038 hcl:74454a iyr:2023
pid:3556412378 byr:2007

Here are some valid passports:

pid:087499704 hgt:74in ecl:grn iyr:2012 eyr:2030 byr:1980
hcl:#623a2f

eyr:2029 ecl:blu cid:129 byr:1989
iyr:2014 pid:896056539 hcl:#a97842 hgt:165cm

hcl:#888785
hgt:164cm byr:2001 iyr:2015 cid:88
pid:545766238 ecl:hzl
eyr:2022

iyr:2010 hgt:158cm hcl:#b6652a ecl:blu byr:1944 eyr:2021 pid:093154719

Count the number of valid passports - those that have all required fields and valid values. Continue to treat cid as optional. In your batch file, how many passports are valid?
=#

function generate_dict_for_each_passport(passports::Array{String})::Array{Dict{String,String}}

    passports_dict::Array{Dict{String,String}} = []
    for passport in passports

        passport_dict::Dict{String,String} = Dict()
        fields = split(passport," ")
        for field in fields
            key_value = split(field,":")
            passport_dict[key_value[1]] = key_value[2]
        end
        push!(passports_dict, passport_dict)
    end
    return passports_dict
end

function validate_n_digit_number(value::String, n::Int, lower_bound::Int, upper_bound::Int)::Bool
    if length(value) != n
        return false end

    value_int = tryparse(Int32,value)
    if value_int == nothing
        return false end

    if (value_int < lower_bound) || (value_int > upper_bound)
        return false end

    return true
end

function validate_byr(passport_dict::Dict{String,String})::Bool
    #byr (Birth Year) - four digits; at least 1920 and at most 2002.
    if !haskey(passport_dict,"byr")
        return false end
    byr = passport_dict["byr"]
    return validate_n_digit_number(byr, 4, 1920, 2002)
end

function validate_iyr(passport_dict::Dict{String,String})::Bool
    #iyr (Issue Year) - four digits; at least 2010 and at most 2020.
    if !haskey(passport_dict,"iyr")
        return false end
    iyr = passport_dict["iyr"]
    return validate_n_digit_number(iyr, 4, 2010, 2020)
end

function validate_eyr(passport_dict::Dict{String,String})::Bool
    #eyr (Expiration Year) - four digits; at least 2020 and at most 2030.
    if !haskey(passport_dict,"eyr")
        return false end
    eyr = passport_dict["eyr"]
    return validate_n_digit_number(eyr, 4, 2020, 2030)
end

function validate_pid(passport_dict::Dict{String,String})::Bool
    #pid (Passport ID) - a nine-digit number, including leading zeroes.
    if !haskey(passport_dict,"pid")
        return false end
    pid = passport_dict["pid"]
    return validate_n_digit_number(pid, 9, 000000000, 999999999)
end

function validate_hgt(passport_dict::Dict{String,String})::Bool
    # hgt (Height) - a number followed by either cm or in:
    #  - If cm, the number must be at least 150 and at most 193.
    #  - If in, the number must be at least 59 and at most 76.

    if !haskey(passport_dict,"hgt")
        return false end
    hgt = passport_dict["hgt"]
    metric = hgt[end-1:end]

    if metric == "cm"
        return validate_n_digit_number(hgt[1:end-2], 3, 150, 193)
    elseif metric == "in"
        return validate_n_digit_number(hgt[1:end-2], 2, 59, 76)
    else
        return false
    end
end

function validate_hcl(passport_dict::Dict{String,String})::Bool
    #hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f.
    if !haskey(passport_dict,"hcl")
        return false end
    hcl = passport_dict["hcl"]

    if hcl[1] != '#'
        return false end

    for hex_unit in hcl[2:end]
        hex_unit = Int(hex_unit)
        if !(hex_unit>=48 && hex_unit<=57) && !(hex_unit>=97 && hex_unit<=102)
            return false end
    end
    return true
end

function validate_ecl(passport_dict::Dict{String,String})::Bool
    #ecl (Eye Color) - exactly one of: amb blu brn gry grn hzl oth.
    valid_ecls::Array{String} = ["amb","blu","brn","gry","grn","hzl","oth"]
    if !haskey(passport_dict,"ecl")
        return false end
    ecl = passport_dict["ecl"]
    if !(ecl in valid_ecls)
        return false end
    return true
end

function validate_pass_fields(passport::Dict{String,String})::Bool

    validation_foos = [ validate_byr,validate_ecl,validate_eyr,validate_hcl,
                        validate_hgt,validate_pid,validate_iyr]

    for foo in validation_foos
        if !foo(passport)
            return false end
    end
    return true
end

function count_valid_passports(passports::Array{Dict{String,String}})::Int32
    valid_count = 0
    for passport in passports
        if validate_pass_fields(passport)
            valid_count +=1 end
    end
    return valid_count
end

passports = generate_dict_for_each_passport(passports)
valid_count = count_valid_passports(passports)
println("Valid passports considering field rules: ", valid_count)