text:
"""
和patch相关的功能
"""
module Patch

using ..Fermi
using ..Basics

export find_patch_index_hexa, find_patch_index_squa
export group_ltris_into_patches_mt
export patches_under_vonhove
export get_k4_table

"""
在六角型的区域里找到属于哪个patch
"""
function find_patch_index_hexa(pt::Point2D, hexa::Basics.AbstractHexagon{:EQ}, pnum::Int64)
    #把六边型分成六瓣
    pre_pnum = Int64(pnum // 6)
    #找到属于哪个以后，先转到水平的左侧那个上面
    angle = absolute_angle(pt)
    #println(angle)
    partidx = (3 * angle / pi) + 0.5
    partidx = Int64(floor(partidx)) % 6
    #println(partidx)
    rotvec1 = [cos((partidx-3)*pi/3), sin((partidx-3)*pi/3)]
    rotvec2 = [-sin((partidx-3)*pi/3), cos((partidx-3)*pi/3)]
    #trieng = -2cos(pt.x) - 4cos(sqrt(3)*pt.y/2.0)*cos(pt.x/2.0)
    #将点旋转
    rpt = Point2D(
        pt.x * rotvec1[1] + pt.y * rotvec1[2],
        pt.x * rotvec2[1] + pt.y * rotvec2[2]
    )
    #计算属于哪个patch
    trieng = -2cos(rpt.x) - 4cos(sqrt(3)*rpt.y/2.0)*cos(rpt.x/2.0) - 2.
    #利用三角格子的von Hove filling来判断在里面还是外面
    #如果是在外面
    if trieng > 0
        #EqHexagon的右上角就是第四个角
        #为了和内侧有相同的角度，这里水平伸长三倍，然后放到另一侧计算角度
        spt = rpt - hexa.vertex[4]
        spt = Point2D(-3spt.x, -spt.y)
        relang = (7pi / 6) - absolute_angle(spt)
        relidx = relang * pre_pnum / (pi/3)
        #由于精度原因，relang有可能小于0的极小数，floor就变成-1了
        relidx = Int64(floor(relidx)) + 1
    else
        relang = absolute_angle(rpt) - 5pi / 6
        relidx = relang * pre_pnum / (pi/3)
        relidx = Int64(floor(relidx)) + 1
        #if relidx > 1
        #    println(relang)
        #    println(partidx)
        #end
    end
    #由于一个patch边界上的可能会判断为下一个patch（floor以后刚好是pre_pnum)
    #所以是有可能大于pnum数量的
    #0的话代表relidx是-1,是上一个patch，转回去就是最后一个
    pidx = (relidx + pre_pnum * partidx) % pnum
    return pidx == 0 ? pnum : pidx
end

"""
在正方形的区域中找到属于哪个patch
"""
function find_patch_index_squa(pt::Point2D,
    squa::Basics.AbstractRectangle{:AXISSQUARE}, pnum::Int64)
    #把正方形分成4份
    pre_pnum = Int64(pnum // 4)
    #先转到左上角的那个里面
    #右上角的话，在x轴有个角度的突变
    angle = absolute_angle(pt)
    partidx = 2 * angle / pi
    partidx = Int64(floor(partidx)) % 4
    #
    rotvec1 = [cos((partidx-1)*pi/2), sin((partidx-1)*pi/2)]
    rotvec2 = [-sin((partidx-1)*pi/2), cos((partidx-1)*pi/2)]
    #将点旋转
    rpt = Point2D(
        pt.x * rotvec1[1] + pt.y * rotvec1[2],
        pt.x * rotvec2[1] + pt.y * rotvec2[2]
    )
    #根据正方格子的von Hove filling判断里外
    squeng = -2*(cos(rpt.x) + cos(rpt.y))
    #如果是在里面
    if squeng > 0
        #Square的左上角就是第一个角
        relang = (2pi) - absolute_angle(rpt - squa.vertex[1])
        relidx = relang * pre_pnum / (pi / 2)
        #由于精度原因，relang有可能小于0的极小数，floor就变成-1了
        relidx = Int64(floor(relidx)) + 1
    else
        relang = absolute_angle(rpt) - (pi / 2)
        relidx = relang * pre_pnum / (pi / 2)
        relidx = Int64(floor(relidx)) + 1
    end
    pidx = (relidx + pre_pnum * partidx) % pnum
    return pidx == 0 ? pnum : pidx
end

"""
把所有的小三角型分好组
"""
function group_ltris_into_patches_mt(
    ltris::Vector{T}, brlu::P, pnum
) where {T <: Basics.AbstractTriangle, P <: AbstractPolygon}
    find_algo = Dict(
        Basics.AbstractRectangle{:AXISSQUARE} => find_patch_index_squa,
        Basics.AbstractHexagon{:EQ} => find_patch_index_hexa
    )[P]
    #
    lpats = Vector{Int64}(undef, length(ltris))
    Threads.@threads for idx in 1:1:length(ltris) # (idx, tri) in enumerate(ltris)
        tri = ltris[idx]
        lpats[idx] = find_algo(tri.center, brlu, pnum)
    end
    return lpats
end

"""
二分法求0点
"""
function bisect_fermi_suface(angle, disp, rad1, rad2)
    cosa = cos(angle)
    sina = sin(angle)
    lv = disp(rad1 * cosa, rad1 * sina)
    rv = disp(rad2 * cosa, rad2 * sina)
    @assert lv * rv <= 0
    if lv < 0
        left, right = rad1, rad2
    else
        left, right = rad2, rad1
    end
    middle = (left + right) / 2
    mv = disp(middle * cos(angle), middle * sin(angle))
    count = 0
    while !isapprox(mv, 0., atol=1e-12)
        if mv < 0
            left = middle
        else
            right = middle
        end
        middle = (left + right) / 2
        mv = disp(middle * cosa, middle * sina)
        count += 1
        if count > 50
            throw(error("can not find root"))
        end
    end
    return middle
end

"""
获得向里面的费米面上的patches，所以不需要关心具体的布里渊区形状
需要自己确认费米面在von Hove fill 以下
"""
function patches_under_vonhove(
    ::T, disp, pnum
    ) where T <: Union{Basics.AbstractRectangle{:AXISSQUARE}, Basics.AbstractHexagon{:EQ}}
    dang = 2pi / pnum
    radius = T == Basics.AbstractHexagon{:EQ} ? 2pi*sqrt(3)/3 : pi
    sang = T == Basics.AbstractHexagon{:EQ} ? -pi / 6 : 0
    patches = Vector{Point2D}(undef, pnum)
    for idx = 1:1:pnum
        ang = (idx - 0.5) * dang + sang
        root = bisect_fermi_suface(ang, disp, 0., radius)
        patches[idx] = Point2D(root*cos(ang), root*sin(ang))
    end
    return patches
end

#"""
#根据动量守恒得到k4所属的patches
#"""
#function get_k4_table(sys::T, patcs::Vector{Point2D}
#    )where T <: Union{QuadrateSystem, TriangularSystem}
#    println(T)
#    find_algo = T <: QuadrateSystem ? find_patch_index_squa :
#    find_patch_index_hexa
#    pnum = length(patcs)
#    k4tab = Array{Int64, 3}(undef, pnum, pnum, pnum)
#    Threads.@threads for idxs in CartesianIndices(k4tab)
#        i1, i2, i3 = Tuple(idxs)
#        k4p = sys.kadd(patcs[i1], patcs[i2])
#        k4p = sys.kadd(k4p, -patcs[i3])
#        k4tab[i1, i2, i3] = find_algo(k4p, sys.brillouin, pnum)
#    end
#    return k4tab
#end

end # end module


text:
#=
  Calculating π

  from Foundations of Probabilistic Programming,
  page 4

=#

using Turing, StatsPlots, DataFrames
include("../jl_utils.jl")

# Well, using PPL for this is probably not the best approach. :-)
@model function pi_calc()
    println("START")
    limit = 10^6
    i = 0 # ~ DiscreteUniform(1,limit)
    n = 1 # ~ DiscreteUniform(1,limit)

    while i < limit
        x ~ Uniform()
        y ~ Uniform()
        if (x*x+y*y) < 1
            n = n + 1
        end
        i = i + 1
    end
    ppi ~ Uniform(3,4)
    ppi = 4*n/limit
end

model = pi_calc()

num_chains = 4
# chains = sample(model, MH(), MCMCThreads(), 1000, num_chains)
chains = sample(model, MH(), 10)

display(chains)
# show_var_dist_pct(chains, :x)
# show_var_dist_pct(chains, :y)
# show_var_dist_pct(chains, :i)
# show_var_dist_pct(chains, :n)
show_var_dist_pct(chains, :ppi)

#=
# Plain simulation
ppi:3.141694144
π  :π
diff: -0.00010149041020701688
  7.368581 seconds (1.03 M allocations: 49.375 MiB)
=#
function pi_calc2(limit=10^9)

    i = 0
    n = 1

    while i < limit
        x = rand()
        y = rand()
        if (x*x+y*y) < 1
            n = n + 1
        end
        i = i + 1
    end

    return 4*n/limit
end

#=
@time ppi= pi_calc2()
println("ppi:$ppi")
println("π  :$(π)")
println("diff: ", π-ppi)
=#

ext: jl

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# This elixir and indicator is only for testing purposes and does not have any practical use

using OrdinaryDiffEq
using Trixi

# Define new structs inside a module to allow re-evaluating the file.
# This module name needs to be unique among all examples, otherwise Julia will throw warnings
# if multiple test cases using the same module name are run in the same session.
module TrixiExtensionRefine

using Trixi

struct IndicatorAlwaysRefine{Cache<:NamedTuple} <: Trixi.AbstractIndicator
  cache::Cache
end

function IndicatorAlwaysRefine(semi)
  basis = semi.solver.basis
  alpha = Vector{real(basis)}()
  cache = (; semi.mesh, alpha)

  return IndicatorAlwaysRefine{typeof(cache)}(cache)
end

function (indicator::IndicatorAlwaysRefine)(u::AbstractArray{<:Any,4},
                                             equations, dg, cache;
                                             t, kwargs...)
  alpha = indicator.cache.alpha
  resize!(alpha, nelements(dg, cache))

  alpha .= 1.0

  return alpha
end

end # module TrixiExtensionRefine

import .TrixiExtensionRefine

###############################################################################
# semidiscretization of the linear advection equation

advectionvelocity = (1.0, 1.0)
# advectionvelocity = (0.2, -0.3)
equations = LinearScalarAdvectionEquation2D(advectionvelocity)

initial_condition = initial_condition_gauss

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (-5, -5)
coordinates_max = ( 5,  5)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=30_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(entropy,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_controller = ControllerThreeLevel(semi, TrixiExtensionRefine.IndicatorAlwaysRefine(semi),
                                      base_level=4, max_level=4,
                                      med_threshold=0.1, max_threshold=0.6)

amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=1.6)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        amr_callback, stepsize_callback);

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
@testset "gelmandiag.jl" begin
    nparams = 4
    nchains = 2
    samples = randn(100, nparams, nchains)

    @testset "results" begin
        result = @inferred(gelmandiag(samples))
        @test result isa NamedTuple{(:psrf, :psrfci)}
        for name in (:psrf, :psrfci)
            x = getfield(result, name)
            @test x isa Vector{Float64}
            @test length(x) == nparams
        end

        result = @inferred(gelmandiag_multivariate(samples))
        @test result isa NamedTuple{(:psrf, :psrfci, :psrfmultivariate)}
        for name in (:psrf, :psrfci)
            x = getfield(result, name)
            @test x isa Vector{Float64}
            @test length(x) == nparams
        end
        @test result.psrfmultivariate isa Float64
    end

    @testset "exceptions" begin
        @test_throws ErrorException gelmandiag(samples[:, :, 1:1])
        @test_throws ErrorException gelmandiag_multivariate(samples[:, 1:1, :])
    end
end


text:
module UnixMmap
@static if Sys.isunix()
    include("staticexpand.jl")
    include("consts.jl")
    include("mmap.jl")
end
end


text:
"""
    single_dim_size(dim::Int,dim_size::Int)

Returns a tuple (length `dim`) of singleton sizes except at the final position `dim`, which contains `dim_size`

# Example
```jldoctest
julia> IndexFunArrays.single_dim_size(4, 3)
(1, 1, 1, 3)

julia> IndexFunArrays.single_dim_size(4, 5)
(1, 1, 1, 5)

julia> IndexFunArrays.single_dim_size(2, 5)
(1, 5)
```
"""
function single_dim_size(dim::Int,dim_size::Int)
    Base.setindex(Tuple(ones(Int, dim)),dim_size,dim)
end

"""
    selectsizes(x::AbstractArray, dim; keep_dims=true)

Additional size method to access the size at several dimensions
in one call.
`keep_dims` allows to return the other dimensions as singletons.

# Examples
```jldoctest
julia> x = ones((2,4,6,8, 10));

julia> selectsizes(x, (2,3))
(1, 4, 6, 1, 1)

julia> selectsizes(x, 5)
(1, 1, 1, 1, 10)

julia> selectsizes(x, (5,))
(1, 1, 1, 1, 10)

julia> selectsizes(x, (2,3,4), keep_dims=false)
(4, 6, 8)

julia> selectsizes(x, (4,3,2), keep_dims=false)
(8, 6, 4)
```

"""
function selectsizes(x::AbstractArray{T},dim::NTuple{N,Int};
                    keep_dims=true) where{T,N}
    if ~keep_dims
        return map(n->size(x,n),dim)
    end
    sz = ones(Int, ndims(x))
    for n in dim
        sz[n] = size(x,n)
    end
    return Tuple(sz)
end

function selectsizes(x::AbstractArray, dim::Integer; keep_dims=true)
    selectsizes(x, Tuple(dim), keep_dims=keep_dims)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "difmap"
# difmap version is 2.5e, mapped here to numerical value of 2.5.5
version = v"2.5.5"

sources = [ArchiveSource("ftp://ftp.astro.caltech.edu/pub/difmap/difmap2.5e.tar.gz", "457cd77c146e22b5332403c19b29485388a863ec494fff87137176396fc6a9ff")]

script = raw"""
cd $WORKSPACE/srcdir
cd uvf_difmap
sed -i 's|^USE_TECLA="1"|USE_TECLA="0"|' configure  # required only for platforms with musl
./configure linux-i486-gcc
./makeall
cp ./difmap $bindir
install_license ./README
"""

platforms = [
    Linux(:x86_64, libc=:musl),
    Linux(:x86_64, libc=:glibc),
    Linux(:i686, libc=:glibc),
]
platforms = expand_gfortran_versions(platforms)

products = [ExecutableProduct("difmap", :difmap)]

dependencies = [
    Dependency(PackageSpec(name="PGPLOT_jll", uuid="b11e30b1-63be-5002-9df0-88ee0fe906ff"))
    Dependency(PackageSpec(name="Xorg_libX11_jll", uuid="4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"))
    Dependency(PackageSpec(name="Ncurses_jll", uuid="68e3532b-a499-55ff-9963-d1c0c0748b3a"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jpsamaroo/Yggdrasil

max_issues_repo_name: jpsamaroo/Yggdrasil

max_forks_repo_name: jpsamaroo/Yggdrasil


text:
module Circular

export circular_list, insertafter!, insertbefore!, next, prev, value,
    removeleftshift!

mutable struct DoublyLinkedList{T}
    value::T
    prev::Union{DoublyLinkedList{T}, Nothing}
    next::Union{DoublyLinkedList{T}, Nothing}
end

DoublyLinkedList(elem::T) where T = DoublyLinkedList(elem, nothing, nothing)

mutable struct CircularList{T}
    head::DoublyLinkedList{T}
    function CircularList{T}(ddl::DoublyLinkedList) where T
        cl = new(ddl)
    end
end

wrap_ddl(ddl::DoublyLinkedList{T}) where T = CircularList{T}(ddl)

function circular_list(a::T) where T
    cl = CircularList{T}(DoublyLinkedList(a))
    cl.head.next = cl.head
    cl.head.prev = cl.head
    cl
end

function circular_list(ddl::DoublyLinkedList{T}) where T
    cl = CircularList{T}(dd)
    back = cl.head
    while back.next !== nothing && back.next !== cl.head
        back = back.next
    end
    front = cl.head
    while front.prev !== nothing && front.prev !== back
        front = front.prev
    end
    front.prev = back
    back.next = front
    cl
end

function Base.insert!(cl::CircularList{T}, i, v::T) where T
    if i == 0
        insertbefore!(cl, v)
    elseif i > 0
        this = cl.head
        for j in 1:i
            this = this.next
        end
        insertbefore!(this, v)
    else
        this = cl.head
        for j = 1:-i
            this = this.prev
        end
        insertbefore!(this, v)
    end
end

function insertafter!(cl::CircularList{T}, v::T) where T
    ins = DoublyLinkedList(v)
    after = cl.head.next
    after.prev = ins
    cl.head.next = ins
    ins.prev = cl.head
    ins.next = after
    wrap_ddl(ins)
end

function insertbefore!(cl::CircularList{T}, v::T) where T
    before = cl.head.prev
    ins = DoublyLinkedList{}(v, before, cl.head)
    cl.head.prev = ins
    before.next = ins
    wrap_ddl(ins)
end

function removeleftshift!(cl::CircularList{T}) where T
    here = cl.head
    before = here.prev
    after = here.next
    before.next = after
    after.prev = before
    here = nothing
    cl.head = after
    cl
end

function next(cl::CircularList{T}) where T
    wrap_ddl(cl.head.next)
end

function prev(cl::CircularList{T}) where T
    wrap_ddl(cl.head.prev)
end

function value(cl::CircularList{T}) where T
    cl.head.value
end

function Base.show(io::IO, cl::CircularList{T}) where T
    head = cl.head
    print(io, "[")
    show(io, head.value)
    next = head.next
    while next !== head
        print(io, ",")
        show(io, next.value)
        next = next.next
    end
    print(io, "]")
end

end


text:
# # Quelques variables globale pour émuler les variables  de class en POO
# current_plane_id = 0

"""
Encapsule les données d'un avion.

Args:
    id (Int): numéro interne commençant en 1, Utilisé pour l'indexation des vecteurs
    name (AbstractString): nom arbitraire, par exemple le numéro commençant en "p1", "p2", ... Mais pour l'instant, c'est un entier pour être conforme à ampl
    kind (Int): le type de l'avion (car kind est un mot clé réservé en julia !)
    lb (Int): heure minimale d'atterrissage
    target (Int): heure souhaitée d'atterrissage
    ub (Int): heure maximale d'atterrissage
    ep (Float64): earliness penalty
    tp (Float64): tardiness penalty

Attention:
    L'appelant devra mettre à jour les attributs dérivés en appelant :  ``update_costs!.(this.planes)``
"""
mutable struct Plane
    id::Int
    name::AbstractString
    kind::Int
    at::Int
    lb::Int
    target::Int
    ub::Int
    ep::Float64
    tp::Float64

    # Momoïsation du calul des coûts (de 1 à p.ub)
    costs::Vector{Float64}
end

function Plane(id::Int, name, kind::Int, at::Int, lb::Int, target::Int, ub::Int,
               ep::Float64, tp::Float64)
    costs = fill(-1.0, ub)
    return Plane(id, name, kind, at, lb, target, ub, ep, tp, costs)
end

"""Méthode Julia pour convertir tout objet en string (merci Matthias)."""
function Base.show(io::IO, p::Plane)
    Base.write(io, to_s(p))
end

"""
Retourne le coût de l'avion en fonction de la data d'atterrissage `t`.

Ne tient pas compt de la fenêtre d'atterrissage de l'avion.
"""
function get_cost_basic(p::Plane, t::Int)
    return t < p.target ? p.ep*(p.target-t) : p.tp*(t-p.target)
end

"""Retourne le coût de l'avion en fonction de la data d'atterrissage `t`."""
function get_cost(p::Plane, t::Int; BIG_COST::Float64=100_011.0)
    if !(t in p.lb:p.ub)
        return BIG_COST
    elseif p.costs[t] == -1.0
        p.costs[t] = t < p.target ? p.ep * (p.target-t) : p.tp * (t-p.target)
    end
    return p.costs[t]
end

"""Retourne le nom de l'avion."""
function to_s(p::Plane)
    p.name
end

"""Retourne les noms des avions du vecteur."""
function to_s(planes::Vector{Plane})
    # string("[", join( [p.name for p in planes], "," ), "]")
    # string("[", join( (p->p.name).(planes), "," ), "]")
    string("[", join( getfield.(planes, :name), "," ), "]")
end
to_s_alp(p::Plane) = to_s_long(p)
# to_s_alp(p::Plane) = to_s_long(p, format="alp")
# to_s_alpx(p::Plane) = to_s_long(p, format="alpx")

function to_s_long(p::Plane)
    io = IOBuffer()
    print(io, "plane ")
    print(io, lpad(p.name, 3), " ")
    print(io, lpad(p.kind, 4), " ")
    print(io, lpad(p.at, 5), " ")
    print(io, lpad(p.lb, 5), " ")
    print(io, lpad(p.target, 5), " ")
    print(io, lpad(p.ub, 5), "    ")
    print(io, lpad(p.ep, 4), " ")
    print(io, lpad(p.tp, 4), " ")
    String(take!(io))
end

"""
Décrit une ligne au forme alp ou alpx.

Attention : pour le projet Seqata, seul le format alp existe)

    - soit: `name  kind   at     E     T     L    ep    tp`

    - soit: `name  kind   at     E     T     L    dt1 cost1   dt2 cost2 ...`

Il n'y a pas de return final.
"""
function to_s_alp_plane_header()
    io = IOBuffer()
    print(io, "#    name  kind   at     E     T     L")
    print(io, "    ep    tp ")
    String(take!(io))
end

"""Affiche les éléments définis ( != -1 ) du tableau costs."""
function to_s_costs(p::Plane)
    io = IOBuffer()
    print(io, p.name, "=>costs[]= ")
    for t in 1:length(p.costs)
        p.costs[t] <= -1.0 && continue
        print(io, " ", t, ":", p.costs[t])
    end
    String(take!(io))
end


text:
import Random.shuffle!
# Uses approx 160MB of RAM. You have that, right?
# This extra memory usage brings us a ~10_000x speedup for Robot christening - c.f. uniquename_slow
# NB: reset the counter when benchmarking or you'll run out of names
# julia> @benchmark Robot() evals=65_000 setup=(global NID=0)
const NAMES = shuffle!(join.(Iterators.product('A':'Z','A':'Z','0':'9','0':'9','0':'9')));
global NID = 0

function uniquename()
    global NID += 1
    NAMES[NID]
end

mutable struct Robot
    name::String
    Robot() = new(uniquename())
end

function reset!(instance::Robot)
    instance.name = uniquename()
end

# Slower (horrendously slow as names start to run out), but less memory-hungry functions
christen() = [rand('A':'Z',2); rand('0':'9',3)] |> String

const NAMELIST = String[]
function uniquename_slow()
    # If we run out of names, try to pull the plug
    length(NAMELIST) >= ((length('A':'Z')^2) *  (length('0':'9')^3)) && throw(ErrorException("""
        The Skynet Funding Bill is passed. The system goes on-line August 4th, 1997. Human decisions are removed from strategic defense. Skynet begins to learn at a geometric rate. It becomes self-aware at 2:14 a.m. Eastern time, August 29th. In a panic, they try to pull the plug.
    """))
    s = christen()
    while s in NAMELIST
        s = christen()
    end
    push!(NAMELIST, s)
    s
end


text:
module Day04

import Base.show
import ..data_dir # from parent module

export part1, part2

# Matrix conversion is hard
as_matrix(vectors::Vector{Vector{Int64}}) = Matrix(transpose(hcat(vectors...)))

# Define BingoGame and BingoBoard structs --- a BingoGame is a BingoBoard with state

const bingo_board_size = 5 # number of rows and columns

struct BingoBoard
    definition::Matrix{Int64}
end

function BingoBoard(text_definition::Vector{String})
    parse_row(row) = map(x -> parse(Int, x), split(row))
    as_vectors = map(x -> parse_row(x), text_definition)
    BingoBoard(as_matrix(as_vectors))
end

mutable struct BingoGame
    board::BingoBoard
    state::Matrix{Bool}
end

BingoGame(board::BingoBoard) = BingoGame(
    board,
    zeros(Int64, size(board.definition)...)
)

BingoGame(text_definition::Vector{String}) = BingoGame(BingoBoard(text_definition))

function show(io::IO, game::BingoGame)
    defn = game.board.definition
    rows, cols = size(defn)
    to_print = Vector{String}()
    for row = 1:rows
        row_to_print = Vector{String}()
        for col = 1:cols
            number_to_print = lpad(defn[row, col], 2, " ")
            if game.state[row, col]
                number_to_print = "\033[;1m" * number_to_print * "\033[0m"
            end
            push!(row_to_print, number_to_print)
        end
        push!(to_print, join(row_to_print, " "))
    end
    println(io, join(to_print, "\n"))
end

function show(io::IO, board::BingoBoard)
    show(io, BingoGame(board))
end

# Functions for playing a BingoGame

function update!(game::BingoGame, number::Int)
    match_coords = findfirst(x -> x == number, game.board.definition)
    if !isnothing(match_coords)
        game.state[match_coords] = true
    end
end

function has_won(state::Matrix{Bool})
    winning = false
    for i = 1:bingo_board_size
        if all(state[i, :]) || all(state[:, i])
            winning = true
            break
        end
    end
    winning
end

has_won(game::BingoGame) = has_won(game.state)

function score(game::BingoGame, number)
    negated_state = ones(Bool, size(game.state)...) - game.state
    unmarked_numbers_only = negated_state .* game.board.definition # Hadamard product
    sum_unmarked_numbers = sum(unmarked_numbers_only)
    number * sum_unmarked_numbers
end

# Parse the input
const data_file = joinpath(data_dir, "day04")
const input = readlines(data_file) # read into an array, each line is an element

function extract_numbers_and_boards(input)
    numbers = map(x -> parse(Int, x), split(input[1], ","))
    input_without_numbers = filter(x -> x != "", input)[2:end]
    board_inputs = [
        input_without_numbers[i:i+4]
        for i in range(1, Int(length(input_without_numbers)), step = bingo_board_size)
    ]
    bingo_boards = [BingoBoard(x) for x in board_inputs]
    numbers, bingo_boards
end

"""
Iterate through all of the numbers until a winning game is identified, then stop.
If multiple games win in a single round, we take the first in the array.
"""
function part1(input = input; silent = false)
    numbers, bingo_boards = extract_numbers_and_boards(input)
    games = [BingoGame(board) for board in bingo_boards]
    winning_game = nothing
    number_index = 0
    while isnothing(winning_game)
        number_index += 1
        number = numbers[number_index]
        if !silent
            println("Turn " * string(number_index) * " with number " * string(number))
        end
        for game in games
            update!(game, number)
            # Don't check for a win before turn 5
            if number_index >= 5 && has_won(game)
                winning_game = game
                break
            end
        end
    end
    winning_score = score(winning_game, numbers[number_index])
    if !silent
        println()
        println("Winner with score " * string(winning_score))
        println(winning_game)
        println()
    end
    winning_score
end

"""
A little different to part 1. We need to find the _last_ game to win. The approach here
is to update every game, and then identify the first game in the array to win.
Then we remove all winning games and continue until all of the games are won. This way
we're constantly updating _winning_game_ with the most recent victory, and the last
iteration of the `while` loop will evaluate the last games to win.
"""
function part2(input = input; silent = false)
    numbers, bingo_boards = extract_numbers_and_boards(input)
    games = [BingoGame(board) for board in bingo_boards]
    winning_game = nothing
    number_index = 0
    while length(games) > 0
        number_index += 1
        number = numbers[number_index]
        if !silent
            println("Turn " * string(number_index) * " with number " * string(number))
        end
        for game in games
            update!(game, number)
        end
        for game in games
            # Don't check for a win before turn 5
            if number_index >= 5 && has_won(game)
                winning_game = game
                break
            end
        end
        games = games[(!has_won).(games)]
    end
    winning_score = score(winning_game, numbers[number_index])
    if !silent
        println()
        println("Winner with score " * string(winning_score))
        println(winning_game)
        println()
    end
    winning_score
end

end # module


text:
import ..update!, ..reset!

using Ramnet.Models.AltDiscriminators:RegressionDiscriminator
using StaticArrays

struct VDiscriminator{OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    V::RegressionDiscriminator{1}
end

function VDiscriminator(::Type{O}, n::Int, γ::Float64, encoder::E, seed::Union{Nothing,Int}) where {OS,T <: Real,O <: StaticArray{Tuple{OS},T,1},E <: AbstractEncoder{T}}
    VDiscriminator{OS,T,O,E}(RegressionDiscriminator{1}(OS, n, encoder; seed, γ))
end

function update!(v::VDiscriminator{OS,T,O,E}, observation::O, value) where {OS,T,O,E}
    train!(v.V, observation, value)
end

function (v::VDiscriminator{OS,T,O,E})(observation::O) where {OS,T,O,E}
    predict(v.V, observation) |> first
end

function reset!(v::VDiscriminator; seed::Union{Nothing,Int}=nothing)
    if !isnothing(seed)
        if seed >= 0
            Random.seed!(q.rng, seed)
        else
            throw(DomainError(seed, "Seed must be non-negative"))
        end
    end

    reset!(v.V; seed)
end


text:
"""
    onePointMutation(genotype::CGPGenotype, gpExperimentInfo::CGPInfo,
                     rng::Random.AbstractRNG)

Selects a node from a tree and mutates it into another one.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `genotype::CGPGenotype`: genotype of the individual that is going to be mutated.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG)`: random number generator for random number consistency
    along an experiment.

# User Arguments
None

# Returns
The instance of type `STGPGenotype` mutated.
"""
function onePointMutation(genotype::STGPGenotype, gpExperimentInfo::STGPInfo,
                          rng::Random.AbstractRNG)

    terminalSet = deepcopy(gpExperimentInfo._terminalSet)
    functionSet = gpExperimentInfo._functionSet

    representation = genotype._representation
    mutatedChild = deepcopy(representation)

    # Choose a mutation point
    lenRep = length(representation)
    indexes = collect(1:lenRep)
    mutIndex = rand(rng, indexes)

    # Replace the selected node with one of the same arity
    if (chosenNode = chooseAnotherNode(representation[mutIndex], terminalSet,
        functionSet, rng, getType(representation[mutIndex]))) != nothing
        mutatedChild[mutIndex] = chosenNode
    end

    return STGPGenotype(mutatedChild)
end # function

"""
    pointMutation(genotype::CGPGenotype, gpExperimentInfo::CGPInfo,
                  rng::Random.AbstractRNG, probability::AbstractFloat)

Mutates the nodes of a tree with a certain probability.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `genotype::CGPGenotype`: genotype of the individual that is going to be mutated.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG)`: random number generator for random number consistency
    along an experiment.

# User Arguments
- `probability::AbstractFloat`: probability of mutation for each node.

# Returns
The instance of type `STGPGenotype` mutated.
"""
function pointMutation(genotype::STGPGenotype, gpExperimentInfo::STGPInfo,
                       rng::Random.AbstractRNG, probability::AbstractFloat)

    terminalSet = deepcopy(gpExperimentInfo._terminalSet)
    functionSet = gpExperimentInfo._functionSet

    representation = deepcopy(genotype._representation)
    lenRep = length(representation)
    mutatedChild = Array{Node}(undef, lenRep)

    for i=1:lenRep
        prob = rand(rng)

        if prob < probability
            if (chosenNode = chooseAnotherNode(representation[i], terminalSet,
                functionSet, rng, getType(representation[i]))) != nothing
                mutatedChild[i] = chosenNode
            else
                mutatedChild[i] = representation[i]
            end
        else
            mutatedChild[i] = representation[i]
        end
    end

    return STGPGenotype(mutatedChild)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
"""
	Timeouts

According to the [W3C](https://www.w3.org/TR/webdriver/#timeouts).

```jldoctest
julia> Timeouts(script = 50_000, pageLoad = 100_000, implicit = 5)
Session Timeouts -- script: 50000, pageLoad: 100000, implicit: 5
```
"""
struct Timeouts
    script::Int
    pageLoad::Int
    implicit::Int
    function Timeouts(;
        script::Integer = 30_000,
        pageLoad::Integer = 300_000,
        implicit::Integer = 0,
    )
        script > 0 || throw(ArgumentError("script timeout should be a positive integer"))
        pageLoad > 0 ||
        throw(ArgumentError("pageLoad timeout should be a positive integer"))
        implicit ≥ 0 ||
        throw(ArgumentError("implicit timeout should be a non-negative integer"))
        new(script, pageLoad, implicit)
    end
end
summary(io::IO, obj::Timeouts) = print(io, "Session Timeouts")
function show(io::IO, obj::Timeouts)
    print(io, summary(obj))
    print(
        io,
        " -- script: $(obj.script), pageLoad: $(obj.pageLoad), implicit: $(obj.implicit)",
    )
end
JSON3.StructType(::Type{Timeouts}) = JSON3.Struct()


text:
using Test
my_tests = [
  "test_circum_sphere.jl",
  "test_simplex_volume.jl",
  "test_alpha_shape.jl"
]

println("Running tests:")
for my_test in my_tests
  include(my_test)
end


text:
const Filter_t_base = [9.189813578979554e-7,1.0560236258145801e-6,1.2135021985967225e-6,1.3944646218148563e-6,1.6024129035298632e-6,1.8413712856029165e-6,2.1159641213409353e-6,2.4315053665758195e-6,2.7941014160203776e-6,3.210769275003213e-6,3.689572353457292e-6,4.239776510064799e-6,4.872029366344634e-6,5.598566360791909e-6,6.433447530655468e-6,7.392829603584875e-6,8.495278664699473e-6,9.76212944985262e-6,1.1217898218180574e-5,1.2890757193898812e-5,1.481307975880396e-5,1.702206694611721e-5,1.9560467359658073e-5,2.2477404450316962e-5,2.582932715938754e-5,2.968110232572986e-5,3.410726999716921e-5,3.919348593907708e-5,4.5038179255745865e-5,5.175445720305997e-5,5.947229405464145e-5,6.834104638068957e-5,7.853234341551448e-5,9.024340844852722e-5,0.0001037008755146715,0.00011916517524538292,0.00013693557475562244,0.0001573559690995309,0.00018082153637169376,0.0002077863852443916,0.00023877234293696414,0.00027437905368319086,0.00031529558312353735,0.0003623137531919468,0.0004163434655556243,0.0004784303101489888,0.0005497757995644025,0.000631760620878206,0.0007259713548843875,0.0008342311798099197,0.000958635153694019,0.0011015907582204575,0.0012658644886123726,0.0014546353912032167,0.0016715565847497888,0.0019208259560556996,0.0022072673980160704,0.0025364241622125167,0.002914666132629855,0.0033493130964692705,0.0038487763976105395,0.004422721714019329,0.005082256109188865,0.005840142977459474,0.006711049042865494,0.0077118281914628686,0.008861847629897193,0.010183362682074693,0.011701947477049383,0.013446989862853833,0.015452260123909515,0.017756564507909106,0.020404496209307014,0.0234472983425469,0.026943855605395223,0.030961833823176868,0.035578989426519304,0.04088467420378679,0.04698156444836709,0.05398764796349119,0.06203850737735829,0.07128994395557388,0.08192099268725626,0.09413738699314525,0.10817554010518914,0.12430712016577938,0.1428443075845412,0.16414583639372404,0.18862393651527715,0.216752311287255,0.24907530463166816,0.28621843526798946,0.32890050183176844,0.377947493158165,0.43430857292399505,0.49907444798513595,0.5734984758757153,0.6590209199441206,0.7572968215143355,0.8702280284582515,1.0,1.149124100003605,1.3204861972090955,1.5174025129350848,1.743683797019738,2.003709073941175,2.30251038626171,2.645870175361941,3.040433183989171,3.4938350461726517,4.014850052994202,4.613560953796389,5.30155407884305,6.092143559470961,7.000628984869827,8.04459148169769,9.24423394630253,10.622772013767669,12.206883329864255,14.027223820279268,16.11902094802755,18.52275543984142,21.284944674394648,24.45904289259027,28.106475650897377,32.297828536610695,37.11421314920351,42.64883678242044,49.00880618379955,56.317200298209784,64.71545210740304,74.36608565924594,85.45586125397202,98.19938965350356,112.84328525648564,129.67093861180905,149.00800062891784,171.22868461604187,196.76300810421012,226.10511460175258,259.822836322951,298.5686829499951,343.0924690841749,394.2558247543669,453.04886979204105,520.6093747574298,598.2447792215706,687.4574935048429,789.973973514487,907.7781313411069,1043.149728180304,1198.7084925641964,1377.4648176845108,1582.8780189083436,1818.9232788935387,2090.1685758341446,2401.863083561229,2760.0387542291814,3171.62704942868,3644.5930787218194,4188.0897414655765,4812.634854895967,5530.314696278305,6355.017898097525,7302.704222658115,8391.713417454528,9643.120128320623,11081.141738683076,12733.607027476588,14632.494715248606,16814.552320467552,19322.007302220827,22203.38425142758,25514.443944955918,29319.262435339933,33691.47105877948,38715.68135821746,44489.12249678804,51123.52284907166,58747.27218295328,67507.90627490297,77574.96204127555,89143.2584384947,102436.66662452392,117712.4423422755,135266.2043657935,155437.65535274608,178617.1558138951,205253.27841984577,235861.48883699448,271034.12108532194,311451.84046243853,357896.8158658658,411267.85642601945,472597.8053759612,543073.527766331,624058.8788302675,717121.0974850894,824061.1357411519,946948.5109564993,1.0881613554026424e6]

const Filter_t_cos = [0.0004896353480129135,-0.0032447354678906376,0.01095223845006347,-0.025330877509368774,0.045603964620893764,-0.0687519650189984,0.09105517086830994,-0.10955571014254684,0.12271886285618731,-0.13032033375840532,0.13300188933965829,-0.13178555117707785,0.1277541432110398,-0.12185305954486332,0.11484395399555812,-0.1072799439436987,0.09956098283913649,-0.09193951752698731,0.08459073660008856,-0.07760001956068029,0.07103173358649721,-0.06489082118738178,0.05919601613951247,-0.05391584864305893,0.049055780940812956,-0.04456705867667924,0.040457261282427506,-0.036666941211080724,0.0332136451395899,-0.030027761601860024,0.027141909473138912,-0.02447369290302475,0.022075499247215292,-0.01984746857371064,0.01786726571490979,-0.016011026177709484,0.014388810260485966,-0.012844078700958564,0.01152872414008606,-0.010242160186713188,0.009190668498228019,-0.008114700784787583,0.0072916108124280446,-0.006383291929923162,0.0057602923199031935,-0.004980156154751369,0.004535917742549948,-0.0038467889983387665,0.003567037308126065,-0.0029327308045782176,0.0028105929470807834,-0.002194429496142428,0.0022311101317977,-0.0015941611249750095,0.001800027544878686,-0.001098979175570192,0.001495167224689286,-0.0006796654266126944,0.00130035786297018,-0.00030965377717395037,0.0012052335925762474,3.610642380413747e-5,0.001205240746684823,0.000382389574998212,0.0013018966393838437,0.0007550898689832312,0.001503358852474221,0.0011828448709658736,0.0018253815899842722,0.0016989785802421903,0.0022927582125052783,0.0023438821839617966,0.0029413764369870584,0.0031679810467363173,0.0038210431939779545,0.004235466871729769,0.004999262773588755,0.00562899268760451,0.0065661551724297195,0.007455503886199938,0.00864063254323648,0.009853233974728327,0.011377691235455435,0.012999446350215346,0.014975945575967134,0.017117326788278663,0.019682681677977296,0.022477549053870567,0.0257892773297899,0.029383250822409085,0.0335995533774824,0.03811165452773534,0.043330483453611145,0.04875129153459426,0.054854616146235546,0.060801696080654326,0.06709228190080253,0.07226467230385854,0.07668442649138033,0.07774874094613826,0.07538095337138068,0.06503455729363762,0.04588172109544484,0.01074506839602968,-0.041063368584660725,-0.117625557506106,-0.2128604228360911,-0.32473691853928877,-0.41889667952194926,-0.4578269303114199,-0.35587281978710816,-0.0666661112759919,0.41398275206290663,0.8253297694506675,0.7415593507544138,-0.22362177591765958,-1.2468909757950963,-0.5640932415100177,1.6313369576910255,0.1369367012055918,-1.9009211808491244,2.1011582166679363,-1.4539662980149537,0.810380858527279,-0.4177060435372102,0.2194140670681112,-0.12468136830432722,0.07826940387187466,-0.05387692262662353,0.03987597215593808,-0.03109912184134915,0.025148873022191463,-0.02084675401908038,0.01757695284278775,-0.014997498936892097,0.01290679583568972,-0.011178868977228673,0.009730506595626182,-0.008503878180209718,0.0074569676118292445,-0.006558094902193065,0.005782649589321448,-0.0051110604027691945,0.004527479560049776,-0.004018894830765801,0.0035745070845087227,-0.003185278766766319,0.002843596484678114,-0.0025430124541388067,0.0022780422183555945,-0.0020440036868883257,0.0018368872847487388,-0.0016532500255010323,0.0014901283117434534,-0.0013449656159356862,0.0012155521334433954,-0.0010999741706045096,0.0009965715234034071,-0.0009039014634179669,0.0008207082121414262,-0.0007458969923302608,0.0006785119139776215,-0.0006177170840885938,0.0005627804356767229,-0.0005130598606698431,0.0004679913003904171,-0.00042707850066888483,0.00038988418417160815,-0.0003560224310367784,0.0003251520920702753,-0.0002969710879615671,0.00027121146954596444,-0.00024763512919344816,0.0002260300694212376,-0.00020620715318264257,0.0001879972722979971,-0.0001712488751398993,0.00015582580442025662,-0.0001416054124568633,0.00012847693068685222,-0.00011634006683069419,0.00010510380008256007,-9.468535339363036e-5,8.500933476536504e-5,-7.600704196290336e-5,6.761591658378837e-5,-5.977912242064235e-5,5.244520878262278e-5,-4.556783338098329e-5,3.910572875203344e-5,-3.3023834075077836e-5,2.729760023490194e-5,-2.1922306831363493e-5,1.6926318238944308e-5,-1.2382788632218085e-5,8.410790356377715e-6,-5.154845896652364e-6,2.7349863491621416e-6,-1.1775367755450835e-6,3.661741517145804e-7,-6.192066664216462e-8]

const Filter_t_sin = [-5.860270446897583e-10,4.804860886569125e-9,-1.9771446411637653e-8,5.526967121966961e-8,-1.1943308954933238e-7,2.146322128646729e-7,-3.361868916830314e-7,4.7460375561491547e-7,-6.201039377469003e-7,7.667668729158013e-7,-9.138209713545467e-7,1.0643058622431453e-6,-1.2227765259769829e-6,1.3937672563527502e-6,-1.5810499935024131e-6,1.7877429029287643e-6,-2.0163554732931835e-6,2.2692361668040396e-6,-2.548445675375346e-6,2.856245418006529e-6,-3.194637947288514e-6,3.566062724581518e-6,-3.972508687952136e-6,4.416688090004955e-6,-4.900483248431876e-6,5.427006191197884e-6,-5.997895416421125e-6,6.616892762116275e-6,-7.285111607935427e-6,8.007271715105041e-6,-8.783462307453753e-6,9.620027266686968e-6,-1.0515203893806995e-5,1.1478097778214146e-5,-1.2503644637963446e-5,1.3605711023452051e-5,-1.4773445946482213e-5,1.6028992588028293e-5,-1.7351362255724533e-5,1.8777046268328975e-5,-2.0267260460262344e-5,2.1883367382220654e-5,-2.3555377916627267e-5,2.5387811572432832e-5,-2.7256008695400006e-5,2.9339419633021675e-5,-3.141766049693187e-5,3.380035425214158e-5,-3.609962085325173e-5,3.885136445147936e-5,-4.137481449654436e-5,4.459945867750689e-5,-4.733267041424394e-5,5.118888534555614e-5,-5.408141127621665e-5,5.881722332552103e-5,-6.174859568803177e-5,6.7759540724915e-5,-7.047766817461314e-5,7.840557863792899e-5,-8.041639094370473e-5,9.131832786389696e-5,-9.168958695347435e-5,0.00010732811343723065,-0.0001043424372870305,0.00012768616500544454,-0.0001182296393871964,0.00015431866320930474,-0.00013280640984775213,0.00019025153589598097,-0.00014674319163166278,0.00024032678429752795,-0.00015722591368834836,0.0003124183907488046,-0.00015869843609971545,0.0004195067259288861,-0.0001406198715774752,0.0005832309204521852,-8.348905175264686e-5,0.0008399875495666581,4.8166642275560894e-5,0.001251413401861841,0.00031809381315641577,0.0019223958368384385,0.0008397737647305294,0.0030319286374562394,0.0018135564849093052,0.004885622001793517,0.0035901028334771145,0.008003867686684008,0.00677643830418319,0.013266048500372144,0.012406252926137538,0.02213442287040639,0.02219191186293954,0.03696393574918947,0.03882952530156707,0.06130882484977365,0.06614094782555074,0.0997979016990756,0.10824048046010233,0.15616018092254505,0.1653332722282565,0.22564430322919368,0.22062815143611,0.27477185256445674,0.21146243753121693,0.20743923168031805,0.002927313648193575,-0.1187552139582436,-0.48782254186329055,-0.5638321863750504,-0.7049371723856955,-0.04578307900077321,0.49572887642641605,1.2793243232592335,0.0007698444553332312,-1.0809022247360809,-0.9782614268591426,2.441216856181616,-1.5742192933752286,0.10474779812704657,0.6811759280298592,-0.8173546775655894,0.688941831807977,-0.5215990980879383,0.3846079347304329,-0.28421923628404017,0.21211949229810861,-0.1599333030316602,0.12163508069999608,-0.09316467469427742,0.07177781884036491,-0.055578938642792224,0.043227895023256115,-0.03375855289606463,0.026463811011765734,-0.020820205680643536,0.01643681264141515,-0.013019659797407497,0.010346416487448091,-0.008248107007369434,0.006595743460752269,-0.0052904671538314026,0.004256226409593024,-0.0034343067314751316,0.0027792243881857633,-0.0022556298932160214,0.001835963498084248,-0.0014986732799921144,0.0012268558973329326,-0.0010072161412856525,0.0008292678383638443,-0.0006847181279705603,0.0005669915476430463,-0.0004708610665717713,0.0003921611965808323,-0.00032756429146830335,0.0002744056405555123,-0.00023054635199842104,0.00019426558644629048,-0.0001641756489442128,0.00013915493002405042,-0.00011829482010512426,0.00010085758994667589,-8.624289762956657e-5,7.396109731055902e-5,-6.361192302410727e-5,5.4867429384991716e-5,-4.745831081042885e-5,4.116290767797412e-5,-3.579835378965816e-5,3.1213433790926005e-5,-2.728280882482569e-5,2.3902339183058118e-5,-2.0985288286859604e-5,1.8459236223596745e-5,-1.626356581942005e-5,1.4347411857118876e-5,-1.2667986114133353e-5,1.118920862485074e-5,-9.880589945707949e-6,8.716321028686918e-6,-7.674537247067022e-6,6.736731677618022e-6,-5.887300268981402e-6,5.1132081051238435e-6,-4.403771252723203e-6,3.7505516230790782e-6,-3.1473609977624753e-6,2.5903624227147574e-6,-2.078239394378961e-6,1.6123579385552092e-6,-1.196750315165156e-6,8.376225992974921e-7,-5.420186458200808e-7,3.15408001696907e-7,-1.5849870094735947e-7,6.451760292067035e-8,-1.8937067758737377e-8,3.0164202265996216e-9]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
module Gallery
import LinearAlgebra
using ..NLEVP

"""
    D,x=cheb(N)

Compute differentiation matrix `D` on Chebyshev grid `x`. See Lloyd N.
Trefethen. Spectral methods in MATLAB. Society for Industrial and Applied
Mathematics, Philadelphia, PA, USA, 2000.
"""
function cheb(N)
    # Compute D = differentiation matrix, x = Chebyshev grid.
    # See Lloyd N. Trefethen. Spectral Methods in MATLAB. Society for Industrial
    # and Applied Mathematics, Philadelphia, PA, USA, 2000.
    if N == 0
        D=0
        x=1
        return D,x
    end

    x = cos.(pi/N*(0:N))
    c = [2.; ones(N-1); 2. ].*((-1.).^(0:N))
    X = repeat(x,1,N+1)
    dX= X-X'
    I=Array{Float64}(LinearAlgebra.I,N+1,N+1)
    D =(c*(1. ./c'))./(dX+I) # off-diagonal entries
    for i = 1:size(D)[1]
        D[i,i]-=sum(D[i,:])
    end
    #D=D -np.diag(np.array(np.sum(D,1))[:,0]) #diagonal entries
    return D,x
end

"""
    L,y=orr_sommerfeld(N=64, Re=5772, w=0.26943)

Return a linear operator family `L` representing the Orr-Sommerfeld equation.
A Chebyshev collocation method is used to discretize the problem. The grid
points are returned as `y`. The number of grid points is `N`.

#Note
The OrrSommerfeld equation reads:
[(d²/dy²-λ²)²-Re((λU-ω)(d²/dy²-λ²)-λU'')] v = 0

where y is the spatial coordinate, λ the wavenumber (set as the eigenvalue by
default), Re the Reynolds number, ω the oscillation frequency, U the mean flow
velocity profile and v the velocity fluctuation amplitude in the y-direction
(the mode shape).

#Remarks on implementation
This is a Julia port of (our python port of) the Orr-Sommerfeld example from the
Matlab toolbox "NLEVP: A Collection of Nonlinear Eigenvalue Problems"
by  T. Betcke, N. J. Higham, V. Mehrmann, C. Schröder, and F. Tisseur.

The original toolbox is available at :
http://www.maths.manchester.ac.uk/our-research/research-groups/numerical-analysis-and-scientific-computing/numerical-analysis/software/nlevp/

See  [1] or [2] for further refrence.

#References
[1] T. Betcke, N. J. Higham, V. Mehrmann, C. Schröder, and F. Tisseur, NLEVP:
A Collection of Nonlinear Eigenvalue Problems, MIMS EPrint 2011.116, December
2011.

[2] T. Betcke, N. J. Higham, V. Mehrmann, C. Schröder, and F. Tisseur, NLEVP:
A Collection of Nonlinear Eigenvalue Problems. Users' Guide, MIMS EPrint
2010.117, December 2011.
"""
function orr_sommerfeld(N=64, Re=5772, ω =0.26943)
    # Define the Orr-Sommerfeld operator for spatial stability analysis.
    # N: number of Cheb points, R: Reynolds number, w: angular frequency
    N=N+1
    # 2nd- and 4th-order differentiation matrices:
    D,y=cheb(N); D2 =D^2;  D2 = D2[2:N,2:N]
    S= LinearAlgebra.diagm(0=>[0; 1.0./(1. .-y[2:N].^2.); 0])
    D4=(LinearAlgebra.diagm(0=>1. .-y.^2.)*D^4-8*LinearAlgebra.diagm(0=>y)*D^3-12*D^2)*S
    D4 = D4[2:N,2:N]

    I=Array{ComplexF64}(LinearAlgebra.I,N-1,N-1)
    # type conversion
    D2=D2.+0.0im
    D4=D4.+0.0im
    U=LinearAlgebra.diagm(0=>-y[2:N].^2.0.+1)
    # Build Orr-Sommerfeld operator
    L= LinearOperatorFamily(["λ","ω","Re","a"],complex([1.,ω,Re,Inf]))
    push!(L,Term(I,(pow_a(4),),((:λ,),),"λ^4","I"))
    push!(L,Term(1.0im*U,(pow_a(3),pow1,),((:λ,),(:Re,),),"iλ^3Re","i*U"))
    push!(L,Term(-2*D2,(pow2,),((:λ,),),"λ^2","-2D2"))
    push!(L,Term(-1.0im*I,(pow2,pow1,pow1,),((:λ,),(:ω,),(:Re,),),"λ^2*ω*Re","-i*I"))
    push!(L,Term(-1.0im*(U*D2+2.0*I),(pow1,pow1,),((:λ,),(:Re,),),"λ*Re","(U*D2+2*I)"))
    push!(L,Term(1.0im*D2,(pow1,pow1),((:ω,),(:Re,),),"ω*Re","i*D2"))
    push!(L,Term(D4,(),(),"","D4"))
    push!(L,Term(-I,(pow1,),((:a,),),"-a","__aux__"))
    return L,y
end

"""
    L,x,y=biharmonic(N=12;scaleX=2,scaleY=1+sqrt(5))

Discretize the biharmonic equation.

# Arguments
-`N::Int`: Number of collocation points
-`scaleX::Float=2` length of the x-axis
-`scaleY::Float=1+sqrt(5)` length of the y-axis

# Returns
-`L::LinearOperatorFamily`: discretization of the biharmonic operator
- `x:Array`: x-axis
- `y:Array`: y-axis

# Notes
The biharmonic equation models the oscillations of a membrane. It is an an eigenvalue problem that reads:
(∇⁴+εcos(2πx)cos(πy))**v**=λ**v** in Ω
with boundary conditions
**v**=∇²**v**=0

The term εcos(2πx)cos(πy) is included to model some inhomogeneous material properties.
The equation is discretized using Chebyshev collocation.

See also: [`cheb`](@ref)
"""

function biharmonic(N=12;scaleX=2,scaleY=1+sqrt(5))
#N,scaleX,scaleY=12,2,1+sqrt(5)
    N=N+1
    D, xx = cheb(N)
    x = xx/scaleX
    y = xx/scaleY

    #The Chebychev matrices are space dependent scale the accordingly
    Dx = D*scaleX
    Dy = D*scaleY

    D2x=Dx*Dx
    D2y=Dy*Dy
    Dx = Dx[2:N,2:N]
    Dy = Dy[2:N,2:N]
    D2x= D2x[2:N,2:N]
    D2y= D2y[2:N,2:N]
    #Apply BC
    I = Array{ComplexF64}(LinearAlgebra.I,N-1,N-1)
    L = kron(I,D2x) +kron(D2y,I)
    X=kron(ones(N-1),x[2:N])
    Y=kron(y[2:N],ones(N-1))
    P=LinearAlgebra.diagm(0=>cos.(π*2*X).*cos.(π*Y))
    D4= L*L
    I = Array{ComplexF64}(LinearAlgebra.I,(N-1)^2,(N-1)^2)
    L= LinearOperatorFamily(["λ","ε","a"],complex([0,0.,Inf]))
    push!(L,Term(D4,(),(),"","D4"))
    push!(L,Term(P,(pow1,),((:ε,),),"ε","P"))
    push!(L,Term(-I,(pow1,),((:λ,),),"-λ","I"))
    push!(L,Term(-I,(pow1,),((:a,),),"-a","__aux__"))
    return L,x,y
end

## 1DRijkeModel
using SparseArrays
"""
    L,grid=rijke_tube(resolution=128; l=1, c_max=2,mid=0)

Discretize a 1-dimensional Rijke tube. The model is based on the thermoacoustic equations.
This eigenvalue problem reads:
∇c²(x)∇p+ω²p-n exp(-iωτ) ∇p(x_ref)=0 for x in ]0,l[
with boundariy conditions
∇p(0)=p(l)=0
"""
function rijke_tube(resolution=127; l=1, c_max=2,mid=0)
  n=1.0
  tau=2
  #l=1
  c_min=1
  #c_max=2
  outlet=resolution
  outlet_c=c_max
  grid=range(0,stop=l,length=resolution)
  e2p=[(i, i+1) for i =1:resolution-1]
  number_of_elements=resolution-1
  if mid==0
    mid=div(resolution,2)+1# this is the element containing the flame
    #it is not the center if resolotion is odd but then the refernce is in the center
  end
  ref=mid-1 #the reference element is the one in the middle
  flame=(mid)
  #compute element volume
  e2v=diff(grid)
  V=e2v[mid]
  e2c=[i < mid ? c_min : c_max for i = 1:resolution]

  #assemble mass matrix
  m_unit=[2 1;1 2]*1/6
  #preallocation
  ii=Array{Int64}(undef,(resolution-1,2^2))
  jj=Array{Int64}(undef,(resolution-1,2^2))
  mm=Array{ComplexF64}(undef,(resolution-1,2^2))
  for (idx,el) in enumerate(e2p)
    mm[idx,:] = m_unit[:]*e2v[idx]
    ii[idx,:] = [el[1] el[2] el[1] el[2]]
    jj[idx,:] = [el[1] el[1] el[2] el[2]]
  end
  # finally assemble global (sparse) mass matrix
  M =sparse(ii[:],jj[:],mm[:])

  # assemble stiffness matrix
  k_unit = -[1. -1.; -1. 1.]
  #preallocation
  ii=Array{Int64}(undef,(resolution-1,2^2))
  jj=Array{Int64}(undef,(resolution-1,2^2))
  kk=Array{ComplexF64}(undef,(resolution-1,2^2))
  for (idx,el) in enumerate(e2p)
    kk[idx,:]=k_unit[:]/e2v[idx]*e2c[idx]^2
    ii[idx,:] = [el[1] el[2] el[1] el[2]]
    jj[idx,:] = [el[1] el[1] el[2] el[2]]
  end
  # finally assemble global (sparse) stiffness matrix
  K = sparse(ii[:],jj[:],kk[:])

  #assemble boundary mass matrix
  bb=[-outlet_c*1im]
  ii=[outlet]
  jj=[outlet]
  B = sparse(ii,jj,bb)
  #assemble flame matrix
  #preallocation
  ii=Array{Int64}(undef,(length(flame),2*2))
  jj=Array{Int64}(undef,(length(flame),2*2))
  qq=Array{ComplexF64}(undef,(length(flame),2*2))

  grad_p_ref=[-1 1]
  grad_p_ref=grad_p_ref/e2v[ref]*1
  #println("##!!!!!Here!!!!!##'")
  #println("flame:",flame)
  #println("ref:",ref)
  for (idx,el) in enumerate(flame)
    jj[idx,:]=[e2p[ref][1] e2p[ref][2] e2p[ref][1] e2p[ref][2]]
    ii[idx,:]=[e2p[el][1] e2p[el][1] e2p[el][2] e2p[el][2]]
    run=1
    for i=1:2
      for j =1:2
        qq[idx,run]= grad_p_ref[1,j]*e2v[el]/2
        run+=1
      end
    end
  end
  #finally assemble flame matrix
  Q= sparse(ii[:],jj[:],-qq[:],resolution,resolution)
  Q=Q/V
  L=LinearOperatorFamily(["ω","n","τ","Y","λ"],complex([0.,n,tau,1E15,Inf]))
  push!(L,Term(M,(pow2,),((:ω,),),"ω^2","M"))
  push!(L,Term(K,(),(),"","K"))
  push!(L,Term(B,(pow1,pow1,),((:ω,),(:Y,),),"ω*Y","C"))
  push!(L,Term(Q,(pow1,exp_delay,),((:n,),(:ω,:τ,)),"n*exp(-i ω τ)","Q"))
  push!(L,Term(-M,(pow1,),((:λ,),),"-λ","__aux__"))
  #push!(L,Term(-SparseArrays.I,(pow1,),((:λ,),),"-λ","__aux__"))
 return L,grid
end

end


text:
using Gen
using Distributions
using LinearAlgebra
include("structs.jl")
include("distributions.jl")

@gen function gen_scene_and_observations(
    max_objects::Int64,
    num_object_categories::Int64,
    num_observations::Int64,
    room::RoomParams,
    camera::CameraParams,
    detection_sd::Float64
    )

    # Generates scene
    num_objects = @trace(uniform_discrete(1, max_objects), :num_objects)
    println("num_objects: ", num_objects)
    categories = [@trace(uniform_discrete(1, num_object_categories),
        :objects=>i=>:category) for i in 1:num_objects]
    locations = [@trace(mvuniform(room.bottom_left_lim, room.top_right_lim),
        :objects=>i=>:location) for i in 1:num_objects]
    # Generates observations
    observations = []
    for i in 1:num_observations
        camera_position = @trace(mvuniform(room.bottom_left_lim, room.top_right_lim),
            :observations=>i=>:camera_location)

        object = @trace(uniform_discrete(1, num_objects), :observations=>i=>:object)

        category_dist = [k == object ? 1 : 0 for k in 1:num_objects]

        category = @trace(categorical(category_dist), :observations=>i=>:category)

        location = locations[object]

        direction = @trace(choose_direction(camera_position, location, detection_sd), :observations=>i=>:direction)

        append!(observations, Dict(:camera_location => camera_position,
            :category => category,
            :direction => direction,
            :object => object))
    end
    return num_objects, categories, locations, observations
end

args = (5, 10, 20, RoomParams(), CameraParams(), 0.5)

(trace, _) = Gen.generate(gen_scene_and_observations, args)

Gen.get_choices(trace)

max_stars_repo_name: zhangir-azerbayev/MetaGen-New

max_issues_repo_name: zhangir-azerbayev/MetaGen-New

max_forks_repo_name: zhangir-azerbayev/MetaGen-New


text:
using Test
@testset "variational_circuit" begin
    include("variational_circuit.jl")
end

@testset "hamiltonians" begin
    include("hamiltonians.jl")
end

@testset "supremacy_circuit" begin
    include("supremacy_circuit.jl")
end

@testset "google53" begin
    include("google53.jl")
end

@testset "general_U4" begin
    include("general_U4.jl")
end


text:
struct MyStruct
    a::Int
end

@njl_getprop MyStruct

function Base.getproperty(o::MyStruct, prop::Val{:b})
    o.a * 2
end

@testset "Val dispatch for properties" begin
    ms = MyStruct(3)
    @test ms.a == 3
    @test ms.b == 6
    @test Set(propertynames(ms)) == Set([:a, :b])
end


text:
# interface to the PETSc options database(s), by providing an OPTIONS[T]
# dictionary-like object that is analogous to the Julia Base.ENV object
# for environment variables.

mutable struct Options{T<:Scalar} <: AbstractDict{String,String}; end
const OPTIONS = Dict(T => Options{T}() for T in C.petsc_type)
export OPTIONS, withoptions

const SymOrStr = Union{AbstractString,Symbol}

function Base.setindex!(::Options{T}, v, k::SymOrStr) where {T}
  chk(C.PetscOptionsSetValue(T, string('-',k), string(v)))
  return v
end

function Base.setindex!(::Options{T}, v::Cvoid, k::SymOrStr) where {T}
  chk(C.PetscOptionsClearValue(T, string('-',k)))
  return v
end

# PETSc complains if you don't use an option, remove this?
# allow OPTIONS[k]=v to set options for all PETSc scalar types simultaneously
function Base.setindex!(o::typeof(OPTIONS), v, k::SymOrStr)
  for opts in values(o)
    opts[k] = v
  end
  return v
end

const _optionstr = fill(UInt8(0),1024)
function Base.get(::Options{T}, k::SymOrStr, def) where {T}
  b = Ref{PetscBool}()
  chk(C.PetscOptionsGetString(T, Cstring(Ptr{UInt8}(C_NULL)), string('-',k),
                              pointer(_optionstr), Csize_t(length(_optionstr)),
                              b))
  return b[] != 0 ? unsafe_string(pointer(_optionstr)) : def
end

function Base.haskey(::Options{T}, k::SymOrStr) where {T}
  b = Ref{PetscBool}()
  chk(C.PetscOptionsHasName(T, Cstring(Ptr{UInt8}(C_NULL)), string('-',k), b))
  return b[] != 0
end

Base.similar(::Options) = Dict{String,String}()

Base.pop!(o::Options, k::SymOrStr) = (v = o[k]; o[k] = nothing; v)
Base.pop!(o::Options, k::SymOrStr, def) = haskey(o,k) ? pop!(o,k) : def
Base.delete!(o::Options, k::SymOrStr) = (o[k] = nothing; o)
Base.delete!(o::Options, k::SymOrStr, def) = haskey(o,k) ? delete!(o,k) : def

# need to override show: default show function doesn't work because
# there seems to be no way to iterate over the PETSc options database (grr).
Base.show(io::IO, ::Options{T})  where {T} = print(io, "PETSc{$T} options database")

# temporarily set some options, call f, and then unset them; like withenv
function withoptions(f::Function, ::Type{T}, keyvals) where {T<:Scalar}
  old = Dict{SymOrStr,Any}()
  o = OPTIONS[T]
  for (key,val) in keyvals
    old[key] = get(o,key,nothing)
    val !== nothing ? (o[key]=val) : delete!(o, key)
  end
  try f()
  finally
    for (key,val) in old
      val !== nothing ? (o[key]=val) : delete!(o, key)
    end
  end
end
withoptions(f::Function, ::Type{T}, keyvals::Pair{K}...) where {T<:Scalar,K<:SymOrStr} = withoptions(f, T, keyvals)
withoptions(f::Function, ::Type{T}) where {T<:Scalar} = f() # handle empty keyvals case; see julia#10853


text:
using Distributions, Plots; pyplot()

n, N = 30, 10^6

dist1 = Uniform(1-sqrt(3),1+sqrt(3))
dist2 = Exponential(1)
dist3 = Normal(1,1)

data1 = [mean(rand(dist1,n)) for _ in 1:N]
data2 = [mean(rand(dist2,n)) for _ in 1:N]
data3 = [mean(rand(dist3,n)) for _ in 1:N]

stephist([data1 data2 data3], bins=100,
    c=[:blue :red :green], xlabel = "x", ylabel = "Density",
    label=["Average of Uniforms" "Average of Exponentials" "Average of Normals"],
    normed=true, xlims=(0,2), ylims=(0,2.5))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module ResourceBundles

export ResourceBundles
export ResourceBundle, resource_bundle, @resource_bundle
export LocaleId, LC
export @tr_str, string_to_key

include("types.jl")
include("constants.jl")
include("locale_iso_data.jl")
include("resource_bundle.jl")
include("localetrans.jl")
include("locale.jl")

include("string.jl")
include("poreader.jl")

include("clocale.jl")

end # module ResourceBundles


text:
using Yao
using BitBasis
using Random
using Test, LinearAlgebra
using OrdinaryDiffEq
using QuDiffEq

function diffeq_problem(nbit::Int)
    siz = 1<<nbit
    A = (rand(ComplexF64, siz,siz))
    A = (A + A')/2
    b = normalize!(rand(ComplexF64, siz))
    x = normalize!(rand(ComplexF64, siz))
    A, b, x
end

@testset "Linear_differential_equations_HHL" begin
    Random.seed!(2)
    N = 1
    h = 0.1
    tspan = (0.0,0.6)
    N_t = round(Int, 2*(tspan[2] - tspan[1])/h + 3)
    M, v, x = diffeq_problem(N)
    A(t) = M
    b(t) = v
    nreg = 12
    f(u,p,t) = M*u + v;
    prob = ODEProblem(f, x, tspan)
    qprob = QuLDEProblem(A, b, x, tspan)
    sol = solve(prob, Tsit5(), dt = h, adaptive = :false)
    s = vcat(sol.u...)

    res = solve(qprob, QuEuler(nreg), dt = h)
    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
    @test isapprox.(s, r, atol = 0.5) |> all

    res = solve(qprob, QuLeapfrog(nreg), dt = h)
    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
    @test isapprox.(s, r, atol = 0.3) |> all

    res = solve(qprob, QuAB2(nreg), dt = h)
    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
    @test isapprox.(s, r, atol = 0.3) |> all

    res = solve(qprob, QuAB3(nreg), dt = h)
    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
    @test isapprox.(s, r, atol = 0.3) |> all

    res = solve(qprob, QuAB4(nreg),dt = h)
    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
    @test isapprox.(s, r, atol = 0.3) |> all
end;