text:
function _generate_2d_nodes!(nodes, nx, ny, LL, LR, UR, UL)
      for i in 0:ny-1
        ratio_bounds = i / (ny-1)

        x0 = LL[1] * (1 - ratio_bounds) + ratio_bounds * UL[1]
        x1 = LR[1] * (1 - ratio_bounds) + ratio_bounds * UR[1]

        y0 = LL[2] * (1 - ratio_bounds) + ratio_bounds * UL[2]
        y1 = LR[2] * (1 - ratio_bounds) + ratio_bounds * UR[2]

        for j in 0:nx-1
            ratio = j / (nx-1)
            x = x0 * (1 - ratio) + ratio * x1
            y = y0 * (1 - ratio) + ratio * y1
            push!(nodes, Node(Tensors.Vec{2}((x, y))))
        end
    end
end

# Check edge orientation consistency
function _check_node_data(nodes, n1,n2,n3)
    a = nodes[n2].x-nodes[n1].x
    b = nodes[n3].x-nodes[n1].x
    if (a[1]*b[2]-a[2]*b[1]) < 0
        #swap nodes 2 and 3
        return (n1,n3,n2)
    end
    return (n1,n2,n3)
end

@inline _mapToGlobalIdx(cells,cellidx,localnodeidx) = cells[cellidx].nodes[localnodeidx]

function _get_nodeset_from_edges(cells,edgeset,CellType)
    nodes = Set{Int}()
    for edge in edgeset
        push!(nodes, _mapToGlobalIdx(cells, edge.cellidx, reference_edge_nodes(CellType)[edge.idx][1]))
        push!(nodes, _mapToGlobalIdx(cells, edge.cellidx, reference_edge_nodes(CellType)[edge.idx][2]))
    end
    return nodes
end

#########################
# Line 1D
###########################
# Line
line_mesh(nel::Int, left::Number=-1.0, right::Number=1.0) =
    line_mesh(LineCell, (nel,), Tensors.Vec{1}((left,)), Tensors.Vec{1}((right,)))
"""
    line_mesh(celltype::Cell, nel::NTuple, [left::Vec, right::Vec)
Return a `mesh` for a line in 1 dimension. `celltype` defined the type of cell,
e.g. `LineCell`. `nel` is a tuple of the number of elements in each direction.
`left` and `right` are optional endpoints of the domain. Defaults to -1 and 1 in all directions.
"""
function line_mesh(::Type{LineCell}, nel::NTuple{1,Int}, left::Tensors.Vec{1,T}=Tensors.Vec{1}((-1.0,)), right::Tensors.Vec{1,T}=Tensors.Vec{1}((1.0,))) where {T}
    nel_x = nel[1]
    n_nodes = nel_x + 1

    # Generate nodes
    coords_x = collect(range(left[1], stop=right[1], length=n_nodes))
    nodes = Node{1,T}[]
    for i in 1:n_nodes
        push!(nodes, Node(Tensors.Vec{1}((coords_x[i],))))
    end

    # Generate cells
    cells = LineCell[]
    for i in 1:nel_x
        push!(cells, LineCell((i, i+1)))
    end

    # Cell node sets
    nodesets = Dict("left"  => Set{Int}([1]),
                    "right" => Set{Int}([nel_x+1]))
    return PolytopalMesh(cells, nodes; nodesets = nodesets)
end

#########################
# Triangle Cells 2D   #
#########################
"""
rectangle_mesh(::Type{RefTetrahedron}, ::Type{Val{2}}, nel::NTuple{2,Int}, LL::Vec{2,T}, UR::Vec{2,T})
Generate a rectangular mesh with triangular cells, where `LL` is the low left node
and `UR` is the upper right one. `nel` is a tuple with the number of partions to be
used in each dimension.
"""
function rectangle_mesh(::Type{TriangleCell}, nel::NTuple{2,Int}, LL::Tensors.Vec{2,T}, UR::Tensors.Vec{2,T}) where {T}
    LR = Tensors.Vec{2}((UR[1],LL[2]))
    UL = Tensors.Vec{2}((LL[1],UR[2]))
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = 2*nel_x*nel_y
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = TriangleCell[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, TriangleCell((node_array[i,j], node_array[i+1,j], node_array[i,j+1]))) # ◺
        push!(cells, TriangleCell((node_array[i+1,j], node_array[i+1,j+1], node_array[i,j+1]))) # ◹
    end

    # Cell edges
    cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
    boundary = EdgeIndex[[EdgeIndex(cl, 3) for cl in cell_array[1,:,1]];
                           [EdgeIndex(cl, 3) for cl in cell_array[2,end,:]];
                           [EdgeIndex(cl, 1) for cl in cell_array[2,:,end]];
                           [EdgeIndex(cl, 2) for cl in cell_array[1,1,:]]]
    offset = 0
    edgesets = Dict{String,Set{EdgeIndex}}()
    edgesets["bottom"] = Set{EdgeIndex}(boundary[(1:length(cell_array[1,:,1]))   .+ offset]); offset += length(cell_array[1,:,1])
    edgesets["right"]  = Set{EdgeIndex}(boundary[(1:length(cell_array[2,end,:])) .+ offset]); offset += length(cell_array[2,end,:])
    edgesets["top"]    = Set{EdgeIndex}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
    edgesets["left"]   = Set{EdgeIndex}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
    edgesets["boundary"] = union(edgesets["bottom"],edgesets["right"],edgesets["top"],edgesets["left"])
    nodesets = Dict{String,Set{Int}}()
    for set in edgesets
        nodesets[set.first] = _get_nodeset_from_edges(cells,set.second, TriangleCell)
    end
    return PolytopalMesh(cells, nodes; edgesets = edgesets, nodesets = nodesets)
end

#########################
# Rectangle Cells 2D   #
#########################
@inline _build_cell(::Type{RectangleCell}, el_nodes, el_faces) = RectangleCell(el_nodes,(el_faces[1],el_faces[2],el_faces[3],el_faces[4]))
"""
rectangle_mesh(::Type{RectangleCell}, nel::NTuple{2,Int}, LL::Vec{2,T}, UR::Vec{2,T}) where {T}
Generate a rectangular mesh with triangular cells, where `LL` is the low left node
and `UR` is the upper right one. `nel` is a tuple with the number of partions to be
used in each dimension.
"""
function rectangle_mesh(::Type{RectangleCell}, nel::NTuple{2,Int}, LL::Tensors.Vec{2,T}, UR::Tensors.Vec{2,T}) where {T}
    LR = Tensors.Vec{2}((UR[1],LL[2]))
    UL = Tensors.Vec{2}((LL[1],UR[2]))
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = nel_x*nel_y
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = RectangleCell[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, RectangleCell((node_array[i,j], node_array[i+1,j], node_array[i+1,j+1], node_array[i,j+1])))
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
    boundary = EdgeIndex[[EdgeIndex(cl, 1) for cl in cell_array[:,1]];
                      [EdgeIndex(cl, 2) for cl in cell_array[end,:]];
                      [EdgeIndex(cl, 3) for cl in cell_array[:,end]];
                      [EdgeIndex(cl, 4) for cl in cell_array[1,:]]]
    # Cell face sets
    offset = 0
    edgesets = Dict{String, Set{EdgeIndex}}()
    edgesets["bottom"] = Set{EdgeIndex}(boundary[(1:length(cell_array[:,1]))   .+ offset]); offset += length(cell_array[:,1])
    edgesets["right"]  = Set{EdgeIndex}(boundary[(1:length(cell_array[end,:])) .+ offset]); offset += length(cell_array[end,:])
    edgesets["top"]    = Set{EdgeIndex}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
    edgesets["left"]   = Set{EdgeIndex}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
    edgesets["boundary"] = union(edgesets["bottom"],edgesets["right"],edgesets["top"],edgesets["left"])
    nodesets = Dict{String,Set{Int}}()
    for set in edgesets
        nodesets[set.first] = _get_nodeset_from_edges(cells,set.second, RectangleCell)
    end
    return PolytopalMesh(cells, nodes; edgesets = edgesets, nodesets = nodesets)
end


text:
if Base.HOME_PROJECT[] !== nothing
    # JuliaLang/julia/pull/28625
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, Optim

# use include("Rosenbrock.jl") etc
# Generate examples
include("generate.jl")

cp(joinpath(@__DIR__, "..", "LICENSE.md"),
   joinpath(@__DIR__, "src", "LICENSE.md"); force = true)

#run('mv ../CONTRIBUTING.md ./dev/CONTRIBUTING.md') # TODO: Should we use the $odir/CONTRIBUTING.md file instead?
makedocs(
    doctest = false,
    sitename = "Optim"
)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-windmill"),
    repo = "github.com/JuliaNLSolvers/Optim.jl.git",
    julia = "1.0"
)

ext: jl


text:
"""
    slope(spl::Spline,dλ)

Returns the pair of `λ, s'` (1st derivative of the spline function).
"""
function slope(spl::Spline,dλ)
    xs=collect(dλ)#,(x[end]-x[1])/dλ)  # spline x vector
    s=zeros(length(xs))         # spline spectral vector
    j=1
    @inbounds for i in 1 : length(xs)
        if j < length(spl.x)
            while xs[i] > spl.x[j+1]
                j+=1
            end
            delta = xs[i]-spl.x[j]
            if typeof(spl) == CubicSpline
                s[i]= 3.0 * spl.d[j] * delta * delta + 2.0 * spl.c[j] * delta + spl.b[j]
            elseif typeof(spl) == HermiteSpline
                xb = (xs[i] - spl.x[j]) / (spl.x[j+1] - spl.x[j])
                xb2 = xb * xb
                #Φ0_ = 6 * xb^2 - 6 * xb
                #Φ1_ = 3 * xb^2 - 4 * xb + 1
                #Ψ0_ = -6 * xb^2 + 6 * xb
                #Ψ1_ = 3 * xb^2 - 2 * xb
                h = spl.x[j+1] - spl.x[j]
                #s[i] = (spl.a[j] * Φ0_ + h * spl.b[j] * Φ1_ + spl.a[j+1] * Ψ0_ + h * spl.b[j+1] * Ψ1_) / h
                s[i] = (spl.a[j] * (6 * xb2 - 6 * xb) + h * spl.b[j] * (3 * xb2 - 4 * xb + 1) + spl.a[j+1] * (-6 * xb2 + 6 * xb) + h * spl.b[j+1] * (3 * xb2 - 2 * xb)) / h
            elseif typeof(spl) == QuadraticSpline
                s[i]= 2.0 * spl.c[j] * delta + spl.b[j]
            elseif typeof(spl) == LinearSpline
                s[i] = spl.b[j]
            elseif typeof(spl) == NearestNeighborSpline
                s[i] = 0.0
            end
        end
    end
    return xs,s
end

# alias
const deriv1 = slope

"""
    curvature(spl::Spline,dλ)

Returns the 2nd derivative of the spline function (s'')
"""
function curvature(spl::Spline,dλ)
    xs=collect(dλ)#,(x[end]-x[1])/dλ)  # spline x vector
    s=zeros(length(xs))         # spline spectral vector
    j=1
    @inbounds for i in 1 : length(xs)
        if j < length(spl.x)
            while xs[i] > spl.x[j+1]
                j+=1
            end
            delta = xs[i]-spl.x[j]
            if typeof(spl) == CubicSpline
                s[i] = 6.0 * spl.d[j] * delta + 2.0 * spl.c[j]
            elseif typeof(spl) == HermiteSpline
                xb = (xs[i] - spl.x[j]) / (spl.x[j+1] - spl.x[j])
                #Φ0__ = 12 * xb - 6
                #Φ1__ = 6 * xb - 4
                #Ψ0__ = -12 * xb + 6
                #Ψ1__ = 6 * xb - 2
                h = spl.x[j+1] - spl.x[j]
                #s[i] = (spl.a[j] * Φ0__ + h * spl.b[j] * Φ1__ + spl.a[j+1] * Ψ0__ + h * spl.b[j+1] * Ψ1__) / h
                s[i] = (spl.a[j] * (12 * xb - 6) + h * spl.b[j] * (6 * xb - 4) + spl.a[j+1] * (-12 * xb + 6) + h * spl.b[j+1] * (6 * xb - 2)) / h
                #s[i] = 2.0 * spl.b[j]
            elseif typeof(spl) == QuadraticSpline
                s[i] = 2.0 * spl.c[j]
            elseif typeof(spl) == LinearSpline
                s[i] = 0.0
            elseif typeof(spl) == NearestNeighborSpline
                s[i] = 0.0
            else
                nothing
            end
        end
    end
    return xs,s
end

# alias
const deriv2 = curvature

"""
    curvrate(spl::Spline,dλ)

Returns the 3rd derivative of the spline function (s''').
"""
function curvrate(spl::Spline,dλ)
    xs=collect(dλ)#,(x[end]-x[1])/dλ)  # spline x vector
    s=zeros(length(xs))         # spline spectral vector
    j=1
    @inbounds for i in 1 : length(xs)
        if j < length(spl.x)
            while xs[i] > spl.x[j+1]
                j+=1
            end
            typeof(spl) == CubicSpline ? s[i]= 6.0 * spl.d[j] :
            typeof(spl) == QuadraticSpline ? s[i] = 0.0 :
            typeof(spl) == LinearSpline ? s[i] = 0.0 :
            typeof(spl) == NearestNeighborSpline ? s[i] = 0.0 : nothing
        end
    end
    return xs,s
end

# alias
const deriv3 = curvrate


text:
module JLBOX_MODULE

end # module JLBOX_MODULE

max_stars_repo_name: jbnicolai/jlbox

max_issues_repo_name: jbnicolai/jlbox

max_forks_repo_name: jbnicolai/jlbox


text:
#!/usr/bin/env julia

# For the files given in ARGS, create all the needed continuation jobs
# for the second phase. That includes all the accessible T2 and all
# the tw₁ starts.

const T2iters = 1e8
const availableT2 = [0.9, 0.8, 0.7, 0.6, 0.5]
const tw1list = map(Int64,[1e4,1e5,1e6,1e7,1e8])

const param = "-m 33 -t 1e4 -T 1e8 -n 10 -L"

function get_fields(file)
    regex = r"quench_T_([0-9.]+)_L_([0-9]+)_([a-z]+)_R([AB]).net"
    fields = match(regex, basename(file)).captures
    length(fields) == 0 && error("Couldn't get metadata from $file")
    T = parse(Float64, fields[1])
    L = parse(Int64, fields[2])
    id = fields[3] * "-R" * fields[4]
    return T, L, id
end

# Return strings containing the continuation jobs for .net file
function get_job_list(file)
    T1, L, id = get_fields(file)
    jobs = String[]
    T2list = [t for t ∈ availableT2 if t<T1]
    for T2 in T2list
        for tw1 in tw1list
            outname = "continue_at_T2_$(T2)_tw1_$(tw1)_for_T1_$(T1)_L_$(L)_id_$(id).net"
            cmd = "./exe/annealer -c $tw1  $param -p -l $L -i $T2iters $T2"
            job = "cat $file | $cmd \"$(outname)\""
            cleanup = "lzma $(outname)"
            push!(jobs, job * "\n" * cleanup * "\n" * "echo 'All done! ($outname)'")
        end
    end
    return jobs
end

jobs = String[]
for file in ARGS
    for j in get_job_list(file)
        push!(jobs,j)
    end
end
info("$(length(jobs)) jobs processed.")

n = 1
for job in jobs
    header=
    """
    #!/bin/sh
    #\$ -o \$HOME/glass_continue_$(n).out -j y
    #\$ -N glass_$(n)

    cd tfg/

    """
    open("continueT2_$n.sh", "w") do io
        println(io,header)
        println(io,job)
    end
    n+=1
end

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


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

@doc raw"""APIServiceSpec contains information for locating and communicating with a server. Only https is supported, though you are able to disable certificate verification.

    IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec(;
        caBundle=nothing,
        group=nothing,
        groupPriorityMinimum=nothing,
        insecureSkipTLSVerify=nothing,
        service=nothing,
        version=nothing,
        versionPriority=nothing,
    )

    - caBundle::Vector{UInt8} : CABundle is a PEM encoded CA bundle which will be used to validate an API server's serving certificate. If unspecified, system trust roots on the apiserver are used.
    - group::String : Group is the API group name this server hosts
    - groupPriorityMinimum::Int32 : GroupPriorityMininum is the priority this group should have at least. Higher priority means that the group is preferred by clients over lower priority ones. Note that other versions of this group might specify even higher GroupPriorityMininum values such that the whole group gets a higher priority. The primary sort is based on GroupPriorityMinimum, ordered highest number to lowest (20 before 10). The secondary sort is based on the alphabetical comparison of the name of the object.  (v1.bar before v1.foo) We'd recommend something like: *.k8s.io (except extensions) at 18000 and PaaSes (OpenShift, Deis) are recommended to be in the 2000s
    - insecureSkipTLSVerify::Bool : InsecureSkipTLSVerify disables TLS certificate verification when communicating with this server. This is strongly discouraged.  You should use the CABundle instead.
    - service::IoK8sKubeAggregatorPkgApisApiregistrationV1beta1ServiceReference : Service is a reference to the service for this API server.  It must communicate on port 443 If the Service is nil, that means the handling for the API groupversion is handled locally on this server. The call will simply delegate to the normal handler chain to be fulfilled.
    - version::String : Version is the API version this server hosts.  For example, \"v1\"
    - versionPriority::Int32 : VersionPriority controls the ordering of this API version inside of its group.  Must be greater than zero. The primary sort is based on VersionPriority, ordered highest to lowest (20 before 10). Since it's inside of a group, the number can be small, probably in the 10s. In case of equal version priorities, the version string will be used to compute the order inside a group. If the version string is \"kube-like\", it will sort above non \"kube-like\" version strings, which are ordered lexicographically. \"Kube-like\" versions start with a \"v\", then are followed by a number (the major version), then optionally the string \"alpha\" or \"beta\" and another number (the minor version). These are sorted first by GA > beta > alpha (where GA is a version with no suffix such as beta or alpha), and then by comparing major version, then minor version. An example sorted list of versions: v10, v2, v1, v11beta2, v10beta3, v3beta1, v12alpha1, v11alpha2, foo1, foo10.
"""
mutable struct IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec <: SwaggerModel
    caBundle::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: caBundle
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    groupPriorityMinimum::Any # spec type: Union{ Nothing, Int32 } # spec name: groupPriorityMinimum
    insecureSkipTLSVerify::Any # spec type: Union{ Nothing, Bool } # spec name: insecureSkipTLSVerify
    service::Any # spec type: Union{ Nothing, IoK8sKubeAggregatorPkgApisApiregistrationV1beta1ServiceReference } # spec name: service
    version::Any # spec type: Union{ Nothing, String } # spec name: version
    versionPriority::Any # spec type: Union{ Nothing, Int32 } # spec name: versionPriority

    function IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec(;caBundle=nothing, group=nothing, groupPriorityMinimum=nothing, insecureSkipTLSVerify=nothing, service=nothing, version=nothing, versionPriority=nothing)
        o = new()
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("caBundle"), caBundle)
        setfield!(o, Symbol("caBundle"), caBundle)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("groupPriorityMinimum"), groupPriorityMinimum)
        setfield!(o, Symbol("groupPriorityMinimum"), groupPriorityMinimum)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("insecureSkipTLSVerify"), insecureSkipTLSVerify)
        setfield!(o, Symbol("insecureSkipTLSVerify"), insecureSkipTLSVerify)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("service"), service)
        setfield!(o, Symbol("service"), service)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec, Symbol("versionPriority"), versionPriority)
        setfield!(o, Symbol("versionPriority"), versionPriority)
        o
    end
end # type IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec

const _property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec = Dict{Symbol,Symbol}(Symbol("caBundle")=>Symbol("caBundle"), Symbol("group")=>Symbol("group"), Symbol("groupPriorityMinimum")=>Symbol("groupPriorityMinimum"), Symbol("insecureSkipTLSVerify")=>Symbol("insecureSkipTLSVerify"), Symbol("service")=>Symbol("service"), Symbol("version")=>Symbol("version"), Symbol("versionPriority")=>Symbol("versionPriority"))
const _property_types_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec = Dict{Symbol,String}(Symbol("caBundle")=>"Vector{UInt8}", Symbol("group")=>"String", Symbol("groupPriorityMinimum")=>"Int32", Symbol("insecureSkipTLSVerify")=>"Bool", Symbol("service")=>"IoK8sKubeAggregatorPkgApisApiregistrationV1beta1ServiceReference", Symbol("version")=>"String", Symbol("versionPriority")=>"Int32")
Base.propertynames(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }) = collect(keys(_property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec))
Swagger.property_type(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec[name]))}
Swagger.field_name(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }, property_name::Symbol) =  _property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec[property_name]

function check_required(o::IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec)
    (getproperty(o, Symbol("groupPriorityMinimum")) === nothing) && (return false)
    (getproperty(o, Symbol("service")) === nothing) && (return false)
    (getproperty(o, Symbol("versionPriority")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1beta1APIServiceSpec }, name::Symbol, val)
    if name === Symbol("caBundle")
    end
end


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

mutable struct IoK8sApiCoreV1PodAntiAffinity <: SwaggerModel
    preferredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1WeightedPodAffinityTerm} } # spec name: preferredDuringSchedulingIgnoredDuringExecution
    requiredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodAffinityTerm} } # spec name: requiredDuringSchedulingIgnoredDuringExecution

    function IoK8sApiCoreV1PodAntiAffinity(;preferredDuringSchedulingIgnoredDuringExecution=nothing, requiredDuringSchedulingIgnoredDuringExecution=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodAntiAffinity, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        validate_property(IoK8sApiCoreV1PodAntiAffinity, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        o
    end
end # type IoK8sApiCoreV1PodAntiAffinity

const _property_map_IoK8sApiCoreV1PodAntiAffinity = Dict{Symbol,Symbol}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>Symbol("preferredDuringSchedulingIgnoredDuringExecution"), Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>Symbol("requiredDuringSchedulingIgnoredDuringExecution"))
const _property_types_IoK8sApiCoreV1PodAntiAffinity = Dict{Symbol,String}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>"Vector{IoK8sApiCoreV1WeightedPodAffinityTerm}", Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>"Vector{IoK8sApiCoreV1PodAffinityTerm}")
Base.propertynames(::Type{ IoK8sApiCoreV1PodAntiAffinity }) = collect(keys(_property_map_IoK8sApiCoreV1PodAntiAffinity))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodAntiAffinity }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1PodAntiAffinity[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodAntiAffinity }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodAntiAffinity[property_name]

function check_required(o::IoK8sApiCoreV1PodAntiAffinity)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodAntiAffinity }, name::Symbol, val)
end


text:
tests = ["parse", "create", "cdata"]

for t in tests
    fpath = "$t.jl"
    @printf("running %s ...\n", fpath)
    include(fpath)
end


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

import Benchmarks

using HypothesisTests

noToken = [lowDiversityWithChallenge(1000, 50, 70, Benchmarks.registryMean) for i in 1:20]
v = var(noToken)
m = mean(noToken)

println("Var: $v Mean: $m")

token = [binaryTokenChallenge(1000, 50, 70, Benchmarks.registryMean) for i in 1:20]
vt = var(token)
mt = mean(token)

println("Var: $vt Mean: $mt")

t = UnequalVarianceTTest(noToken, token)

println(t)

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
"""
minmod(a::Float64,b::Float64,c::Float64)

If the signs of `a, b, c` are the same return the smallest (closest to 0.0),
otherwise, return 0.0.
"""
function minmod(a::Float64,b::Float64,c::Float64)
    s = (Int.(sign(a))+Int.(sign(b))+Int.(sign(c)))÷3

    if (abs(s)===1)
        return sign(a)*min(abs(a),abs(b),abs(c))
    else
        return 0.0
    end
end

function minmodB(a::Float64,b::Float64,c::Float64,M,h)
    b = b+M*h^2*sign(b)
    c = c+M*h^2*sign(c)

    return minmod(a,b,c)
end

"""
linear(cell_coeffs::Vector{Float64}, V::Matrix{Float64}, Vinv::Matrix{Float64},
w::Float64, Δ::Float64)

Take coefficients legendre basis polynomials and project them down to a linear polynomial.
"""
function linear(cell_coeffs::Vector{Float64}, V::Matrix{Float64}, Vinv::Matrix{Float64},
    D::Matrix{Float64}, w::Vector{Float64}, Δ::Float64)

    # reweight to canonical legendre basis
    # convert to lagrange basis
    cell_coeffs = Vinv*(cell_coeffs*2.0./(Δ*w))
    # project to order 1
    if length(cell_coeffs)>2
        cell_coeffs[3:end] .= 0.0
    end
    # while we're here, also get the slope of the linear function
    slope = (D*cell_coeffs)[1]*2.0/Δ
    # convert back to original basis
    cell_coeffs = (V*cell_coeffs)./(2.0./(Δ*w))

    return cell_coeffs, slope
end

function limit(coeffs::Vector{Float64}, V::Matrix{Float64}, Vinv::Matrix{Float64},
    D::Matrix{Float64}, w::Vector{Float64}, Δvec::AbstractVector{Float64},
    cell_nodes_matrix::Matrix{Float64}, num_phases::Int, n₋::Int, n₊::Int)

    # V, Vinv, D, w = vandermonde(n_bases_per_cell(d.dq))

    limited_coeffs = copy(coeffs)

    limited_coeffs = limited_coeffs[n₋+1:end-n₊]

    limited_coeffs = reshape(limited_coeffs,length(w),num_phases,length(Δvec))

    cell_averages = sum(limited_coeffs,dims=1)
    cell_averages = reshape(cell_averages,num_phases,length(Δvec))./transpose(Δvec)

    poly_reconstruct_left = limited_coeffs[1,:,:]*2.0./(Δvec[:]'*w[1])
    poly_reconstruct_right = limited_coeffs[end,:,:]*2.0./(Δvec[:]'*w[end])

    for phase in 1:num_phases
        for cell in 1:length(Δvec)
            a = cell_averages[phase,cell] - poly_reconstruct_left[phase,cell]
            b = cell_averages[phase,cell] - cell_averages[phase,max(1,cell-1)]
            c = cell_averages[phase,min(length(Δvec),cell+1)] - cell_averages[phase,cell]

            left_limited_flux = cell_averages[phase,cell] - minmod(a,b,c)

            a = poly_reconstruct_right[phase,cell] - cell_averages[phase,cell]

            right_limited_flux = cell_averages[phase,cell] + minmod(a,b,c)

            needs_limiting = !(isapprox(left_limited_flux,poly_reconstruct_left[phase,cell];atol=1e-8) &&
                isapprox(right_limited_flux,poly_reconstruct_right[phase,cell];atol=1e-8))
            if needs_limiting
                cell_coeffs = limited_coeffs[:,phase,cell]
                cell_coeffs, slope = linear(cell_coeffs,V,Vinv,D,w,Δvec[cell])

                nodes = cell_nodes_matrix[:,cell] .- (cell_nodes_matrix[1,cell]+cell_nodes_matrix[end,cell])/2.0

                limited_coeffs[:,phase,cell] = (cell_averages[phase,cell] .+
                    (nodes*minmod(slope,b/Δvec[cell],c/Δvec[cell])))./(2.0./(Δvec[cell]*w))
            end
        end
    end
    return [coeffs[1:n₋]; limited_coeffs[:]; coeffs[end-n₊+1:end]]
end

function limit(d::SFMDistribution{DGMesh{T}}) where T
    V, Vinv, D, w = vandermonde(n_bases_per_cell(d.dq))

    limited_coeffs = limit(d.coeffs[:],V,Vinv,D,w,Δ(d.dq),cell_nodes(d.dq),n_phases(d.dq),N₋(d.dq),N₊(d.dq))
    return SFMDistribution(limited_coeffs,d.dq)
end

# abstract type AbstractLimiter end

struct Limiter #<: AbstractLimiter
    fun::Function               # the limiter function
    generate_params::Function   # a function to map a DiscretisedFluidQueue to extra parameters for the limiter
end

muscl_params(dq::DiscretisedFluidQueue) =
    (vandermonde(n_bases_per_cell(dq))..., Δ(dq), cell_nodes(dq),n_phases(dq),N₋(dq),N₊(dq))

GeneralisedMUSCL = Limiter(limit,muscl_params)

NoLimiter = Limiter(identity,x->())


text:
abstract type SolverResult end
abstract type SolverType end
abstract type SolverSubtype end

struct Displacement <: SolverResult end

struct Direct <: SolverType end
struct CG  <: SolverType end

struct MatrixFree <: SolverSubtype end
struct Assembly <: SolverSubtype end

Utilities.getpenalty(solver::AbstractFEASolver) = solver.penalty
function Utilities.setpenalty!(solver::AbstractFEASolver, p)
    solver.prev_penalty = deepcopy(solver.penalty)
    if p isa AbstractPenalty
        solver.penalty = p
    elseif p isa Number
        setpenalty!(solver.penalty, p)
    else
        throw("Unsupported penalty value $p.")
    end
    return solver
end
Utilities.getprevpenalty(solver::AbstractFEASolver) = solver.prev_penalty

function FEASolver(::Type{Displacement}, ::Type{Direct}, problem; kwargs...)
    DirectDisplacementSolver(problem; kwargs...)
end

function FEASolver(::Type{Displacement}, ::Type{CG}, problem; kwargs...)
    FEASolver(Displacement, CG, MatrixFree, problem; kwargs...)
end

function FEASolver(::Type{Displacement}, ::Type{CG}, ::Type{MatrixFree}, problem; kwargs...)
    StaticMatrixFreeDisplacementSolver(problem; kwargs...)
end

function FEASolver(::Type{Displacement}, ::Type{CG}, ::Type{Assembly}, problem; kwargs...)
    PCGDisplacementSolver(problem; kwargs...)
end

function default_quad_order(problem)
    if TopOptProblems.getdim(problem) == 2 && TopOptProblems.nnodespercell(problem) in (3, 6) || TopOptProblems.getdim(problem) == 3 && TopOptProblems.nnodespercell(problem) in (4, 10)
        return 3
    end
    if TopOptProblems.getgeomorder(problem) == 2
        return 6
    else
        return 4
    end
end


text:
#=
try_catch_finally:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-08
=#

module TryCatchFinally
    # valuables
    str = "Hello Julia"

    # methods
    function get_char(idx)
        return str[idx]
    end

    # 1. find exception at try phase
    # 2. catch phase
    # 3. finally phase
    function catch_error()
        try
            println("get char[0]")
            char = get_char(0)
        catch
            println("Idx is out of bounds")
        finally
            println("get char[5]")
            char = get_char(5)
            println("char is $(char)")
            return char
        end
    end

    function main()
        catch_error()
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TryCatchFinally
    TryCatchFinally.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using Distributions
using Statistics
using LinearAlgebra

export gillespie, gillespie_tspan, gillespie_full, gillespie_avg, gillespie_avg_v2

"""
    pfsample(w)

    Select a single random integer betwenn 1 and length(w) with
    probabilities proportional to the weights given in w.

    Based on the function StatsBase.sample(w)

"""
function pfsample(w::Array{Float64,1})
    t = rand() * sum(w)
    n = length(w)
    i = 1
    cw = w[1]
    while cw < t && i < n
        i += 1
        @inbounds cw += w[i]
    end
    return i
end

"""
	calc_nu(num_species, s_i, r_i)

	Calculates the stochiometric matrix nu

	Uses as input
	num_species		number of species
	s_i				stochiometric matrix for product molecules
					(only interaction reaction)
	r_i				stochiometric matrix for reactant molecules
					(only interaction reaction)

	Returns the stochiometric matrix for the systems inclusive the baseline
	creation and annihilation reaction

"""
function calc_nu(num_species::Int64, s_i::Matrix{Int64}, r_i::Matrix{Int64})
	eye = Matrix{Int}(I, num_species, num_species)
	return [eye; -eye; (s_i-r_i)]
end

"""
	prop_fct(x, k, r_i)

	Calculates the propensity function (reaction probability vector)

	Uses as input
	x		vector of molecule numbers
	params	reaction rate parameters
	r_i		stochiometric matrix for reactant molecules
			(only interaction reaction)

	Returns the reaction probability vector

"""
function prop_fct(x, params::Vector{Vector{Float64}}, r_i::Matrix{Int64})
	num_int = length(params[3])
	Hint = ones(num_int)
	for m in 1:num_int
		for k in 1:length(params[1])
			for r in 0:r_i[m,k]-1
				@inbounds Hint[m] *= (x[k]-r)
			end
		end
	end
	return [params[1]; params[2].*x'; params[3].*Hint]
end

"""
	regular_timeGrid(res, ta)

	Takes the output from the irregular time discretization
	of the Gillespie algorithm and makes the time grid regular

	Uses as input
	res		Result struct from the gillespie function
	ta		Desired time discretization

	Returns the evenly spaced time discretized Result

"""
function regular_timeGrid(res::Result, ta::Vector{Float64})

	num_species = size(res.data, 2)		# number of species
	l_ta = length(ta)					# length of time grid

	erg = zeros(l_ta, num_species)		# result
	i = 1
	j = 1

	# loop
	while i <= l_ta
		if ta[i] <= res.time[j+1]
			erg[i,:] = res.data[j,:]
		else
			i -= 1
			j += 1
		end
		i += 1
	end

	return Result(ta, erg)
end

"""
    gillespie(x0, params, nu, r_i, tf)

    Gillespie simulation of one specific trajectory and
    irregular time discretization.

    Uses as input
    x0      Vector of molecule numbers
    params  Vector of reaction constants
    nu		Stoichiometric matrix, reactions in rows
    r_i		Reactant coefficients
    tf      Final time when the reaction has to stop

    Returns the result in a struct with time and copy
    number of molecules.

"""
function gillespie(x0::Vector{Int64}, params::Vector{Vector{Float64}}, nu::Matrix{Int64}, r_i::Matrix{Int64}, tf::Float64)

    ta = Vector{Float64}()
    t = 0.0
    push!(ta,t)

    # Set up initial x
    nstates = length(x0)
    x = copy(x0')
    xa = copy(x0)

    # Number of reactions
    numpf = size(nu, 1)

    nsteps = 0
    while t <= tf
        pf = prop_fct(x,params,r_i)
        sumpf = sum(pf)
        if sumpf == 0.0
            # break condition
            t += 2*tf
        end
        dt = rand(Exponential(1/sumpf))
        t += dt
        push!(ta,t)
        # update event
        ev = pfsample(pf)
        deltax = view(nu,ev,:)
        for i in 1:nstates
            @inbounds x[1,i] += deltax[i]
        end
        for xx in x
            push!(xa,xx)
        end
        # update nsteps
        nsteps += 1
    end
    xar = transpose(reshape(xa,length(x),nsteps+1))
    return Result(ta,xar)
end

"""
    gillespie_tspan(x0, params, nu, r_i, time_grid)

    Gillespie simulation of one specific trajectory and
    even spaced time discretization.

    Uses as input
    x0      Vector of molecule numbers
    params  Vector of reaction constants
    nu      Stoichiometric matrix, reactions in rows
    r_i		Reactant coefficients
    time_grid	Vector of time discretization

    Returns the result in a struct with time and copy
    number of molecules.

"""
function gillespie_tspan(x0::Vector{Int64}, params::Vector{Vector{Float64}}, nu::Matrix{Int64}, r_i::Matrix{Int64}, time_grid::Vector{Float64})

    # Set up initial x
    nstates = length(x0)
    x = copy(x0')
    xnew = copy(x)
    xa = Vector{Int64}()

    # Number of reactions
    numpf = size(nu, 1)

    tm = time_grid[1]
    tf = time_grid[end]

    for t in time_grid

		while tm <= t

			x = copy(xnew)
		    pf = prop_fct(x,params,r_i)
		    sumpf = sum(pf)

		    if sumpf == 0.0
		    	# break condition
		        tm += 2*length(time_grid)
		    end

		    # update event
	        ev = pfsample(pf)
	        deltax = view(nu,ev,:)
	        for i in 1:nstates
	            @inbounds xnew[1,i] += deltax[i]
	        end

		    tm += rand(Exponential(1/sumpf))

		end

		for xx in x
			push!(xa,xx)
		end

	end

	return transpose(reshape(xa,length(x),length(time_grid)))

end

"""
    gillespie_full(w0, params, s_i, r_i, ta, N)

    Full N Gillespie Simulations with saved copy numbers.

    Uses as input
    w0      Initial conditions drawn from a Poisson
            distribution with means w0
    params  Vector of reaction constants
    s_i,r_i Stoichiometric matrix, reactions in rows
    ta		Vector of time discretization
    num_runs   Number of runs/repetitions

    Returns the full Gillespie Simulation of N runs

"""
function gillespie_full(w0::Vector{Float64}, params::Vector{Vector{Float64}}, s_i::Matrix{Int64}, r_i::Matrix{Int64}, ta::Vector{Float64}, num_runs::Int64)

	# Stochiometric matrix
	nu = calc_nu(length(params[1]), s_i, r_i)

    # Array to store results
    copyN = zeros(Int32, length(ta), size(s_i,2), num_runs)

    # loop over runs
    @simd for i in 1:num_runs
	    x0 = rand.(Poisson.(w0))
	    erg = gillespie_tspan(x0, params, nu, r_i, ta)
	    copyN[:,:,i] = erg[:,:]
    end

	return copyN
end

function gillespie_full(p::Parameters, num_runs::Int64)
	return gillespie_full(p.x0, p.k, p.s_i, p.r_i, tspan(p), num_runs)
end

"""
    gillespie_avg(w0, params, s_i, r_i, ta, N)

    Average of N Gillespie simulations.

    Same as gillespie_full, but averages over the N simulations.

    Returns the mean result in a struct with time and
    mean copy number of molecules.

"""
function gillespie_avg(w0::Vector{Float64}, params::Vector{Vector{Float64}}, s_i::Matrix{Int64}, r_i::Matrix{Int64}, ta::Vector{Float64}, num_runs::Int64)

    copyN = gillespie_full(w0, params, s_i, r_i, ta, num_runs)
    return Result(ta, reshape(mean(copyN, dims=3), length(ta), size(s_i,2))')

end

function gillespie_avg(p::Parameters, num_runs::Int64)
	return gillespie_avg(p.x0, p.k, p.s_i, p.r_i, tspan(p), num_runs)
end

"""
	gillespie_avg_v2(w0, params, s_i, r_i, ta, num_runs)

	Average of N Gillespie simulations.

	Same as gillespie_avg, but the regular time grid is taken after
	the gillespie simulation and not within.

	Returns the mean result in a struct with time and
	mean copy number of molecules.

"""
function gillespie_avg_v2(w0::Vector{Float64}, params::Vector{Vector{Float64}}, s_i::Matrix{Int64}, r_i::Matrix{Int64}, ta::Vector{Float64}, num_runs::Int64)

	# Stochiometric matrix
	nu = calc_nu(length(params[1]), s_i, r_i)

    # Array to store results
    copyN = zeros(Int32, length(ta), size(s_i,2), num_runs)

    # loop over runs
    @simd for i in 1:num_runs
	    x0 = rand.(Poisson.(w0))
	    gil = gillespie(x0, params, nu, r_i, ta[end])
	    erg = regular_timeGrid(gil, ta).data
	    copyN[:,:,i] = erg[:,:]
    end

	return Result(ta, reshape(mean(copyN, dims=3), length(ta), size(s_i,2))')

end

function gillespie_avg_v2(p::Parameters, num_runs::Int64)
	return gillespie_avg_v2(p.x0, p.k, p.s_i, p.r_i, tspan(p), num_runs)
end

max_stars_repo_name: mkurjahn/CRN

max_issues_repo_name: mkurjahn/CRN

max_forks_repo_name: mkurjahn/CRN


text:
abstract type CTmStruct end

const z32 = zero(Int32)

mutable struct TmCStruct <: CTmStruct
    sec  ::Int32 # 0..
    min  ::Int32 # 0..
    hour ::Int32 # 0..
    mday ::Int32 # 1..
    month::Int32 # 0..
    year ::Int32 # 1900=>0..
    wday ::Int32 # 0==Sunday..
    yday ::Int32 # 0..
    isdst::Int32
    # on some platforms the struct is 14 wordsc_gettimeofday_sec, , even though 9 are specified
    _10::Int32
    _11::Int32
    _12::Int32
    _13::Int32
    _14::Int32

    TmCStruct() = new(z32,z32,z32,z32,z32,z32,z32,z32,z32,z32,z32,z32,z32,z32)
end

function (==)(a::TmCStruct, b::TmCStruct)
    fields = (:sec,:min, :hour, :mday, :month, :year, :wday, :yday, :isdst)
    res = true
    for f in fields
       if (a).(f) != (b).(f)
          print(f)
          res = false
       end
    end
    res
end

(==)(a::CTmStruct, b::CTmStruct) = (==)(a.tm, b.tm)

mutable struct LclTmStruct <: CTmStruct
    tm::TmCStruct

    LclTmStruct() = new(TmCStruct())
    LclTmStruct(tm::TmCStruct) = new(tm)
    function LclTmStruct(t::Real)
        t = floor(t)
        tm = TmCStruct()
        ccall(:localtime_r, Ptr{TmCStruct}, (Ptr{Int}, Ptr{TmCStruct}), &t, &tm)
        new(tm)
    end
end

mutable struct UtcTmStruct <: CTmStruct
    tm::TmCStruct

    UtcTmStruct() = new(TmCStruct())
    UtcTmStruct(tm::TmCStruct) = new(tm)
    function UtcTmStruct(t::Real)
        t = floor(t)
        tm = TmCStruct()
        ccall(:gmtime_r, Ptr{TmCStruct}, (Ptr{Int}, Ptr{TmCStruct}), &t, &tm)
        new(tm)
    end
end

function LclTmStruct(year::T,month=1::T,day=1::T,hour=0::T,min=0::T,sec=0::T;isdst=(-1)::T) where T<:Integer
  lts = LclTmStruct()
  lts.tm.year  = Int32(year-1900)
  lts.tm.month = Int32(month-1)
  lts.tm.mday  = Int32(day)
  lts.tm.hour  = Int32(hour)
  lts.tm.min   = Int32(min)
  lts.tm.sec   = Int32(sec)
  lts.tm.isdst = Int32(isdst)
  lts
end

function UtcTmStruct(year::T,month=1::T,day=1::T,hour=0::T,min=0::T,sec=0::T;isdst=0::T) where T<:Integer
  uts = UtcTmStruct()
  uts.tm.year  = Int32(year-1900)
  uts.tm.month = Int32(month-1)
  uts.tm.mday  = Int32(day)
  uts.tm.hour  = Int32(hour)
  uts.tm.min   = Int32(min)
  uts.tm.sec   = Int32(sec)
  uts.tm.isdst = Int32(isdst)
  uts
end

#=
  unbreak_{ut,lcl,std,dst}time: CTmStruct --> seconds after POSIX Epoch
  breakdn_{ut,lcl,std,dst}time: seconds after POSIX Epoch --> CTmStruct
=#

function breakdn_lcltime(sec::Int64)
    s   = sec
    lts = LclTmStruct()
    ccall(:localtime_r, Ptr{TmCStruct}, (Ptr{Int64}, Ptr{TmCStruct}), &s, &lts.tm)
    return lts
end
breakdn_lcltime(t::Real) = breakdn_lcltime(floor(Int64,t))

breakdn_stdtime(sec::Int64) = breakdn_lcltime(unbreak_stdtime(breakdn_lcltime(sec)))
breakdn_stdtime(t::Real) = breakdn_stdtime(floor(Int64,t))

breakdn_dsttime(sec::Int64) = breakdn_lcltime(unbreak_dsttime(breakdn_lcltime(sec)))
breakdn_dsttime(t::Real) = breakdn_dsttime(floor(Int64,t))

function breakdn_utmtime(sec::Int64)
    s  = sec
    uts = UtcTmStruct()
    ccall(:gmtime_r,Ptr{TmCStruct}, (Ptr{Int64}, Ptr{TmCStruct}), &s, &uts.tm)
    return uts
end
breakdn_utmtime(t::Real) = breakdn_utmtime(floor(Int64,t))

function unbreak_lcltime(cts::CTmStruct)
    tcs = cts.tm
    tcs.isdst = -one(Int32)
    s = ccall(:mktime, Int64, (Ptr{TmCStruct},), &tcs)
    if s == typemax(Int64) - one(Int64)
       throw(ArgumentError("unrepresentable mktime(tm=$(tm))"))
    end
    return s
end

function unbreak_stdtime(cts::CTmStruct)
    tcs = cts.tm
    tcs.isdst = zero(Int32)
    s = ccall(:mktime, Int64, (Ptr{TmCStruct},), &tcs)
    if s == typemax(Int64) - one(Int64)
       throw(ArgumentError("unrepresentable mktime(tm=$(tm))"))
    end
    return s
end

function unbreak_dsttime(cts::CTmStruct)
    tcs = cts.tm
    tcs.isdst = one(Int32)
    s = ccall(:mktime, Int64, (Ptr{TmCStruct},), &tcs)
    if s == typemax(Int64) - one(Int64)
       throw(ArgumentError("unrepresentable mktime(tm=$(tm))"))
    end
    return s
end

#=
unbreak_utmtime adapted from public domain source code
by Eric S. Raymond <esr@thyrsus.com>
at http://www.catb.org/esr/time-programming
=#
const cumdays = Int64[ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 ];

function unbreak_utmtime(cts::CTmStruct)
    mon   = Int64(cts.tm.month)
    mon12 = mon % 12
    year  = 1900 + Int64(cts.tm.year) + div(mon,12)
    ans = (year - 1970) * 365 + cumdays[1+mon12];
    ans += div(year - 1968, 4)
    ans -= div(year - 1900, 100)
    ans += div(year - 1600, 400)
    if ((year % 4) == 0 && ((year % 100) != 0 || (year % 400) == 0) && (mon12 < 2))
        ans -= 1
    end
    ans += Int64(cts.tm.mday) - 1
    ans *= 24
    ans += Int64(cts.tm.hour)
    ans *= 60;
    ans += Int64(cts.tm.min)
    ans *= 60;
    ans += Int64(cts.tm.sec)
    if (cts.tm.isdst == one(Int32))
        ans -= 3600
    end
    ans
end

# TECHNICALITY
#   these two functions are mutually inversive
#   their implementations may not be inverses

utm_from_lcltime(sec::Int64) = unbreak_lcltime(breakdn_utmtime(sec))
lcl_from_utmtime(sec::Int64) = unbreak_utmtime(breakdn_stdtime(sec))

utm_offset(sec::Int64) = sec - utm_from_lcltime(sec)
lcl_offset(sec::Int64) = -utm_offset(sec)

function utm_from_lcltime(sec::Float64)
    isec = floor(Int64, sec)
    utm = unbreak_lcltime(breakdn_utmtime(isec))

    fsec = signbit(lcl_offset(isec)) ? isec-sec : sec-isec
    Base.Dates.unix2datetime(Base.Dates.datetime2unix(utm) + fsec)
end

function lcl_from_utmtime(sec::Float64)
    isec = floor(Int64, sec)
    lcl = unbreak_utmtime(breakdn_stdtime(isec))
    lcl # adjust for fractional secs +/- wrt local timezone
end

utm_unwrapped() =
    Base.Dates.unix2datetime(
        utm_from_lcltime(
            floor(Int64,Base.Dates.datetime2unix(Base.Dates.now())) ))

# dt as localtime --> ut (without leap seconds)
function utm_unwrapped(dt::DateTime)
    sec = floor(Int64, Base.Dates.datetime2unix(dt))
    sec = utm_from_lcltime(sec)
    Base.Dates.unix2datetime( sec )
end


text:
###############################################################################
#
# Create DataFrame to save simulation output
#
###############################################################################
"""
    function create_dataframe(FT, weather::DataFrame)

Create a data frame to store simulation output, given
- `FT` Floating number type
- `weather` Weather profile in a growing season
"""
function create_dataframe(
            FT,
            weather::DataFrame
)
    df    = DataFrame();

    # climatic info
    df[!, "Time"  ]  = (weather).Day + (weather).Hour / FT(24);
    df[!, "T_air" ] .= FT(0);
    df[!, "D_air" ] .= FT(0);
    df[!, "Wind"  ] .= FT(0);
    df[!, "Rain"  ] .= FT(0);
    df[!, "Ca"    ] .= FT(0);

    # whole tree level
    df[!, "SWC"   ] .= FT(0);
    df[!, "P_soil"] .= FT(0);
    df[!, "H_sun" ] .= FT(0);
    df[!, "A_net" ] .= FT(0);
    df[!, "E_crit"] .= FT(0);

    # sunlit layers
    df[!, "LAI_sl"] .= FT(0);
    df[!, "PAR_sl"] .= FT(0);
    df[!, "RAD_sl"] .= FT(0);
    df[!, "E_sl"  ] .= FT(0);
    df[!, "P_sl"  ] .= FT(0);
    df[!, "An_sl" ] .= FT(0);
    df[!, "Ag_sl" ] .= FT(0);
    df[!, "C_sl"  ] .= FT(0);
    df[!, "G_sl"  ] .= FT(0);
    df[!, "T_sl"  ] .= FT(0);

    # shaded layers
    df[!, "LAI_sh"] .= FT(0);
    df[!, "PAR_sh"] .= FT(0);
    df[!, "RAD_sh"] .= FT(0);
    df[!, "E_sh"  ] .= FT(0);
    df[!, "P_sh"  ] .= FT(0);
    df[!, "An_sh" ] .= FT(0);
    df[!, "Ag_sh" ] .= FT(0);
    df[!, "C_sh"  ] .= FT(0);
    df[!, "G_sh"  ] .= FT(0);
    df[!, "T_sh"  ] .= FT(0);

    return df
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module ModuleMoments
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators
using StaticArrays
export moments
include("moments_fundamentalmatrix.jl")
include("moments_homographymatrix.jl")
end


text:
using Printf
using GLMakie
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity

Nx = 128
Ny = 128
Nz = 32

Lx = Ly = 100kilometers
Lz = 200

grid = RectilinearGrid(size=(Nx, Ny, Nz), halo=(3, 3, 3),
                       x=(0, Lx), y=(0, Ly), z=(-Lz, 0),
                       topology=(Periodic, Bounded, Bounded))

@show grid

Qᵇ(x, y, t) = 2e-8 * y / Ly
Qᵘ(x, y, t) = -1e-5
b_top_bc = FluxBoundaryCondition(Qᵇ)
u_top_bc = FluxBoundaryCondition(Qᵘ)
b_bcs = FieldBoundaryConditions(top=b_top_bc)
u_bcs = FieldBoundaryConditions(top=u_top_bc)

Δh = Ly / Ny
κ₄ = Δh^4 / 1day
κ₂ = Δh^2 / 10days

horizontal_closure = HorizontalScalarBiharmonicDiffusivity(ν=κ₄, κ=κ₄)
#horizontal_closure = HorizontalScalarDiffusivity(ν=κ₂, κ=κ₂)
boundary_layer_closure = CATKEVerticalDiffusivity()
#boundary_layer_closure = ConvectiveAdjustmentVerticalDiffusivity(convective_κz=0.1)

closure = (boundary_layer_closure, horizontal_closure)
#closure = boundary_layer_closure

model = HydrostaticFreeSurfaceModel(; grid, closure,
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    tracers = (:b, :e),
                                    boundary_conditions = (; b=b_bcs, u=u_bcs),
                                    buoyancy = BuoyancyTracer())

N² = 1e-5
h = Lz / 3
bᵢ(x, y, z) = N² * z #N² * h * exp(z / h) # + 1e-8 * rand()
set!(model, b=bᵢ)

simulation = Simulation(model, Δt=1minutes, stop_iteration=10000)

#=
slice_indices = (
    west   = (1,  :, :),
    east   = (Nx, :, :),
    south  = (:, 1,  :),
    north  = (:, Ny, :),
    bottom = (:, :, 1),
    top    = (:, :, Nz),
)

for side in keys(indices)
    indices = slice_indices[side]

    simulation.output_writers[:side] =
        JLD2OutputWriter(model, merge(model.velocities, model.tracers);
                         schedule = TimeInterval(10minutes),
                         prefix = "heterogeneous_cooling_" * side,
                         indices,
                         overwrite_existing = true)
end
=#

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers);
                     schedule = TimeInterval(5minutes),
                     prefix = "heterogeneous_cooling",
                     overwrite_existing = true)

function progress(sim)
    u, v, w = sim.model.velocities
    e = sim.model.tracers.e

    msg = @sprintf("Iter: %d, t: %s, Δt: %s, max|u|: (%6.2e, %6.2e, %6.2e) m s⁻¹",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w))

    msg *= @sprintf(", max(e): %6.2e m² s⁻²", maximum(abs, e))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(1))

run!(simulation)

filepath = "heterogeneous_cooling.jld2"
b_ts = FieldTimeSeries(filepath, "b")
e_ts = FieldTimeSeries(filepath, "e")
u_ts = FieldTimeSeries(filepath, "u")
v_ts = FieldTimeSeries(filepath, "v")
Nt = length(b_ts.times)

fig = Figure(resolution=(1600, 1200))

ax_bxy = Axis(fig[1, 1], aspect=1)
ax_exy = Axis(fig[1, 2], aspect=1)
ax_syz = Axis(fig[2, 1], aspect=2)
ax_eyz = Axis(fig[2, 2], aspect=2)
slider = Slider(fig[3, :], range=1:Nt, startvalue=1)
n = slider.value

b_xy = @lift interior(b_ts[$n], :, :, Nz)
b_xz = @lift interior(b_ts[$n], :, 1, :)
b_yz = @lift interior(b_ts[$n], 1, :, :)

e_xy = @lift interior(e_ts[$n], :, :, Nz)
e_xz = @lift interior(e_ts[$n], :, 1, :)
e_yz = @lift interior(e_ts[$n], 1, :, :)

x, y, z = nodes(b_ts)

shear_yz = @lift begin
    u = u_ts[$n]
    v = v_ts[$n]
    shear_u = compute!(Field(∂z(u)^2))
    shear_v = compute!(Field(∂z(v)^2))
    shear_op = @at (Center, Center, Center) sqrt(shear_u + shear_v)
    shear = compute!(Field(shear_op))
    interior(shear, 1, :, :)
end

heatmap!(ax_bxy, x, y, b_xy)
heatmap!(ax_exy, x, y, e_xy)

heatmap!(ax_syz, y, z, shear_yz)
contour!(ax_syz, y, z, b_yz, levels=15)

heatmap!(ax_eyz, y, z, e_yz)
contour!(ax_eyz, y, z, b_yz, levels=15)

display(fig)

#=
ax = Axis3(fig[1, 1], aspect=:data)

b = model.tracers.b
bxy = interior(b, :, :, Nz)
bxz = interior(b, :, 1, :)
byz = interior(b, 1, :, :)

x, y, z = nodes(b)

x_xz = repeat(x, 1, Nz)
y_xz = 0.995 * Ly * ones(Nx, Nz)
z_xz = repeat(reshape(z, 1, Nz), Nx, 1)

x_yz = 0.995 * Lx * ones(Ny, Nz)
y_yz = repeat(y, 1, Nz)
z_yz = repeat(reshape(z, 1, Nz), Ny, 1)

x_xy = x
y_xy = y
z_xy = -0.001 * Lz * ones(Nx, Ny)

surface!(ax, x_xz, y_xz, z_xz, color=bxz)
surface!(ax, x_yz, y_yz, z_yz, color=byz)
surface!(ax, x_xy, y_xy, z_xy, color=bxy)
=#

display(fig)


text:
# This file calculates scattering from the Helmholtz equation via the adaptive spectral method.
# Δu + k^2u = 0,
# ∂u∂n|Γ = 0,
# ∂u^i/∂n = -im k d[2] e^{im k x⋅d},
# u = u^i + u^s.
# The normal derivative ∂u/∂n of the entire wave is calculated on the sound-hard boundaries.
# The scattered wave is calculated by convolving [u] with the fundamental solution.
# Then, the total wave is obtained by summing the incident and the scattered waves.

using ApproxFun, SingularIntegralEquations
include("Scatteraux.jl")

k = 50.
ω = 2π
d = (1,-1)
d = d[1]/hypot(d[1],d[2]),d[2]/hypot(d[1],d[2])
ui(x,y) = exp(im*k*(d⋅(x,y)))

# The Helmholtz Green's function, split into singular and nonsingular pieces.
g1(x,y) = besselj0(k*abs(y-x))/2
g2(x,y) = x == y ? -k^2/4 : -k*besselj1(k*abs(y-x))./abs(y-x)/2
g3(x,y) = g3neumann(x,y) # In /Scatteraux.jl
g4old(x,y) = im*k/4*hankelh1(1,k*abs(y-x))./abs(y-x).*imag(y-x)
g4(x,y) = im*k/4*besselj1(k*abs(y-x))./abs(y-x).*imag(y-x)  # For linesum
g5(x,y) = -k/2*besselj1(k*abs(y-x))./abs(y-x).*imag(y-x)  # For logkernel
g6(x,y) = k/2*abs(y-x).*(bessely1(k*abs(y-x)) - 2besselj1(k*abs(y-x)).*logabs(y-x)/π) # For Re{Cauchy}

dom = ∪(Interval,[-2.,1.],[-1.,2.])
sp = Space(dom)
cwsp,cwsp2 = CauchyWeight(sp⊗sp,0),CauchyWeight(sp⊗sp,2)
∂ui∂nΓ,⨍ = Fun(t->-im*k*d[2]*ui(real(t),imag(t)),sp),
                DefiniteLineIntegral(PiecewiseSpace(map(d->JacobiWeight(.5,.5,Ultraspherical(1,d)),dom.domains)))#DefiniteLineIntegral(dom)

@time G = GreensFun(g1,cwsp2;method=:Cholesky) + GreensFun(g2,cwsp;method=:Cholesky) + GreensFun(g3,sp⊗sp;method=:Cholesky)

@time u = ⨍[G]\-∂ui∂nΓ
println("The length of u is: ",length(u))
us(x,y) = linesum(g4,u,complex(x,y))+logkernel(g5,u,complex(x,y))+π*real(cauchy(g6,real(u),complex(x,y)))+π*im*real(cauchy(g6,imag(u),complex(x,y)))
#dom += 0im


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# The concept and trivial implementation of graph visitors

abstract type AbstractGraphVisitor end

# trivial implementation

# invoked when a vertex v is encountered for the first time
# this function returns whether to continue search
discover_vertex!(vis::AbstractGraphVisitor, v) = true

# invoked when the algorithm is about to examine v's neighbors
open_vertex!(vis::AbstractGraphVisitor, v) = true

# invoked when a neighbor is discovered & examined
examine_neighbor!(vis::AbstractGraphVisitor, u, v, ucolor::Int, vcolor::Int, ecolor::Int) = true

# invoked when all of v's neighbors have been examined
close_vertex!(vis::AbstractGraphVisitor, v) = true

struct TrivialGraphVisitor <: AbstractGraphVisitor
end

# This is the common base for BreadthFirst and DepthFirst
abstract type AbstractGraphVisitAlgorithm end

const AbstractEdgeMap{T} = Associative{Edge,T}
const AbstractVertexMap{T<:Integer,U} = Union{AbstractVector{T},Associative{T,U}}

struct DummyEdgeMap <: AbstractEdgeMap{Int}
end

getindex(d::DummyEdgeMap, e::Edge) = 0
setindex!(d::DummyEdgeMap, x::Int, e::Edge) = x
get(d::DummyEdgeMap, e::Edge, x::Int) = x
Base.show(io::IO, ::MIME"text/plain", ::DummyEdgeMap) =
    write(io, "LightGraphs.DummyEdgeMap()") # for displaying at the REPL
Base.show(io::IO, ::LightGraphs.DummyEdgeMap) =
    write(io, "LightGraphs.DummyEdgeMap()") # for string(DummyEdgeMap)

###########################################################
#
#   General algorithms based on graph traversal
#
###########################################################

# List vertices by the order of being discovered

struct VertexListVisitor{T<:Integer} <: AbstractGraphVisitor
    vertices::Vector{T}
end

function VertexListVisitor(n::T=0) where T<:Integer
    vs = Vector{T}()
    sizehint!(vs, n)
    return VertexListVisitor(vs)
end

function discover_vertex!(visitor::VertexListVisitor, v::Integer)
    push!(visitor.vertices, v)
    return true
end

function visited_vertices(
    g::AbstractGraph,
    alg::AbstractGraphVisitAlgorithm,
    sources)
    T = eltype(g)
    visitor = VertexListVisitor(nv(g))
    traverse_graph!(g, alg, sources, visitor)
    visitor.vertices::Vector{T}
end

# Print visit log

struct LogGraphVisitor{S<:IO} <: AbstractGraphVisitor
    io::S
end

function discover_vertex!(vis::LogGraphVisitor, v::Integer)
    println(vis.io, "discover vertex: $v")
    return true
end

function open_vertex!(vis::LogGraphVisitor, v::Integer)
    println(vis.io, "open vertex: $v")
    return true
end

function close_vertex!(vis::LogGraphVisitor, v::Integer)
    println(vis.io, "close vertex: $v")
    return true
end

function examine_neighbor!(vis::LogGraphVisitor, u::Integer, v::Integer, ucolor::Int, vcolor::Int, ecolor::Int)
    println(vis.io, "examine neighbor: $u -> $v (ucolor = $ucolor, vcolor = $vcolor, edgecolor= $ecolor)")
    return true
end

function traverse_graph_withlog(
    g::AbstractGraph,
    alg::AbstractGraphVisitAlgorithm,
    sources,
    io::IO = STDOUT
)
    visitor = LogGraphVisitor(io)
    traverse_graph!(g, alg, sources, visitor)
end


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior: parseonly
using Behavior.Gherkin: ParseOptions

argslength = length(ARGS)
if argslength !== 1 && argslength !== 2
    println("Usage: julia parseonly.jl <root-directory> [--experimental]")
    exit(1)
end

use_experimental = argslength == 2 && ARGS[2] == "--experimental"
parseoptions = ParseOptions(allow_any_step_order=true, use_experimental=use_experimental)

results = parseonly(ARGS[1], parseoptions=parseoptions)
num_files = length(results)
num_success = count(x -> x.success === true, results)
if num_success === num_files
    println("All good!")
else
    num_failed = num_files - num_success
    println("Files failed parsing:")
    for rs in results
        print(rs.filename)

        if rs.success
            println(": OK")
        else
            if use_experimental
                println(": NOT OK")
                println(rs.result)
            else
                println()
                println(" reason: ", rs.result.reason)
                println(" expected: ", rs.result.expected)
                println(" actual: ", rs.result.actual)
                println(" line $(rs.result.linenumber): $(rs.result.line)")
            end
        end
    end
    println("Parsing failed: ", num_failed)
    println("Total number of files: ", num_files)
end


text:
using Test
using PolynomialBasis
# using Revise
using CutCellDG
include("../useful_routines.jl")

polyorder = 1
delta = 0.1
numqp = 2

L, W = 4.0, 2.0

solverbasis = LagrangeTensorProductBasis(2, polyorder)
refpoints = interpolation_points(solverbasis)
dim, nf = size(refpoints)
cgmesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [2, 1], nf)
dgmesh = CutCellDG.DGMesh([0.0, 0.0], [L, W], [2, 1], refpoints)

normal = [1.0, 0.0]
x0 = [2.0 + delta, 0.0]
distancefunction(x) = plane_distance_function(x, normal, x0)
levelset = CutCellDG.LevelSet(distancefunction, cgmesh, solverbasis)

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)

refseedpoints, refseedcellids =
    CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialseedpoints =
    CutCellDG.map_to_spatial(refseedpoints, refseedcellids, cutmesh)

cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)
interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)

mergedmesh =
    CutCellDG.MergedMesh!(cutmesh, cellquads, facequads, interfacequads)
@assert CutCellDG.has_merged_cells(mergedmesh)

querypoints = reshape([4.0, 1.0], 2, 1)
tol = 1e4eps()
boundingradius = 1.5 * maximum(CutCellDG.element_size(dgmesh))
closestpoints, closestcellids, flags =
    CutCellDG.closest_points_on_zero_levelset(
        querypoints,
        spatialseedpoints,
        refseedcellids,
        levelset,
        tol,
        boundingradius,
    )

@test all(flags)

parentclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    -1,
    mergedmesh,
)
productclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    +1,
    mergedmesh,
)

@test allapprox(parentclosestrefpoints, [1.1, 0.0], 1e2eps())
@test allapprox(productclosestrefpoints, [-0.9, 0.0], 1e2eps())
################################################################################

################################################################################
levelsetbasis = HermiteTensorProductBasis(2)
solverbasis = LagrangeTensorProductBasis(2, polyorder)
refpoints = interpolation_points(solverbasis)
dim, nf = size(interpolation_points(levelsetbasis))
cgmesh = CutCellDG.CGMesh([0.0, 0.0], [L, W], [2, 1], nf)
dgmesh = CutCellDG.DGMesh([0.0, 0.0], [L, W], [2, 1], refpoints)

normal = [1.0, 0.0]
x0 = [2.0 + delta, 0.0]
distancefunction(x) = plane_distance_function(x, normal, x0)
levelset = CutCellDG.LevelSet(distancefunction, cgmesh, levelsetbasis)

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)

refseedpoints, refseedcellids =
    CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialseedpoints =
    CutCellDG.map_to_spatial(refseedpoints, refseedcellids, cutmesh)

cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)
interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)

mergedmesh =
    CutCellDG.MergedMesh!(cutmesh, cellquads, facequads, interfacequads)
@assert CutCellDG.has_merged_cells(mergedmesh)

querypoints = reshape([4.0, 1.0], 2, 1)
tol = 1e4eps()
boundingradius = 1.5 * maximum(CutCellDG.element_size(dgmesh))
closestpoints, closestcellids, flags =
    CutCellDG.closest_points_on_zero_levelset(
        querypoints,
        spatialseedpoints,
        refseedcellids,
        levelset,
        tol,
        boundingradius,
    )

@test all(flags)

parentclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    -1,
    mergedmesh,
)
productclosestrefpoints = CutCellDG.map_to_reference_on_merged_mesh(
    closestpoints,
    closestcellids,
    +1,
    mergedmesh,
)

@test allapprox(parentclosestrefpoints, [1.1, 0.0], 1e3eps())
@test allapprox(productclosestrefpoints, [-0.9, 0.0], 1e3eps())


text:
Lar = LinearAlgebraicRepresentation

"""
	frag_face_channel(in_chan, out_chan, V, EV, FE, sp_idx)

Parallel fragmentation of faces in `FE` against faces in `sp_idx`.
"""
function frag_face_channel(in_chan, out_chan, V, EV, FE, sp_idx)
    run_loop = true
    while run_loop
        sigma = take!(in_chan)
        if sigma != -1
            put!(out_chan, frag_face(V, EV, FE, sp_idx, sigma))
        else
            run_loop = false
        end
    end
end

"""
	frag_face(V, EV, FE, sp_idx, sigma)

`sigma` face fragmentation against faces in `sp_idx[sigma]`
"""
function frag_face(V, EV, FE, sp_idx, sigma)

    vs_num = size(V, 1)

	# 2D transformation of sigma face
    sigmavs = (abs.(FE[sigma:sigma,:])*abs.(EV))[1,:].nzind
    sV = V[sigmavs, :]
    sEV = EV[FE[sigma, :].nzind, sigmavs]
    M = Lar.Arrangement.submanifold_mapping(sV)
    tV = ([V ones(vs_num)]*M)[:, 1:3]  # folle convertire *tutti* i vertici
    sV = tV[sigmavs, :]
    # sigma face intersection with faces in sp_idx[sigma]
    for i in sp_idx[sigma]
        tmpV, tmpEV = Lar.Arrangement.face_int(tV, EV, FE[i, :])

        sV, sEV = Lar.skel_merge(sV, sEV, tmpV, tmpEV)
    end

    # computation of 2D arrangement of sigma face
    sV = sV[:, 1:2]
    nV, nEV, nFE = planar_arrangement(sV, sEV, sparsevec(ones(Int8, length(sigmavs))))
    if nV == nothing ## not possible !! ... (each original face maps to its decomposition)
        return [], spzeros(Int8, 0,0), spzeros(Int8, 0,0)
    end
    nvsize = size(nV, 1)
    nV = [nV zeros(nvsize) ones(nvsize)]*inv(M)[:, 1:3] ## ????
    return nV, nEV, nFE
end

function merge_vertices(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp, err=1e-4)
    vertsnum = size(V, 1)
    edgenum = size(EV, 1)
    facenum = size(FE, 1)
    newverts = zeros(Int, vertsnum)
    # KDTree constructor needs an explicit array of Float64
    V = Array{Float64,2}(V)
    W = convert(Lar.Points, LinearAlgebra.transpose(V))
    kdtree = KDTree(W)
	# remove vertices congruent to a single representative
    todelete = []
    i = 1
    for vi in 1:vertsnum
        if !(vi in todelete)
            nearvs = Lar.inrange(kdtree, V[vi, :], err)
            newverts[nearvs] .= i
            nearvs = setdiff(nearvs, vi)
            todelete = union(todelete, nearvs)
            i = i + 1
        end
    end
    nV = V[setdiff(collect(1:vertsnum), todelete), :]

    # translate edges to take congruence into account
    edges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    oedges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    for ei in 1:edgenum
        v1, v2 = EV[ei, :].nzind
        edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
        oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))
    end
    nedges = union(edges)
    # remove edges of zero length
    nedges = filter(t->t[1]!=t[2], nedges)
    nedgenum = length(nedges)
    nEV = spzeros(Int8, nedgenum, size(nV, 1))

    etuple2idx = Dict{Tuple{Int, Int}, Int}()
    for ei in 1:nedgenum
    	begin
        	nEV[ei, collect(nedges[ei])] .= 1
        	nEV
        end
        etuple2idx[nedges[ei]] = ei
    end
    for e in 1:nedgenum
    	v1,v2 = findnz(nEV[e,:])[1]
    	nEV[e,v1] = -1; nEV[e,v2] = 1
    end

    # compute new faces to take congruence into account
    faces = [[
        map(x->newverts[x], FE[fi, ei] > 0 ? oedges[ei] : reverse(oedges[ei]))
        for ei in FE[fi, :].nzind
    ] for fi in 1:facenum]

    visited = []
    function filter_fn(face)

        verts = []
        map(e->verts = union(verts, collect(e)), face)
        verts = Set(verts)

        if !(verts in visited)
            push!(visited, verts)
            return true
        end
        return false
    end

    nfaces = filter(filter_fn, faces)

    nfacenum = length(nfaces)
    nFE = spzeros(Int8, nfacenum, size(nEV, 1))

    for fi in 1:nfacenum
        for edge in nfaces[fi]
            ei = etuple2idx[Tuple{Int, Int}(sort(collect(edge)))]
            nFE[fi, ei] = sign(edge[2] - edge[1])
        end
    end

    return Lar.Points(nV), nEV, nFE
end

function spatial_arrangement_1(
		V::Lar.Points,
		copEV::Lar.ChainOp,
		copFE::Lar.ChainOp, multiproc::Bool=false)

	# spaceindex computation
	FV = Lar.compute_FV( copEV, copFE )
	model = (convert(Lar.Points,V'), FV)
	sp_idx = Lar.spaceindex(model)

	# initializations
    fs_num = size(copFE, 1)
    rV = Array{Float64,2}(undef,0,3)
    rEV = SparseArrays.spzeros(Int8,0,0)
    rFE = SparseArrays.spzeros(Int8,0,0)

	# multiprocessing of face fragmentation
    if (multiproc == true)
        in_chan = Distributed.RemoteChannel(()->Channel{Int64}(0))
        out_chan = Distributed.RemoteChannel(()->Channel{Tuple}(0))
        @async begin
            for sigma in 1:fs_num
                put!(in_chan, sigma)
            end
            for p in Distributed.workers()
                put!(in_chan, -1)
            end
        end
        for p in Distributed.workers()
            @async Base.remote_do(
                frag_face_channel, p, in_chan, out_chan, V, EV, FE, sp_idx)
        end
        for sigma in 1:fs_num
            rV, rEV, rFE = Lar.skel_merge(rV, rEV, rFE, take!(out_chan)...)
        end
    else
	# sequential (iterative) processing of face fragmentation
        for sigma in 1:fs_num
            #print(sigma, "/", fs_num, "\r")
            nV, nEV, nFE = Lar.Arrangement.frag_face(V, copEV, copFE, sp_idx, sigma)
			#nV, nEV, nFE = Lar.fragface(V, copEV, copFE, sp_idx, sigma)
			nV = convert(Lar.Points, nV)
            a,b,c = Lar.skel_merge( rV,rEV,rFE,  nV,nEV,nFE )
            rV=a;  rEV=b;  rFE=c
        end
    end
	# merging of close vertices, edges and faces (3D congruence)
	rV, rEV, rFE = merge_vertices(rV, rEV, rFE)
    return rV, rEV, rFE
end

function spatial_arrangement_2(
		rV::Lar.Points,
		rcopEV::Lar.ChainOp,
		rcopFE::Lar.ChainOp, multiproc::Bool=false)

	#rcopCF = Lar.build_copFC(rV, rcopEV, rcopFE)
	rcopCF = minimal_3cycles(rV, rcopEV, rcopFE)
    return rV, rcopEV, rcopFE, rcopCF
end

"""
    spatial_arrangement(V::Points, copEV::ChainOp, copFE::ChainOp; [multiproc::Bool])

Compute the arrangement on the given cellular complex 2-skeleton in 3D.

A cellular complex is arranged when the intersection of every possible pair of cell
of the complex is empty and the union of all the cells is the whole Euclidean space.
The function returns the full arranged complex as a list of vertices V and a chain of borders EV, FE, CF.

## Additional arguments:
- `multiproc::Bool`: Runs the computation in parallel mode. Defaults to `false`.
"""
function spatial_arrangement(
		V::Lar.Points, # by rows
		copEV::Lar.ChainOp,
		copFE::Lar.ChainOp, multiproc::Bool=false)

	# face subdivision
	rV, rcopEV, rcopFE = Lar.Arrangement.spatial_arrangement_1( V, copEV, copFE, multiproc ) # copFE global
	bicon_comps = Lar.Arrangement.biconnected_components(rcopEV)
	#W,bicon_comps = Lar.biconnectedComponent((W,EV))
	#@error "comps# = $(length(bicon_comps))"
	# 3-complex and containment graph
	rV, rEV, rFE, rCF = Lar.Arrangement.spatial_arrangement_2(rV, rcopEV, rcopFE)
end


text:
function update_via_trained_prior!(sfs, dfs, cfs, tfs,
                                   batchprop::Float64,
                                   prior_thin::Int;
                                   batchsizes=nothing,
                                   fix, use_repulsive, Z_marg_lamgam, sb_ibp,
                                   time_updates, temper::Float64=1.0,
                                   verbose::Int=0,
                                   minibatch_update_all_params=false)

  # TODO:
  # Need to chane this to:
  # 1. Sample theta given minibatch
  #     - Note batchsize
  #     - Note thinning factor
  # 2. Accept theta (subset) with metropolis ratio
  # 3. Sample lambda
  # 4. Sample gamma
  # 5. Sample missing data

  # Generate minibatch
  if batchsizes == nothing
    batchsizes = round.(Int, dfs.data.N * batchprop)
  end
  idx_mini, idx_mega, d_mini, d_mega = sample_minibatch(batchsizes, dfs.data.y)
  dfs_mini = DataFS(d_mini, dfs.X)
  sfs_mini = deepcopy(sfs)
  for i in 1:dfs.data.I
    sfs_mini.theta.lam[i] = deepcopy(sfs_mini.theta.lam[i][idx_mini[i]])
    sfs_mini.theta.gam[i] = deepcopy(sfs_mini.theta.gam[i][idx_mini[i], :])
    sfs_mini.theta.y_imputed[i] = deepcopy(
      sfs.theta.y_imputed[i][idx_mini[i], :])
  end

  # Sample theta given minibatch
  for _ in 1:prior_thin
    if minibatch_update_all_params
      update_state_feature_select!(sfs_mini, cfs, dfs_mini, tfs,
                                   ll=zeros(0),
                                   fix=vcat(fix, :y_imputed, :marg_loglike),
                                   use_repulsive=use_repulsive,
                                   Z_marg_lamgam=Z_marg_lamgam,
                                   sb_ibp=sb_ibp, time_updates=time_updates)
    else
      update_Z_marg_lamgam!(sfs_mini.theta, cfs.constants, dfs_mini.data,
                            sb_ibp, use_repulsive=use_repulsive)
    end
  end

  # FIXME!
  # Compute metropolis ratio
  log_accept_ratio = let
    y_mega = [deepcopy(sfs.theta.y_imputed[i][idx_mega[i], :])
              for i in 1:dfs.data.I]
    c = cfs.constants
    ll_prop = compute_marg_loglike(sfs_mini.theta, c, d_mega, 1.0,
                                   y=y_mega, compute_prob_miss=false)
    ll_curr = compute_marg_loglike(sfs.theta, c, d_mega, 1.0,
                                   y=y_mega, compute_prob_miss=false)
    if verbose > 0
      print(" -- ll_prop: $(round(ll_prop, digits=2))")
      print(" -- ll_curr: $(round(ll_curr, digits=2))")
    end
    (ll_prop - ll_curr) / temper
  end

  if verbose > 1
    println()
    println("current mu0: $(-cumsum(sfs.theta.delta[0]))")
    println("current mu1: $(cumsum(sfs.theta.delta[1]))")
    println("current sig2: $(sfs.theta.sig2)")
    println("current w: $(sfs.theta.W)")
    println("proposed mu0: $(-cumsum(sfs_mini.theta.delta[0]))")
    println("proposed mu1: $(cumsum(sfs_mini.theta.delta[1]))")
    println("proposed sig2: $(sfs_mini.theta.sig2)")
    println("proposed w: $(sfs_mini.theta.W)")
  end

  if verbose > 0
    print(" -- log_acceptance_ratio: $(round(log_accept_ratio, digits=2))")
  end

  # Whether or not to accept proposal
  accept_proposal = log_accept_ratio > log(rand())

  if accept_proposal
    if verbose > 0
      print(" (accepted)")

      # Print updated Z
      if verbose > 0.5
        ks_old = findall(vec(sum(sfs.r, dims=1) .> 0))
        ks_new = findall(vec(sum(sfs_mini.r, dims=1) .> 0))
        Z_old = sfs.theta.Z[:, ks_old]
        Z_new = sfs_mini.theta.Z[:, ks_new]
        if Z_old != Z_new
          println("\nOld Z:")
          Base.show(stdout, "text/plain", Z_old)
          println("\nNew Z:")
          Base.show(stdout, "text/plain", Z_new)
        end
      end
    end

    if minibatch_update_all_params
      curr_lam = deepcopy(sfs.theta.lam)
      curr_gam = deepcopy(sfs.theta.gam)
      curr_y = deepcopy(sfs.theta.y_imputed)

      sfs.theta = sfs_mini.theta

      sfs.theta.lam .= curr_lam
      sfs.theta.gam .= curr_gam
      sfs.theta.y_imputed .= curr_y

      sfs.r .= deepcopy(sfs_mini.r)
      sfs.W_star .= deepcopy(sfs_mini.W_star)
      sfs.omega .= deepcopy(sfs_mini.omega)
    else
      sfs.theta.Z .= deepcopy(sfs_mini.theta.Z)
    end
  else
    if verbose > 0
      print(" (rejected)")
    end
  end

  if minibatch_update_all_params
    update_lam!(sfs.theta, cfs.constants, dfs.data)
    update_gam!(sfs.theta, cfs.constants, dfs.data)
    update_y_imputed!(sfs.theta, cfs.constants, dfs.data, tfs.tuners)
  else
    update_state_feature_select!(sfs, cfs, dfs, tfs,
                                 ll=zeros(0),
                                 fix=vcat(fix, :Z),
                                 use_repulsive=use_repulsive,
                                 Z_marg_lamgam=Z_marg_lamgam,
                                 sb_ibp=sb_ibp, time_updates=time_updates)
  end

  if verbose > 1
    println("\nomega: $(sfs.omega)")
    println("\np: $(compute_p(dfs.X, sfs.omega))")
    println("W:")
    Base.show(stdout, "text/plain", sfs.theta.W)
    println()
    println("R: $(sum(sfs.r, dims=2))")
  end

  if size(unique(sfs.theta.Z, dims=2), 2) < cfs.constants.K && use_repulsive
    println("This Z has repeated columns after updating!")
    Base.show(stdout, "text/plain", sfs.theta.Z)
    println("\nFlipping $(cfs.constants.K) bits")
    sfs.theta.Z .= flip_bits(sfs.theta.Z, cfs.constants.K)
  end
end


text:
using Juriba.Radiation
using Juriba.Constants
using StaticArrays, QuadGK, Interpolations

const RD = Juriba.Radiation

@Base.kwdef mutable struct Compton <: AbstractRadiation

    size::Int					        # size of arrays
    en_phot::Vector{Float64}			# array of photon energies
    num_phot::Vector{Float64} 			# array of number of photons in units of erg/s/Hz
    en_phot_obs::Vector{Float64}		# same as above but in observer frame
    num_phot_obs::Vector{Float64} 		# same as above but in observer frame

    seed_ph::Interpolations.Extrapolation
    iter_ph::InterruptException.Extrapolation
    esc_p_sph::Interpolations.Extrapolation
    esc_p_cyl::Interpolation.Extrapolation

    z::Float64 = 0.				            # Dimensions of emitting region
    r::Float64 = 0.
    vol::Float64 = 0.					    # Volume of emitting region
    beta::Float64 = 0.				        # speed of the emitting region
    angle::Float64 = 0.		        # Viewing angle/Doppler factor of emitting region
    dopfac::Float64 = 0.
    dopnum::Float64 = 0.				        # Doppler boosting exponent, depends on geometry
    counterjet::Bool = false			# boolean switch if user wants to include counterjet emission
    geometry::String = ""   	        # string to track geometry of emitting region

    seed_size::Int					# size of seed photon field
    Niter::Int						# number of IC iterations
    tau::Float64                    # optical depth/comtpon Y of emitting region
    ypar::Float64
    rphot::Float64 					# photospheric radius when tau > 1, used to renormalize volume
    escape_corr::Float64				# escape term, used to renormalize our spectra to CompPS

    seed_energ::Float64				# array of seed frequencies in Hz
    seed_urad::Float64 				# array of seed photon number density in log10(#/erg/cm^3)
    iter_urad::Float64 				# array of iterated photon number density in log10(#/erg/cm^3)

end

const Te_table = @SVector [19.999,30.,50.,70.,90.,120.,150.,200.,
                  280.,400.,600.,900.,1300.,1800.,2500.001]
const tau_table	= @SVector [0.04999,0.07,0.1,0.19,0.27,0.38,0.54,0.76,
                   0.94,1.08,1.34,1.75,2.1,2.6,3.0001]

esc_table_sph = SMatrix{15, 15}(
                    0.02,0.07,0.17,0.27,0.33,0.4,0.45,0.45,0.55,0.55,0.54,0.51,0.48,0.55,0.55,
                    0.02,0.07,0.18,0.275,0.33,0.4,0.45,0.47,0.55,0.53,0.56,0.51,0.5,0.55,0.55,
                    0.02,0.07,0.18,0.275,0.34,0.4,0.45,0.47,0.55,0.54,0.52,0.51,0.53,0.55,0.58,
                    0.025,0.085,0.18,0.275,0.35,0.4,0.45,0.47,0.55,0.54,0.59,0.57,0.5,0.55,0.6,
                    0.028,0.085,0.2,0.275,0.35,0.4,0.45,0.49,0.54,0.54,0.6,0.59,0.53,0.55,0.65,
                    0.035,0.09,0.21,0.29,0.35,0.4,0.45,0.49,0.55,0.57,0.58,0.6,0.53,0.55,0.7,
                    0.045,0.12,0.23,0.3,0.35,0.41,0.42,0.49,0.54,0.57,0.58,0.58,0.5,0.55,0.7,
                    0.055,0.135,0.23,0.3,0.35,0.4,0.42,0.49,0.52,0.57,0.59,0.61,0.58,0.57,0.85,
                    0.07,0.14,0.235,0.3,0.35,0.4,0.41,0.49,0.5,0.57,0.58,0.61,0.61,0.57,0.85,
                    0.08,0.145,0.235,0.3,0.35,0.39,0.4,0.49,0.5,0.57,0.58,0.63,0.7,0.6,0.85,
                    0.09,0.15,0.245,0.3,0.35,0.4,0.4,0.45,0.49,0.55,0.62,0.64,0.73,0.6,0.85,
                    0.105,0.17,0.26,0.31,0.35,0.4,0.4,0.45,0.47,0.56,0.62,0.65,0.78,0.65,0.85,
                    0.115,0.19,0.28,0.32,0.35,0.4,0.4,0.45,0.44,0.5,0.67,0.68,0.8,0.7,1.,
                    0.14,0.215,0.3,0.34,0.38,0.4,0.4,0.45,0.42,0.48,0.68,0.73,0.85,0.8,1.,
                    0.16,0.23,0.31,0.35,0.38,0.4,0.4,0.45,0.42,0.48,0.67,0.78,0.9,0.9,1.
                    )

esc_table_cyl = SMatrix{15, 15}(
                    0.02,0.06,0.22,0.3,0.5,0.5,0.57,0.7,0.75,0.75,0.75,0.77,0.75,0.76,0.7,
                    0.02,0.06,0.19,0.28,0.4,0.5,0.57,0.65,0.65,0.75,0.7,0.75,0.7,0.7,0.75,
                    0.017,0.04,0.17,0.27,0.32,0.5,0.57,0.6,0.6,0.7,0.65,0.65,0.73,0.67,0.7,
                    0.01,0.04,0.12,0.17,0.23,0.35,0.4,0.55,0.5,0.55,0.55,0.65,0.7,0.67,0.68,
                    0.01,0.04,0.1,0.17,0.19,0.28,0.32,0.4,0.45,0.5,0.5,0.5,0.68,0.66,0.68,
                    0.015,0.04,0.09,0.17,0.19,0.25,0.25,0.35,0.45,0.4,0.45,0.5,0.64,0.64,0.68,
                    0.02,0.05,0.09,0.15,0.17,0.22,0.24,0.3,0.35,0.35,0.45,0.4,0.62,0.61,0.64,
                    0.02,0.05,0.095,0.15,0.17,0.22,0.22,0.25,0.3,0.35,0.4,0.38,0.4,0.4,0.6,
                    0.024,0.05,0.095,0.14,0.16,0.21,0.21,0.25,0.27,0.32,0.4,0.42,0.4,0.4,0.5,
                    0.025,0.055,0.1,0.15,0.17,0.21,0.21,0.25,0.27,0.3,0.35,0.4,0.35,0.4,0.43,
                    0.03,0.065,0.11,0.15,0.18,0.21,0.21,0.25,0.27,0.3,0.35,0.4,0.4,0.45,0.47,
                    0.045,0.085,0.135,0.17,0.19,0.22,0.24,0.27,0.3,0.35,0.4,0.45,0.5,0.5,0.55,
                    0.06,0.1,0.16,0.2,0.22,0.25,0.26,0.29,0.33,0.35,0.45,0.5,0.5,0.55,0.65,
                    0.09,0.14,0.2,0.23,0.26,0.29,0.29,0.33,0.4,0.4,0.48,0.55,0.6,0.65,0.75,
                    0.11,0.165,0.23,0.26,0.29,0.31,0.33,0.36,0.4,0.45,0.52,0.6,0.65,0.75,0.8
                    )

function compton(s1::Int, s2::Int)

    size = s1
    seed_size = s2

    en_phot = zeros(Float64, size)
    num_phot = similar(en_phot)
    iter_urad = similar(en_phot)
    en_phot_obs = zeros(Float64, 2*size)
    num_phot_obs = similar(en_phot_obs)

    seed_energ = zeros(Float64, seed_size)
    seed_urad = zeros(Float64, seed_size)

    Niter = 20
    ypad = 0
    escape_corr = 1.0

    counterjet = false

    # xsize = 15
    # ysize = 15

    esc_p_sph = CubicSplineInterpolation((Te_table, tau_table),esc_table_sph)
    esc_p_cyl = CubicSplineInterpolation((Te_table, tau_table),esc_table_cyl)

    return Compton(size=size, seed_size=seed_size, en_phot=en_phot,
                   num_phot=num_phot, iter_urad=iter_urad, en_phot_obs=en_phot_obs,
                   num_phot_obs=num_phot_obs, seed_energ=seed_energ,
                   seed_urad=seed_urad, Niter=Niter, ypar=ypar,
                   escape_corr=escape_corr, esc_p_cyl=esc_p_cyl, esc_p_sph=esc_p_sph)
end

# This function is the kernel of eq 2.48 in Blumenthal & Gould(1970), represents the scattered photon spectrum
# for a given electron and includes the Klein-Nishina cross section.
function comfnc(ein, p::RD.comfnc_params)
    game = p.game
    e1 = p.e1

    phodis = p.phodis

    einit = exp(ein)
    btst = einit/(game*emerg)
    eg4 = 3.0*einit*game
    utst = eg4*(emerg + eg4)

    if ((e1 < btst && (btst-e1)/btst >= 4.0e-8) ||
        (e1 > utst && (e1-utst)/utst >= 4.0e-8))
        return 0.0
    else
        biggam = eg4/emerg
        q = e1/(biggam*(1.0 - e1))
        phonum = phodis(einit)
        tm1 = 2.0*q*log(q)
        tm2 = (1.0 + 2.0*q)*(1.0 - q)
        tm3 = 0.5*((biggam*q)^2)*(1.0 - q)/(1.0 + biggam*q)

        return (tm1 + tm2 + tm3)*10^phonum
    end

end

# This function is the integral of comfnc above over the total seed photon distribution
function comint(gam, p::comint_params)

    eph = p.eph
    ephmin = p.ephmin
    ephmax = p.ephmax

    eldis = p.eldis
    phodis = p.phodis

    game = exp(gam)
    e1 = eph/(game*emerg)
    econst = 2π*re0^2*cee
    blim = log(max(eph/(4.0*game*(game-eph/emerg)), ephmin))
    ulim = log(min(eph, ephmax))

    (ulim <= blim) && return 0.0

    F2params = RD.comfnc_params(game, e1, phodis)
    F2func(x) = comfnc(x, F2params)

    result, error = quadgk(F2func, blim, ulim, rtol=1.0, atol=1.0, maxevals=100)

    den = eldis(game)

    return econst*den*result/game
end

"""This integrates the individual electron spectrum from
comint over the total electron distribution"""
function comintegral(com::Compton, it, blim, ulim, enphot, enphmin, enphmax,
                     eldis)

    F1params = RD.comint_params(enphot, enphmin, enphmax,
                                eldis, com.seed_ph)
    F1params_it = RD.comint_params(enphot, enphmin, enphmax,
                                   eldis, com.iter_ph)

    if iszero(it)
        params = F1params
    else
        params = F1params_it
    end

    F1func(gam) = comint(gam, params)

    result, errror = quadgk(F1func, blim, ulim, rtol=1.0, atol=1.0, maxevals=100)
    return result
end

"""
This function calculates the final spectrum in all frequency bins,
including multiple scatters.
Note: the reason the Doppler boosting is a factor of 2 instead of 3 is because the calculations are done for
a conical jet in the Lind&BLanford 1985 prescription
"""
function compton_spectrum(com::Compton, gmin, gmax, eldis)

    ephmin = com.seed_energ[1]
    ephmax = com.seed_energ[end]

    dopfac_cj = com.dopfac*(1.0 - com.beta*cos(com.angle))/
                (1.0 +  com.beta*cos(com.angle))

    for it = 1:com.Niter
        for i = 1:com.size
            blim = log(max(gmin, com.en_phot[i]/emerg))
            ulim = log(gmax)

            if blim >= ulim
                comres = 1e-100
            else
                comres = comintegral(com, it, blim, ulim, com.en_phot[i],
                                  ephmin, ephmax, eldis)
            end

            com.num_phot[i] = com.num_phot[i]+comres*com.vol*com.en_phot[i]*herg
            com.en_phot_obs[i] = com.en_phot[i]*com.dopfac
            com.num_phot_obs[i] = com.num_phot[i]*dopfac_cj^com.dopnum

            if com.counterjet
                com.en_phot_obs[i+com.size] = com.en_phot[i]*dopfac_cj
                com.num_phot_obs[i+com.size] = com.num_phot[i]*dopfac_cj^com.dopnum
            else
                com.en_phot_obs[i+com.size] = 0.0
                com.num_phot_obs[i+com.size] = 0.0
            end

            if iszero(comres)
                com.iter_urad[i] = -50.0
            else
                com.iter_urad[i] = log10(com.escape_corr*comres*com.vol)/(π*com.r^2*cee)
            end

        end

        ephmin = com.en_phot[1]
        ephmax = com.en_phot[end]
        com.iter_ph = LinearInterpolation(com.en_phot, com.iter_urad)
    end
end

"""Method to include cyclosynchrotron array from `Cyclosyn.jl` to the seed field. The two input arrays should be
`get_energ()` and `get_nphot()` methods from `Cyclosyn.jl`. If the seed_urad array wasn't empty, the contribution
is automatically added to the existing one.
note: the second condition has a <= sign to dodge numerical errors when the seed field energy density is
extremely low, which can result in negative values/nan for the seed_urad """
function cyclosyn_seed(com::Compton, seed_arr, seed_lum)
    chr2 = cee*herg*π*com.r^2
    for i in eachindex(seed_arr)
        com.seed_energ[i] = seed_arr[i]

        if !iszero(com.seed_urad[i])
            com.seed_urad[i] = log10(10^com.seed_urad[i] + seed_lum[i]/(chr2*com.seed_energ[i]))
        elseif seed_lum[i]/(chr2*com.seed_energ[i]) <= 0
            com.seed_urad[i] = -100
        else
            com.seed_urad[i] = log10(seed_lum[i]/(chr2*com.seed_energ[i]))
        end

    end

    com.seed_ph = LinearInterpolation(seed_energ, seed_urad)

end

"""
Method to include black body to seed field for IC;
Note: Urad and Tbb need to be passed in the co-moving frame, the function does NOT account for beaming
"""
function bb_seed_k(com::Compton, seed_arr, Urad, Tbb)

    ulim = 100.0*Tbb*kboltz
    bfield_const = herg*cee^2*sbconst*Tbb^4

    seed_freq_array(com, seed_arr)
    for i = 1:com.seed_size
        if (com.seed_energ[i] < ulim)
            bbfield = (2.0*Urad*(seed_energ[i]/herg)^2.)/(bfield_const*(exp(seed_energ[i]/(kboltz*Tbb)) - 1.0))
        else
            bbfield = 1.0e-100
        end

        if !iszero(com.seed_urad[i])
            com.seed_urad[i] = log10(10^com.seed_urad[i] + bbfield)
        elseif bbfield > 0
            com.seed_urad[i] = log10(bbfield)
        else
            com.seed_urad[i] = -100.0
        end

    end

    com.seed_ph = LinearInterpolation(com.seed_energ, com.seed_urad)

end

function bb_seed_kev(com::Compton, seed_arr, Urad, Tbb)

    ulim = 100.0*Tbb*kboltz_kev2erg
    bfield_const = herg*cee^2*sbconst*(Tbb*kboltz_kev2erg/kboltz)^4

    seed_freq_array(com, seed_arr)
    for i = 1:com.seed_size
        if (com.seed_energ[i] < ulim)
            bbfield = (2.0*Urad*(seed_energ[i]/herg)^2.)/(bfield_const*(exp(seed_energ[i]/(kboltz*Tbb)) - 1.0))
        else
            bbfield = 1.0e-100
        end

        if !iszero(com.seed_urad[i])
            com.seed_urad[i] = log10(10^com.seed_urad[i] + bbfield)
        elseif bbfield > 0
            com.seed_urad[i] = log10(bbfield)
        else
            com.seed_urad[i] = -100.0
        end

    end

    com.seed_ph = LinearInterpolation(com.seed_energ, com.seed_urad)

end

"""
Calculates incident photon field for a Shakura-Sunyaev disk with aspect ratio H, inner temperature Tin,
inner radius Rin, at a distance z from the disk, taking beaming into account
"""
function disk_integral(alfa, p)

end

function seed_freq_array(com::Compton, seed_arr)
    com.seed_energ .= seed_arr
end


text:
using ObjectDetector, FileIO
using Profile

Profile.init(n=10000000, delay=0.01)
yolomod = YOLO.v3_416_COCO()

batch = emptybatch(yolomod)
img = load(joinpath(dirname(dirname(pathof(ObjectDetector))),"test","images","dog-cycle-car.png"))

batch[:,:,:,1] .= gpu(resizePadImage(img, yolomod)) # Send resized image to the batch
res = yolomod(batch) # Run the model on the length-1 batch

imgBoxes = drawBoxes(img, res)

save(joinpath(@__DIR__,"test.png"), imgBoxes)

ObjectDetector.benchmark()

save(joinpath(@__DIR__,"test.png"), resizePadImage(img, yolomod))


text:
# Using mutable to allow missing values
@missable mutable struct Player
    name::String
    rating::String
end

@missable mutable struct Stats
    hp::Int
    attack::Int
    defense::Int
    special_attack::Int
    special_defense::Int
    speed::Int
end

@missable mutable struct Pokemon
    id::PokemonID
    name::String
    moves::Vector{MoveID} = MoveID[]
    gender::Gender
    level::Int64
    item::ItemID
    stats::Stats
    ability::AbilityID
    nature::Nature
    sent_out_number::Int  # number of sent out (1/2 for leading pokemons)
    shiny::Yes
    dynamaxed::Yes
end

@missable mutable struct Team
    title::String
    author::String
    notes::String
    rental_code::String
    pokemons::Vector{Pokemon} = Pokemon[]
end

function update_team!(team::Team, i, field, value)
    pokes = team.pokemons
    if length(pokes) ≡ 0
        for _ in 1:6
            push!(pokes, Pokemon())
        end
    end
    poke = pokes[i]
    if field ≡ :moves
        moves = poke.moves
        if !contains(moves, value)
            push!(moves, value)
        end
    elseif field ∈ (:hp, :attack, :defense, :special_attack, :special_defense, :speed)
        if ismissing(poke.stats)
            poke.stats = Stats()
        end
        setfield!(poke.stats, field, value)
    else
        setfield!(poke, field, value)
    end
    team
end


text: Dict{String, Any}("end_esc" => Dict{String, Any}("value" => "String does not end here\" but ends here\\", "type" => "string"), "multiline_end_esc" => Dict{String, Any}("value" => "When will it end? \"\"\"...\"\"\" should be here\"", "type" => "string"), "lit_multiline_not_unicode" => Dict{String, Any}("value" => "\\u007f", "type" => "string"), "multiline_unicode" => Dict{String, Any}("value" => " ", "type" => "string"), "lit_end_esc" => Dict{String, Any}("value" => "String ends here\\", "type" => "string"), "multiline_not_unicode" => Dict{String, Any}("value" => "\\u0041", "type" => "string"), "lit_multiline_end" => Dict{String, Any}("value" => "There is no escape\\", "type" => "string"))


text:
# Unscented Transform

# Test constructors
ut_params = UnscentedTransform()
@test ut_params == UnscentedTransform(1e-3, 2, 0)

# Test linear system augment methods
sys = LinearSystem(eye(2), 2.25*eye(2))
dstate = UncertainDiscreteState(1.5*ones(2), 3.0*eye(2))
astate, asys = augment(dstate, sys)
@test (length(astate.x) == 4) & (size(astate.P) == (4,4)) & (size(asys.A) == (4,4))

# Test consistency between linear prediction and UnscentedTransform
nominal_prediction = predict(dstate, sys, 1)
astate, asys = augment(dstate, sys)
disc_fcn(x) = predict(x, asys, 1)
unscented_prediction = transform(astate, disc_fcn, ut_params)
@test isapprox(nominal_prediction.x, unscented_prediction.x[1:2])
@test isapprox(nominal_prediction.P, unscented_prediction.P[1:2,1:2])

# Test nonlinear system augment methods
F(t,x) = x.^2
dF(t,x) = diagm(2*x)
nsys = NonlinearSystem(F, dF, 2*eye(2))
dstate = UncertainDiscreteState(1.5*ones(2), 3.0*eye(2))
astate, asys = augment(dstate, nsys)
@test length(astate.x) == 4

# Test linear observer augment methods
obs = LinearObserver([1.0, 2.0]', 1.0)
dstate = UncertainDiscreteState(1.5*ones(2), 3.0*eye(2))
astate, aobs = augment(dstate, obs)
@test (length(astate.x) == 3) & (size(astate.P) == (3,3)) & (size(aobs.H) == (1,3))

# Test nonlinear observer augment methods
H(t,x) = [x[1].^2]
dH(t,x) = [2.0*x[1] 0.0]
nobs = NonlinearObserver(H, dH, 2.0)
dstate = UncertainDiscreteState(1.5*ones(2), 3.0*eye(2))
astate, aobs = augment(dstate, nobs)
@test length(astate.x) == 3

# TODO: Test combined augment methods
augment(dstate, sys, obs)
augment(dstate, nsys, obs)
augment(dstate, sys, nobs)
augment(dstate, nsys, nobs)

# Test constructors
UnscentedKalmanFilter(nsys, nobs, dstate)
UnscentedKalmanFilter(nsys, nobs, dstate, UnscentedTransform())
UnscentedKalmanFilter(nsys, nobs, dstate, [1])

# Test linear observer prediction
obs = LinearObserver(eye(2), 2*eye(2))
ut_state = predict(dstate, obs, ut_params)
expected_state = predict(dstate, obs)
@test isapprox(ut_state.x, expected_state.x)
@test isapprox(ut_state.P, expected_state.P)

# Test linear system prediction
sys = LinearSystem(eye(2), 2*eye(2))
ut_state = predict(dstate, sys, ut_params, 2)
expected_state = predict(dstate, sys, 2)
@test isapprox(ut_state.x, expected_state.x)
@test isapprox(ut_state.P, expected_state.P)

# Test linear ukf
srand(1);
linear_sys = LinearSystem(0.5*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainDiscreteState([1.0, 2.0], 0.1*eye(2))
ukf = UnscentedKalmanFilter(linear_sys, linear_obs, initial_est)
kf = KalmanFilter(linear_sys, linear_obs, initial_est)
simulator = make_simulator(kf)
for i = 1:4
    true_state, measurement = simulate(simulator, i)
    process!(ukf, measurement)
    process!(kf, measurement)
    @test isapprox(ukf.estimate.x, kf.estimate.x)
end

srand(1);
linear_sys = LinearSystem(0.5*eye(2), 0.001*eye(2))
linear_obs = LinearObserver(eye(2), 0.001*eye(2))
initial_est = UncertainDiscreteState([1.0, 2.0], 0.1*eye(2))
ukf = UnscentedKalmanFilter(
    NonlinearSystem(linear_sys), NonlinearObserver(linear_obs), initial_est)
kf = KalmanFilter(linear_sys, linear_obs, initial_est)
simulator = make_simulator(kf)
for i = 1:4
    true_state, measurement = simulate(simulator, i)
    process!(ukf, measurement)
    process!(kf, measurement)
    @test isapprox(ukf.estimate.x, kf.estimate.x)
end

# Test nonlinear UKF with linear observer
srand(1);
F(t,x) = [sin(x[2]), -x[1]]
dF_dx(t,x) = [0 cos(x[2]); -1 0]
nonlinear_sys = NonlinearSystem(F, dF_dx, 0.01*eye(2,2))
linear_obs = LinearObserver(eye(2), 0.1*eye(2))
initial_est = UncertainDiscreteState([0.0, 2.0], 0.1*eye(2))
ukf = UnscentedKalmanFilter(nonlinear_sys, linear_obs, initial_est);
simulator = make_simulator(ukf)
for idx = 1:10
    true_state, measurement = simulate(simulator, idx)
    process!(ukf, measurement)
    @test mahalanobis(true_state, ukf.estimate) < 3
end

# Test convenience function
srand(1);
F(t,x) = [sin(x[2]), -x[1]]
dF_dx(t,x) = [0 cos(x[2]); -1 0]
nonlinear_sys = NonlinearSystem(F, dF_dx, 0.01*eye(2,2))
linear_obs = LinearObserver(eye(2), 0.1*eye(2))
estimate = UncertainDiscreteState([0.0, 2.0], 0.1*eye(2))
ut_params = UnscentedTransform()
ukf = UnscentedKalmanFilter(nonlinear_sys, linear_obs, estimate)
simulator = make_simulator(ukf)
for idx = 1:10
    true_state, measurement = simulate(simulator, idx)
    estimate = unscented_kalman_filter(
        estimate, nonlinear_sys, linear_obs, ut_params, measurement)
end

# Test convenience function with archiving
simulator = make_simulator(ukf)
archive = EstimatorHistory(ukf)
for idx = 1:10
    true_state, measurement = simulate(simulator, idx)
    estimate = unscented_kalman_filter(
        estimate, nonlinear_sys, linear_obs, ut_params, measurement, archive)
end


text:
using StatisticalRethinking, DataFrames, NamedArrays, OrderedCollections

na = NamedArray(
    [18 134 152; 3 4 7; 21 138 159],
    (   [:vaccinated, :unvaccinated, :total],
        [:varicella, :noninfected, :total]
    ),
    ("Rows", "Cols")
)

na |> display
println()

# Risk of varicella among vaccinated children = 18 / 152 = 0.118 = 11.8%

risk_vac = na[:vaccinated, :varicella] / na[:vaccinated, :total]

# Risk of varicella among unvaccinated children = 3 / 7 = 0.429 = 42.9%

risk_unvac = na[:unvaccinated, :varicella] / na[:unvaccinated, :total]

# Risk ratio = 0.118 / 0.429 = 0.28

rr = risk_vac / risk_unvac

# EXAMPLE: Calculating Vaccine Effectiveness
# Calculate the vaccine effectiveness from the varicella data in Table 3.13.

# VE = (42.9 – 11.8) / 42.9 = 31.1 / 42.9 = 72%

ve = (risk_unvac - risk_vac) /risk_unvac

# Alternatively, VE = 1 – RR = 1 – 0.28 = 72%

1 - rr |> display

# So, the vaccinated group experienced 72% fewer varicella cases
# than they would have if they had not been vaccinated.

df_covid = DataFrame(
    :completed => [true, true, false, false],
    :company => [:pfizer, :moderna, :jj, :astra],
    :unvaccinated => [18261, 15000, missing, missing],
    :vaccinated => [18262, 15000, missing, missing],
    :inf_unv => [162, 185, missing, missing],
    :inf_vac => [8, 11, missing, missing]
)

df_covid |> display

df = df_covid[df_covid.completed, :]


text:
module JungleHelperRuinsSpikesUp

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Spikes (Up, Ruins) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.SpikesUp,
        "rectangle",
        Dict{String, Any}(
          "type" => "JungleHelper/ruins"
        )
    ),
    "Spikes (Down, Ruins) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.SpikesDown,
        "rectangle",
        Dict{String, Any}(
          "type" => "JungleHelper/ruins"
        )
    ),
	"Spikes (Left, Ruins) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.SpikesLeft,
        "rectangle",
        Dict{String, Any}(
          "type" => "JungleHelper/ruins"
        )
    ),
	"Spikes (Right, Ruins) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.SpikesRight,
        "rectangle",
        Dict{String, Any}(
          "type" => "JungleHelper/ruins"
        )
    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper