text:
using IndividualDisplacements
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/flow_fields.jl"));

"""
    init_global_randn(np ::Int , 𝑃::NamedTuple)

Randomly distribute `np` points over the Earth, within `𝑃.msk`
region, and return position in grid index space (`i,j,subdomain`).
"""
function init_global_randn(np ::Int , 𝑃::NamedTuple)
    (lon, lat) = randn_lonlat(maximum([2*np 10]))
    (_,_,_,_,f,x,y)=InterpolationFactors(𝐷.Γ,lon,lat)
    m=findall( (f.!==0).*((!isnan).(x)) )
    n=findall(nearest_to_xy(𝐷.msk,x[m],y[m],f[m]).==1.0)[1:np]
    return permutedims([x[m[n]] y[m[n]] f[m[n]]])
end

"""
    reset_📌!(𝐼::Individuals,frac::Number,📌::Array)

Randomly select a fraction (frac) of the particles and reset
their positions (𝐼.📌) to a random subset of the specificed 📌.
"""
function reset_📌!(𝐼::Individuals,frac::Number,📌::Array)
    np=length(𝐼.🆔)
    n_reset = Int(round(𝐷.frac*np))
    k_reset = rand(1:np, n_reset)
    l_reset = rand(1:np, n_reset)
    𝐼.📌[k_reset]=deepcopy(📌[l_reset])
    isempty(𝐼.🔴.ID) ? m=maximum(𝐼.🆔) : m=max(maximum(𝐼.🔴.ID),maximum(𝐼.🆔))
    𝐼.🆔[k_reset]=collect(1:n_reset) .+ m
end

##

"""
    isosurface(θ,T,z)

```
isosurface(𝐼.𝑃.θ0,15,Γ["RC"])
```
"""
function isosurface(θ,T,z)
    d=NaN*similar(θ[:,1])
    nr=size(θ,2)
    for j=1:size(d,1)
        for k=1:nr-1
            i=findall(isnan.(d[j]).&(θ[j,k].>T).&(θ[j,k+1].<=T))
            a=(θ[j,k][i] .- T)./(θ[j,k][i] .- θ[j,k+1][i])
            d[j][i]=(1 .- a).*Γ["RC"][k] + a.*Γ["RC"][k+1]
            i=findall(isnan.(d[j]).&(θ[j,k].<=T).&(θ[j,k+1].>T))
            a=(θ[j,k+1][i] .- T)./(θ[j,k+1][i] .- θ[j,k][i])
            d[j][i]=(1 .- a).*Γ["RC"][k+1] + a.*Γ["RC"][k]
        end
    end
    return d
end

"""
    set_up_𝑃(k::Int,t::Float64,Γ::Dict,pth::String)

Define `FlowFields` data structure (𝑃) along with ancillary variables (𝐷)
for the specified grid (`Γ` dictionnary), vertical level (`k`), and
file location (`pth`).

_Note: the initial implementation approximates month durations to
365 days / 12 months for simplicity and sets 𝑃.𝑇 to [-mon/2,mon/2]_
"""
function set_up_FlowFields(k::Int,Γ::Dict,pth::String)
    XC=exchange(Γ["XC"]) #add 1 lon point at each edge
    YC=exchange(Γ["YC"]) #add 1 lat point at each edge
    iDXC=1. ./Γ["DXC"]
    iDYC=1. ./Γ["DYC"]
    γ=Γ["XC"].grid
    mon=86400.0*365.0/12.0
    func=Γ["update_location!"]

    if k==0
        msk=Γ["hFacC"]
        (_,nr)=size(msk)
        𝑃=FlowFields(MeshArray(γ,Float64,nr),MeshArray(γ,Float64,nr),
        MeshArray(γ,Float64,nr),MeshArray(γ,Float64,nr),
        MeshArray(γ,Float64,nr+1),MeshArray(γ,Float64,nr+1),
        [-mon/2,mon/2],func)
    else
        msk=Γ["hFacC"][:, k]
        𝑃=FlowFields(MeshArray(γ,Float64),MeshArray(γ,Float64),
        MeshArray(γ,Float64),MeshArray(γ,Float64),[-mon/2,mon/2],func)
    end

    𝐷 = (🔄 = update_FlowFields!, pth=pth,
         XC=XC, YC=YC, iDXC=iDXC, iDYC=iDYC,
         k=k, msk=msk, θ0=similar(msk), θ1=similar(msk))

    tmp = IndividualDisplacements.dict_to_nt(IndividualDisplacements.NeighborTileIndices_cs(Γ))
    𝐷 = merge(𝐷 , tmp)

    return 𝑃,𝐷
end

"""
    update_FlowFields!(𝑃::𝐹_MeshArray2D,𝐷::NamedTuple,t::Float64)

Update flow field arrays (in 𝑃), 𝑃.𝑇, and ancillary variables (in 𝐷)
according to the chosen time `t` (in `seconds`).

_Note: for now, it is assumed that (1) the time interval `dt` between
consecutive records is diff(𝑃.𝑇), (2) monthly climatologies are used
with a periodicity of 12 months, (3) vertical 𝑃.k is selected_
"""
function update_FlowFields!(𝑃::𝐹_MeshArray2D,𝐷::NamedTuple,t::Float64)
    dt=𝑃.𝑇[2]-𝑃.𝑇[1]

    m0=Int(floor((t+dt/2.0)/dt))
    m1=m0+1
    t0=m0*dt-dt/2.0
    t1=m1*dt-dt/2.0

    m0=mod(m0,12)
    m0==0 ? m0=12 : nothing
    m1=mod(m1,12)
    m1==0 ? m1=12 : nothing

    (U,V)=read_velocities(𝑃.u0.grid,m0,𝐷.pth)
    u0=U[:,𝐷.k]; v0=V[:,𝐷.k]
    u0[findall(isnan.(u0))]=0.0; v0[findall(isnan.(v0))]=0.0 #mask with 0s rather than NaNs
    u0=u0.*𝐷.iDXC; v0=v0.*𝐷.iDYC; #normalize to grid units
    (u0,v0)=exchange(u0,v0,1) #add 1 point at each edge for u and v

    (U,V)=read_velocities(𝑃.u0.grid,m1,𝐷.pth)
    u1=U[:,𝐷.k]; v1=V[:,𝐷.k]
    u1[findall(isnan.(u1))]=0.0; v1[findall(isnan.(v1))]=0.0 #mask with 0s rather than NaNs
    u1=u1.*𝐷.iDXC; v1=v1.*𝐷.iDYC; #normalize to grid units
    (u1,v1)=exchange(u1,v1,1) #add 1 point at each edge for u and v

    𝑃.u0[:]=u0[:]
    𝑃.u1[:]=u1[:]
    𝑃.v0[:]=v0[:]
    𝑃.v1[:]=v1[:]
    𝑃.𝑇[:]=[t0,t1]

end

"""
    update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)

Update flow field arrays (in 𝑃), 𝑃.𝑇, and ancillary variables (in 𝐷)
according to the chosen time `t` (in `seconds`).

_Note: for now, it is assumed that (1) the time interval `dt` between
consecutive records is diff(𝑃.𝑇), (2) monthly climatologies are used
with a periodicity of 12 months, (3) vertical 𝑃.k is selected_
"""
function update_FlowFields!(𝑃::𝐹_MeshArray3D,𝐷::NamedTuple,t::Float64)
    dt=𝑃.𝑇[2]-𝑃.𝑇[1]

    m0=Int(floor((t+dt/2.0)/dt))
    m1=m0+1
    t0=m0*dt-dt/2.0
    t1=m1*dt-dt/2.0

    m0=mod(m0,12)
    m0==0 ? m0=12 : nothing
    m1=mod(m1,12)
    m1==0 ? m1=12 : nothing

    (_,nr)=size(𝐷.Γ["hFacC"])

    (U,V)=read_velocities(𝑃.u0.grid,m0,𝐷.pth)
    u0=U; v0=V
    u0[findall(isnan.(u0))]=0.0; v0[findall(isnan.(v0))]=0.0 #mask with 0s rather than NaNs
    for k=1:nr
        u0[:,k]=u0[:,k].*𝐷.iDXC; v0[:,k]=v0[:,k].*𝐷.iDYC; #normalize to grid units
        (tmpu,tmpv)=exchange(u0[:,k],v0[:,k],1) #add 1 point at each edge for u and v
        u0[:,k]=tmpu
        v0[:,k]=tmpv
    end
    w0=IndividualDisplacements.read_nctiles(𝐷.pth*"WVELMASS/WVELMASS","WVELMASS",𝑃.u0.grid,I=(:,:,:,m0))
    w0[findall(isnan.(w0))]=0.0 #mask with 0s rather than NaNs

    (U,V)=read_velocities(𝑃.u0.grid,m1,𝐷.pth)
    u1=U; v1=V
    u1[findall(isnan.(u1))]=0.0; v1[findall(isnan.(v1))]=0.0 #mask with 0s rather than NaNs
    for k=1:nr
        u1[:,k]=u1[:,k].*𝐷.iDXC; v1[:,k]=v1[:,k].*𝐷.iDYC; #normalize to grid units
        (tmpu,tmpv)=exchange(u1[:,k],v1[:,k],1) #add 1 point at each edge for u and v
        u1[:,k]=tmpu
        v1[:,k]=tmpv
    end
    w1=IndividualDisplacements.read_nctiles(𝐷.pth*"WVELMASS/WVELMASS","WVELMASS",𝑃.u0.grid,I=(:,:,:,m1))
    w1[findall(isnan.(w1))]=0.0 #mask with 0s rather than NaNs

    𝑃.u0[:,:]=u0[:,:]
    𝑃.u1[:,:]=u1[:,:]
    𝑃.v0[:,:]=v0[:,:]
    𝑃.v1[:,:]=v1[:,:]
    for k=1:nr
        tmpw=exchange(-w0[:,k],1)
        𝑃.w0[:,k]=tmpw./𝐷.Γ["DRC"][k]
        tmpw=exchange(-w1[:,k],1)
        𝑃.w1[:,k]=tmpw./𝐷.Γ["DRC"][k]
    end
    𝑃.w0[:,1]=0*exchange(-w0[:,1],1)
    𝑃.w1[:,1]=0*exchange(-w1[:,1],1)
    𝑃.w0[:,nr+1]=0*exchange(-w0[:,1],1)
    𝑃.w1[:,nr+1]=0*exchange(-w1[:,1],1)

    θ0=IndividualDisplacements.read_nctiles(𝐷.pth*"THETA/THETA","THETA",𝑃.u0.grid,I=(:,:,:,m0))
    θ0[findall(isnan.(θ0))]=0.0 #mask with 0s rather than NaNs
    𝐷.θ0[:,:]=θ0[:,:]

    θ1=IndividualDisplacements.read_nctiles(𝐷.pth*"THETA/THETA","THETA",𝑃.u0.grid,I=(:,:,:,m1))
    θ1[findall(isnan.(θ1))]=0.0 #mask with 0s rather than NaNs
    𝐷.θ1[:,:]=θ1[:,:]

    𝑃.𝑇[:]=[t0,t1]
end

ext: jl


text:
# create type MatrixNetwork
type MatrixNetwork{T}
    n::Int64 # number of columns/rows
    rp::Vector{Int64} # row pointers
    ci::Vector{Int64} # column indices
    vals::Vector{T} # corresponding values
end


text:
# =========================================================================== #
## Model
# =========================================================================== #

abstract type AbstractModel end
abstract type AbstractSubModel <: AbstractModel end

mutable struct Model <: AbstractModel
    # add submodel vector
    nodes::Vector{AbstractBlock}
    names::Vector{String}
    signals::Vector{Signal}
    map_node_to_index::Dict{AbstractBlock, Int}
    graph::datastructs.MatrixGraph
    # incidence_matrix::Matrix{Int}

    function Model()
        nodes   = Vector{AbstractBlock}(undef, 0)
        names   = Vector{String}(undef, 0)
        signals = Vector{Signal}(undef, 0)

        map_node_to_index = Dict{AbstractBlock, Int}()

        return new(nodes, names, signals, map_node_to_index)
    end
end

function Base.show(io::IO, model::Model)
    println("Model with properties:")
    println("\t   n nodes: $(length(model.nodes))")
    println("\t   n names: $(length(model.names))")
    println("\t n signals: $(length(model.signals))")

end

function Base.push!(model::Model, node::N, name::S) where {N <: AbstractBlock, S <: AbstractString}

    last_idx = length(model.nodes)
    new_idx  = last_idx + 1
    model.map_node_to_index[node] = new_idx

    push!(model.nodes, node)
    push!(model.names, name)

    return model
end

# =========================================================================== #
## Connection
# =========================================================================== #

function connect!(model::Model, from_pin::OutType, to_pin::InType) where {OutType <: AbstractOutputPin, InType <: AbstractInputPin}
# tenedyn.connect!(model, sensor.y, ctrl_minus.x2)

    signal = Signal(from_pin, to_pin)
    push!(model.signals, signal)

    if has_child(from_pin)
        child_signal = from_pin.left_child_signal

        while has_right_sibling(child_signal)
            child_signal = child_signal.right_sibling_signal
        end
        child_signal.right_sibling_signal = signal
    else
        from_pin.left_child_signal = signal
    end

    return model
end

function find_integrators(model::Model)
    L       = length(model.nodes)
    n_init  = ceil(Int, L/5)
    indices = Vector{Int}(undef, n_init)
    jj      = 1

    for ii = 1:length(model.nodes)
            t = typeof( model.nodes[ii] )
            if t <: tenedyn.AbstractIntegrator
                indices[jj] = ii
                jj += 1
            end
    end

    indices = indices[1:jj-1]
    return indices
end

function compile(model::Model)

        # validation
        validate_model(model)

        # recompute the connection graph
        create_graph!(model)

        # initialize the forest
        integrator_indices  = find_integrators(model)
        n_states = length(integrator_indices)
        forest   = initialize_forest(n_states)

        tree_node_vecvec           = Vector{Vector{Any}}(undef, n_states)
        evalulation_indices_vecvec = Vector{Vector{Int}}(undef, n_states)

        # build the trees
        for ii = 1:n_states
            tree = build_tree(model, integrator_indices[ii])
            forest[ii] = tree

            (tree_nodes, tree_depths, tree_leaf, evaluation_indices) =
                tenedyn.extended_depth_first_search(model, tree)

            tree_node_vecvec[ii]           = tree_nodes
            evalulation_indices_vecvec[ii] = evaluation_indices
        end

        return (forest, integrator_indices, tree_node_vecvec, evalulation_indices_vecvec)
end

function validate_model(model::Model)
        # do nothing yet
end

function create_graph!(model::Model)
        # create incidence matrix

        n_nodes = length(model.nodes)
        n_edges = length(model.signals)
        inc_mat = zeros(Int, n_nodes, n_nodes)

        for ii = 1:n_edges
            edge = model.signals[ii]
            start_node = edge.start_pin.containing_object
            end_node   = edge.end_pin.containing_object

            si = find_block_index(model, start_node)
            ei = find_block_index(model, end_node)

            inc_mat[si, ei] = -1
            inc_mat[ei, si] = +1
        end

        model.graph = datastructs.MatrixGraph(inc_mat)

        return model
end

function find_block_index(model::Model, node)
        return model.map_node_to_index[node]
end

function initialize_forest(n_states::Int)
        forest = Vector{datastructs.Tree}(undef, n_states)
        return forest
end

struct TreeSearchHelper
        depth::Int
        visited::Bool
        last_tree_node::datastructs.TreeNode
end

function build_tree(model::Model, node_index::Int)

        # basically bfs/dfs search
        node = model.nodes[node_index]

        tree = datastructs.Tree(node)

        NodeType     = AbstractBlock
        init_size    = div( length(model.nodes), 2)
        helper       = TreeSearchHelper(1, false, tree.node)
        node_queue   = datastructs.Queue{NodeType}(init_size)
        helper_queue = datastructs.Queue{TreeSearchHelper}(init_size)
        push!(node_queue,   node)
        push!(helper_queue, helper)

        is_first = true
        while !isempty(node_queue)

            node   = pop!(node_queue)
            helper = pop!(helper_queue)

            if helper.visited
                # do nothing ?
            else

                # add the same node once again
                # helper = TreeSearchHelper(helper.depth, true, helper.last_tree_node)
                # push!(node_queue,   node)  # probably not necessary
                # push!(helper_queue, helper) # probably not necessary

                # stop search if:
                #   - the node has no input pin
                #   - input pin is not defined == nil
                #   => all of these can be expressed by the incidence matrix
                #   - OR the child node is an integrator

                # add the children to the queue
                if !is_first && is_integrator(node)
                    continue
                end
                is_first = false

                next_depth = helper.depth + 1
                visited    = false

                node_index  = model.map_node_to_index[node]
                adj_indices = datastructs.find_incoming_nodes(model.graph, node_index)

                for ii = 1:length(adj_indices)
                    child_node = model.nodes[ adj_indices[ii] ]
                    tree_node  = datastructs.TreeNode(child_node)
                    datastructs.insert_after!(tree, helper.last_tree_node, tree_node)

                    push!(node_queue,   child_node)
                    push!(helper_queue, TreeSearchHelper(next_depth, visited, tree_node))
                end

            end
        end

        return tree
end

function extended_depth_first_search(model::Model, tree::datastructs.Tree)

        n        = tree.count
        depths   = zeros(Int, n)
        nodes    = Vector{Any}(undef, n)
        is_leaf  = falses(n)
        evaluation_indices = zeros(Int, n)

        init_size    = div(n, 2)
        helper       = datastructs.TreeSearchHelper(1, false)
        node_stack   = datastructs.Stack{Any}(init_size)

        helper_stack = datastructs.Stack{datastructs.TreeSearchHelper}(init_size)
        push!(node_stack,   tree.node)
        push!(helper_stack, helper)

        index = 0
        while !isempty(node_stack)

            node   = pop!(node_stack)
            helper = pop!(helper_stack)

            if helper.visited

                index         += 1
                nodes[index]   = node
                depths[index]  = helper.depth
                is_leaf[index] = datastructs.is_leaf(node)
                evaluation_indices[index] = model.map_node_to_index[node.data]
            else

                # add the same node once again
                helper = datastructs.TreeSearchHelper(helper.depth, true)
                push!(node_stack,   node)
                push!(helper_stack, helper)

                # add the children
                next_depth = helper.depth + 1
                visited    = false

                if !datastructs.has_child(node)
                    continue
                end

                node = node.left_child
                push!(node_stack,   node)
                push!(helper_stack, datastructs.TreeSearchHelper(next_depth, visited))

                while datastructs.has_right_sibling(node)
                    node = node.right_sibling
                    push!(node_stack,   node)
                    push!(helper_stack, datastructs.TreeSearchHelper(next_depth, visited))
                end

            end

        end

        return (nodes, depths, is_leaf, evaluation_indices)
end

function eval_tree(
        model::Model,
        tree::datastructs.Tree,
        tree_nodes, # vector of any -> any treenode
        evaluation_indices::Vector{Int},
        time::Float64)

        value = tree.node.data.y.value
        L     = length(evaluation_indices)

        for ii = 1:L

            node_idx  = evaluation_indices[ii]
            node      = model.nodes[node_idx]
            tree_node = tree_nodes[ii]

            # eval node
            eval(node, time)
            # at the end of this eval the output of the node is updated
            # the output can be up propagated to its parent

            if !datastructs.is_root(tree_node)
                parent_tree_node = tree_node.parent
                parent_node      = parent_tree_node.data
                # parent_node.x?   = node.y.value

                # maybe a hash table for n1.outputpin - n2.inputpin -> int into edges

                # I could ask:
                symbol_list = get_output_pin_names(node)
                for sym in symbol_list
                    out_pin = getfield(node, sym)

                    # search in the outgoing edges for the parent node
                    signal = out_pin.left_child_signal
                    found = false
                    while true
                        potentially_the_parent = signal.end_pin.containing_object
                        if potentially_the_parent == parent_node
                            # propagate up the tree
                            signal.end_pin.value = out_pin.value
                            found = true
                            break
                        else
                            if has_right_sibling(signal)
                                signal = signal.right_sibling_signal
                            else
                                break
                            end
                        end
                    end
                    if !found
                        error("parent node has not been found")
                    end
                    value = signal.end_pin.value
                end

            end # if

        end # for

        return value
end

function get_output_pin_names(node)
        symbol_list = propertynames(node)
        bool        = map( sym -> isa(getfield(node, sym), OutputPin), symbol_list )
        idx         = findall(bool)
        symbol_list = symbol_list[idx]

        return symbol_list
end

function get_input_pin_names(node)
    symbol_list = propertynames(node)
    bool        = map( sym -> isa(getfield(node, sym), InputPin), symbol_list )
    idx         = findall(bool)
    symbol_list = symbol_list[idx]

    return symbol_list
end

function create_ode_function(
        model::Model,
        forest,
        integrator_indices,
        tree_node_vecvec,
        evalulation_indices_vecvec)

        # eval_tree(
        # model::Model,
        # tree::datastructs.Tree,
        # tree_nodes, # vector of any -> any treenode
        # evaluation_indices::Vector{Int},
        # time::Float64)

        n_states = length(integrator_indices)

        ode_fcn(der_q, time, q) = begin

            # update integrators
            for ii = 1:n_states
                model.nodes[ integrator_indices[ii] ].y.value = q[ii]
            end

            # evaluate the trees
            for ii = 1:n_states

                tree               = forest[ii]
                tree_nodes         = tree_node_vecvec[ii]
                evaluation_indices = evalulation_indices_vecvec[ii]

                value = eval_tree(
                    model,
                    tree,
                    tree_nodes,
                    evaluation_indices,
                    time)

                der_q[ii] = value
            end

            return der_q
        end

        return ode_fcn
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Atmos RHS

"""
This is a slightly modified ClimaAtmos Ekman column model
"""

#=
Ekman column:
    ∂_t ρ =  ∇ (μ ∇ ρ - w ρ)
    ∂_t ρθ =  ∇ (μ ∇ ρθ - w ρθ)
    ∂_t u =  ∇ (μ ∇ u - w u) + (v - v_g)
    ∂_t v =  ∇ (μ ∇ v - w v) - (u - u_g)
    ∂_t w =  ∇ (μ ∇ w - w w) - g - c_p θ ∂_z Π

where
    Π = (p/p_0)^{R/c_p}

top BCs are insulating and impenetrable:
    ∂_t T = 0
    u = u_g
    v = v_g
    w = 0.0
    ∂_t ρ = 0

and bottom BCs use bulk formulae for surface fluxes of heat and momentum:
    ∂_t ρθ = F₃ = -Ch ρ ||u|| (T_sfc - ρθ / ρ)
    ∂_t u  = F₁ = -Cd u ||u||
    ∂_t v  = F₂ = -Cd v ||u||
    w = 0.0
    ∂_t ρ = 0

We also use this model to accumulate fluxes it calculates
    ∂_t F_accum = -(F₁, F₂, F₃)
=#

function ∑tendencies_atm!(dY, Y, (parameters, T_sfc), t)
    @unpack Cd, f, ν, uvg, C_p, MSLP, R_d, R_m, C_v, grav = parameters

    # unpack tendencies and state
    (Yc, Yf, F_sfc) = Y.x
    (dYc, dYf, dF_sfc) = dY.x

    UnPack.@unpack ρ, uv, ρθ = Yc

    w = Yf
    dρ = dYc.ρ
    duv = dYc.uv
    dρθ = dYc.ρθ
    dw = dYf

    # auxiliary calculations
    ρ_1 = parent(ρ)[1]
    ρθ_1 = parent(ρθ)[1]
    uv_1 = Operators.getidx(uv, Operators.Interior(), 1)
    u_wind = LinearAlgebra.norm(uv_1)

    # surface flux calculations
    fluxes = calculate_sfc_fluxes(DryMonin(), parameters, T_sfc[1], ρθ_1 / ρ_1, uv_1, ρ_1, t)

    surface_flux_ρθ = fluxes.SH
    surface_flux_uv = fluxes.τ

    # accumulate in the required right units
    @inbounds begin
        dY.x[3][1] = (surface_flux_ρθ)
        dY.x[3][2] = (surface_flux_uv)[1]
        dY.x[3][2] = (surface_flux_uv)[2]
    end

    # boundary conditions
    bcs_bottom_uv_flux = Operators.SetValue(surface_flux_uv)
    bcs_top_uv_value = Operators.SetValue(uvg)

    bcs_bottom_ρθ_flux = Operators.SetValue(Geometry.Cartesian3Vector(surface_flux_ρθ))
    bcs_top_ρθ_flux = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT)))

    bcs_bottom_ρ_flux = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT)))
    bcs_top_ρ_flux = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT)))

    # density
    If = Operators.InterpolateC2F()
    ∂f = Operators.GradientC2F()
    ∂c = Operators.DivergenceF2C(bottom = bcs_bottom_ρ_flux, top = bcs_top_ρ_flux)
    @. dρ = -∂c(w * If(ρ))

    # potential temperature
    If = Operators.InterpolateC2F()
    ∂f = Operators.GradientC2F()
    ∂c_adv = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
        top = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
    )
    ∂c_diff = Operators.DivergenceF2C(bottom = bcs_bottom_ρθ_flux, top = bcs_top_ρθ_flux)
    # TODO!: Undesirable casting to vector required

    @. dρθ = -∂c_adv(w * If(ρθ)) + ρ * ∂c_diff(Geometry.CartesianVector(ν * ∂f(ρθ / ρ)))

    A = Operators.AdvectionC2C(
        bottom = Operators.SetValue(Geometry.Cartesian12Vector(0.0, 0.0)),
        top = Operators.SetValue(Geometry.Cartesian12Vector(0.0, 0.0)),
    )

    # uv
    ∂c = Operators.DivergenceF2C(bottom = bcs_bottom_uv_flux)
    ∂f = Operators.GradientC2F(top = bcs_top_uv_value)
    duv .= (uv .- Ref(uvg)) .× Ref(Geometry.Cartesian3Vector(f))

    @. duv += ∂c(ν * ∂f(uv)) - A(w, uv)

    # w
    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    ∂f = Operators.GradientC2F()
    ∂c = Operators.GradientF2C()
    Af = Operators.AdvectionF2F()
    divf = Operators.DivergenceC2F()
    B = Operators.SetBoundaryOperator(
        bottom = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
        top = Operators.SetValue(Geometry.Cartesian3Vector(zero(FT))),
    )
    Φ(z) = grav * z
    Π(ρθ) = C_p * (R_d * ρθ / MSLP)^(R_m / C_v)
    zc = Fields.coordinate_field(axes(ρ))
    @. dw = B(Geometry.CartesianVector(-(If(ρθ / ρ) * ∂f(Π(ρθ))) - ∂f(Φ(zc))) + divf(ν * ∂c(w)) - Af(w, w))

    return dY
end

"""
Initialize fields located at cell centers in the vertical.
"""
@inline function init_ekman_column_1d_c(z, params::NamedTuple)
    @unpack grav, C_p, MSLP, R_d, T_surf_atm, T_min_ref, u0, v0, w0 = params

    T_surf = T_surf_atm
    Γ = grav / C_p
    T = max(T_surf - Γ * z, T_min_ref)
    p = MSLP * (T / T_surf)^(grav / (R_d * Γ))
    if T == T_min_ref
        z_top = (T_surf - T_min_ref) / Γ
        H_min = R_d * T_min_ref / grav
        p *= exp(-(z - z_top) / H_min)
    end
    θ = T_surf # potential temperature

    ρ = p / (R_d * θ * (p / MSLP)^(R_d / C_p))

    # velocity
    uv = Geometry.Cartesian12Vector(u0, v0) # u, v components

    # potential temperature
    ρθ = ρ * T_surf

    return (ρ = ρ, uv = uv, ρθ = ρθ)
end

"""
Initialize fields located at cell interfaces in the vertical.
"""
@inline function init_ekman_column_1d_f(z, params::NamedTuple)
    @unpack grav, C_p, MSLP, R_d, T_surf_atm, T_min_ref, u0, v0, w0 = params

    w = Geometry.Cartesian3Vector(w0) # w component

    return (w = w)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Normalize
using Test

@testset "Normalize.jl" begin
    include("test_boxcox.jl")
    include("test_yeojohnson.jl")
end


text:
using FactCheck
include("../src/code.jl")

# this non-test helps during the manual testing of Watcher.jl
# as the normal tests just finish too fast

println("starting runtests.jl")
println(Libc.strftime(time()))
sleep(3)
println("stopping unittest")


text:
export RawAcquisitionData_3DUTE

function RawAcquisitionData_3DUTE(b::BrukerFile)
    T = Complex{acqDataType(b)}

    filename = joinpath(b.path, "fid")
    filenameTraj = joinpath(b.path, "traj")

    N = acqSize(b)
    # The data is padded in case it is not a multiple of 1024
    # For multi-channel acquisition data at concatenate then padded to a multiple of 1024 bytes
    numChannel = pvmEncNReceivers(b)
    profileLength = Int((ceil(N[1]*numChannel*sizeof(T)/1024))*1024/sizeof(T))
    numAvailableChannel = pvmEncAvailReceivers(b)
    phaseFactor = acqPhaseFactor(b)
    numSlices = acqNumSlices(b)
    numEchos = acqNumEchos(b)
    numEncSteps2 = length(N) == 3 ? N[3] : 1
    numRep = acqNumRepetitions(b)

    I = open(filename,"r") do fd
      read!(fd,Array{T,7}(undef, profileLength,
                                     numEchos,
                                     phaseFactor,
                                     numSlices,
                                     div(N[2], phaseFactor),
                                     numEncSteps2,
                                     numRep))[1:N[1]*numChannel,:,:,:,:,:,:]
    end

    traj = open(filenameTraj,"r") do fd
        read!(fd,Array{Float64,3}(undef, 3,N[1],
                                       N[2]))
    end

    objOrd = acqObjOrder(b)
    objOrd = objOrd.-minimum(objOrd)

    gradMatrix = acqGradMatrix(b)

    offset1 = acqReadOffset(b)
    offset2 = acqPhase1Offset(b)
    offset3 = ndims(b) == 2 ? acqSliceOffset(b) : acqPhase2Offset(b)

    profiles = Profile[]
    for nR = 1:numRep
      for nEnc2 = 1:numEncSteps2
        for nPhase2 = 1:div(N[2], phaseFactor)
          for nSl = 1:numSlices
            for nPhase1 = 1:phaseFactor
              for nEcho=1:numEchos
                  counter = EncodingCounters(kspace_encode_step_1=0,
                                             kspace_encode_step_2=0,
                                             average=0,
                                             slice=objOrd[nSl],
                                             contrast=nEcho-1,
                                             phase=0,
                                             repetition=nR-1,
                                             set=0,
                                             segment=0 )

                  G = gradMatrix[:,:,nSl]
                  read_dir = (G[1,1],G[2,1],G[3,1])
                  phase_dir = (G[1,2],G[2,2],G[3,2])
                  slice_dir = (G[1,3],G[2,3],G[3,3])

                  # Not sure if the following is correct...
                  pos = offset1[nSl]*G[:,1] +
                        offset2[nSl]*G[:,2] +
                        offset3[nSl]*G[:,3]

                  position = (pos[1], pos[2], pos[3])

                  head = AcquisitionHeader(number_of_samples=N[1], idx=counter,
                                          read_dir=read_dir, phase_dir=phase_dir,
                                          slice_dir=slice_dir, position=position,
                                          center_sample=0,
                                          available_channels = numAvailableChannel, #TODO
                                          active_channels = numChannel,
                                          trajectory_dimensions = 3,
                                          sample_time_us = parse(Float64,b["PVM_DigDw"])/1000)

                  dat = map(T, reshape(I[:,nEcho,nPhase1,nSl,nPhase2,nEnc2,nR],:,numChannel))
                  push!(profiles, Profile(head,convert.(Float32,traj[:,:,nPhase2]),dat) )
              end
            end
          end
        end
      end
    end

    params = brukerParams(b)
    params["trajectory"] = "custom"
    params["encodedSize"] = pvmMatrix(b)
    return RawAcquisitionData(params, profiles)
end


text:
#######################
# Authorization Types #
#######################

abstract type Authorization end

# TODO: SecureString on 0.7
struct OAuth2 <: Authorization
    token::String
end

struct UsernamePassAuth <: Authorization
    username::String
    password::String
end

struct AnonymousAuth <: Authorization end

struct JWTAuth <: Authorization
    JWT::String
end

####################
# JWT Construction #
####################

function base64_to_base64url(string)
    replace(replace(replace(string, "=" => ""), '+' => '-'), '/' => '_')
end

function JWTAuth(app_id::Int, key::MbedTLS.PKContext; iat = now(Dates.UTC), exp_mins = 1)
    algo = base64_to_base64url(base64encode(JSON.json(Dict(
        "alg" => "RS256",
        "typ" => "JWT"
    ))))
    data = base64_to_base64url(base64encode(JSON.json(Dict(
        "iat" => trunc(Int64, Dates.datetime2unix(iat)),
        "exp" => trunc(Int64, Dates.datetime2unix(iat+Dates.Minute(exp_mins))),
        "iss" => app_id
    ))))
    signature = base64_to_base64url(base64encode(MbedTLS.sign(key, MbedTLS.MD_SHA256,
        MbedTLS.digest(MbedTLS.MD_SHA256, string(algo,'.',data)), RNG[])))
    JWTAuth(string(algo,'.',data,'.',signature))
end

function JWTAuth(app_id::Int, privkey::String; kwargs...)
    JWTAuth(app_id, MbedTLS.parse_keyfile(privkey); kwargs...)
end

###############
# API Methods #
###############

@api_default function authenticate(api::GitHubAPI, token::AbstractString; params = Dict(), options...)
    auth = OAuth2(token)
    params["access_token"] = auth.token
    gh_get(api, "/"; params = params, options...)
    return auth
end

#########################
# Header Authentication #
#########################

authenticate_headers!(headers, auth::AnonymousAuth) = headers

function authenticate_headers!(headers, auth::OAuth2)
    headers["Authorization"] = "token $(auth.token)"
    return headers
end

function authenticate_headers!(headers, auth::JWTAuth)
    headers["Authorization"] = "Bearer $(auth.JWT)"
    return headers
end

function authenticate_headers!(headers, auth::UsernamePassAuth)
    headers["Authorization"] = "Basic $(base64encode(string(auth.username, ':', auth.password)))"
    return headers
end

###################
# Pretty Printing #
###################

function Base.show(io::IO, a::OAuth2)
    token_str = a.token[1:6] * repeat("*", length(a.token) - 6)
    print(io, "GitHub.OAuth2($token_str)")
end


text:
export MtlError

const MTLError = Ptr{NsError}

mutable struct MtlError <: Exception
	code::NsInteger
	domain::String
	userinfo::String

	ptr::MTLError
end

Base.unsafe_convert(::Type{MTLError}, err::MtlError) = err.ptr

function MtlError(err::MTLError)
	code = mtErrorCode(err)
	domain = mtErrorDomain(err) |> unsafe_string
	userinfo = mtErrorUserInfo(err) |> unsafe_string

	obj = MtlError(code, domain, userinfo, err)
	finalizer(unsafe_destroy!, obj)
	return obj
end

function unsafe_destroy!(err::MtlError)
	mtRelease(err.ptr)
end

## properties

Base.propertynames(::MtlError) = (
	# error properties
	:code, :domain, :userInfo,
	# localized error descriptions
	:localizedDescription, :localizedRecoveryOptions,
	:localizedRecoverySuggestion, :localizedFailureReason)

function Base.getproperty(o::MtlError, f::Symbol)
    if f === :code
        mtErrorCode(o)
    elseif f === :domain
        unsafe_string(mtErrorDomain(o))
    elseif f === :userinfo
        ptr = mtErrorUserInfo(o)
        if ptr == C_NULL
			Dict{String,Any}()
		else
			unsafe_string(ptr)
		end
    elseif f === :localizedDescription
        unsafe_string(mtErrorLocalizedDescription(o))
    elseif f === :localizedRecoveryOptions
		count = Ref{Csize_t}(0)
		mtErrorLocalizedRecoveryOptions(o, count, C_NULL)
		options = Vector{String}(undef, count[])
		mtErrorLocalizedRecoveryOptions(o, count, options)
		unsafe_string.(options)
    elseif f === :localizedRecoverySuggestion
        ptr = mtErrorLocalizedRecoverySuggestion(o)
        ptr == C_NULL ? nothing : unsafe_string(ptr)
    elseif f === :localizedFailureReason
        ptr = mtErrorLocalizedFailureReason(o)
        ptr == C_NULL ? nothing : unsafe_string(ptr)
    else
        getfield(o, f)
    end
end

## display

function Base.showerror(io::IO, err::MtlError)
	print(io, "MtlError: $(err.localizedDescription) (code $(err.code), $(err.domain))")

	if err.localizedFailureReason !== nothing
		print(io, "\nFailure reason: $(err.localizedFailureReason)")
	end

	recovery_options = err.localizedRecoveryOptions
	if !isempty(recovery_options)
		print(io, "\nRecovery Options:")
		for option in recovery_options
			print(io, "\n - $(option)")
		end
	end
end


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries

using Test
using LinearAlgebra, SparseArrays

@testset "Tests with mixtures of Taylor1 and TaylorN" begin
    @test TaylorSeries.NumberNotSeries == Union{Real,Complex}
    @test TaylorSeries.NumberNotSeriesN == Union{Real,Complex,Taylor1}

    set_variables("x", numvars=2, order=6)
    xH = HomogeneousPolynomial(Int, 1)
    yH = HomogeneousPolynomial(Int, 2)
    tN = Taylor1(TaylorN{Float64}, 3)

    @test eltype(xH) == Int
    @test eltype(tN) == TaylorN{Float64}
    @test tN.order == 3
    @test string(zero(tN)) == "  0.0 + 𝒪(‖x‖¹) + 𝒪(t⁴)"
    @test string(tN) == " ( 1.0 + 𝒪(‖x‖¹)) t + 𝒪(t⁴)"
    @test string(tN + 3Taylor1(Int, 2)) == " ( 4.0 + 𝒪(‖x‖¹)) t + 𝒪(t³)"
    @test string(xH * tN) == " ( 1.0 x₁ + 𝒪(‖x‖²)) t + 𝒪(t⁴)"

    @test constant_term(xH) == xH
    @test constant_term(tN) == zero(TaylorN([xH]))
    @test linear_polynomial(xH) == xH
    @test linear_polynomial(1+tN+tN^2) == tN
    @test nonlinear_polynomial(1+tN+tN^2) == tN^2

    tN = Taylor1([zero(TaylorN(Float64,1)), one(TaylorN(Float64,1))], 3)
    @test typeof(tN) == Taylor1{TaylorN{Float64}}
    @test string(zero(tN)) == "  0.0 + 𝒪(‖x‖⁷) + 𝒪(t⁴)"
    @test string(tN) == " ( 1.0 + 𝒪(‖x‖⁷)) t + 𝒪(t⁴)"
    @test string(Taylor1([xH+yH])) == "  1 x₁ + 1 x₂ + 𝒪(t¹)"
    @test string(Taylor1([zero(xH), xH*yH])) == " ( 1 x₁ x₂) t + 𝒪(t²)"
    @test string(tN * Taylor1([0,TaylorN([xH+yH])])) == "  0.0 + 𝒪(‖x‖⁷) + 𝒪(t²)"

    t = Taylor1(3)
    xHt = HomogeneousPolynomial(typeof(t), 1)
    @test eltype(xHt) == Taylor1{Float64}
    @test string(xHt) == " ( 1.0 + 𝒪(t¹)) x₁"
    xHt = HomogeneousPolynomial([one(t), zero(t)])
    yHt = HomogeneousPolynomial([zero(t), t])
    @test string(xHt) == " ( 1.0 + 𝒪(t⁴)) x₁"
    @test string(yHt) == " ( 1.0 t + 𝒪(t⁴)) x₂"
    @test string(HomogeneousPolynomial([t])) == " ( 1.0 t + 𝒪(t⁴))"
    @test 3*xHt == HomogeneousPolynomial([3*one(t), zero(t)])
    @test t*xHt == HomogeneousPolynomial([t, zero(t)])
    @test complex(0,1)*xHt == HomogeneousPolynomial([1im*one(t), zero(1im*t)])
    @test eltype(complex(0,1)*xHt) == Taylor1{Complex{Float64}}
    @test (xHt+yHt)(1, 1) == 1+t
    @test (xHt+yHt)([1, 1]) == (xHt+yHt)((1, 1))

    tN1 = TaylorN([HomogeneousPolynomial([t]), xHt, yHt^2])
    @test tN1[0] == HomogeneousPolynomial([t])
    @test tN1(t,one(t)) == 2t+t^2
    @test tN1([t,one(t)]) == tN1((t,one(t)))
    t1N = convert(Taylor1{TaylorN{Float64}}, tN1)
    @test t1N[0] == HomogeneousPolynomial(1)
    ctN1 = convert(TaylorN{Taylor1{Float64}}, t1N)
    @test eltype(xHt) == Taylor1{Float64}
    @test eltype(tN1) == Taylor1{Float64}
    @test eltype(Taylor1([xH])) == HomogeneousPolynomial{Int}
    @test eltype(tN1) == Taylor1{Float64}
    @test get_order(HomogeneousPolynomial([Taylor1(1), 1.0+Taylor1(2)])) == 1
    @test 3*tN1 == TaylorN([HomogeneousPolynomial([3t]),3xHt,3yHt^2])
    @test t*tN1 == TaylorN([HomogeneousPolynomial([t^2]),xHt*t,t*yHt^2])
    @test string(tN1) ==
        " ( 1.0 t + 𝒪(t⁴)) + ( 1.0 + 𝒪(t⁴)) x₁ + ( 1.0 t² + 𝒪(t⁴)) x₂² + 𝒪(‖x‖³)"
    @test string(t1N) ==
        "  1.0 x₁ + 𝒪(‖x‖³) + ( 1.0 + 𝒪(‖x‖³)) t + ( 1.0 x₂² + 𝒪(‖x‖³)) t² + 𝒪(t⁴)"
    @test tN1 == ctN1
    @test tN1+tN1 == 2*tN1
    @test tN1+1im*tN1 == complex(1,1)*tN1
    @test tN1+t == t+tN1
    @test tN1-t == -t+tN1
    zeroN1 = zero(tN1)
    oneN1 = one(tN1)
    @test tN1-tN1 == zeroN1
    @test zero(zeroN1) == zeroN1
    @test zeroN1 == zero.(tN1)
    @test oneN1[0] == one(tN1[0])
    @test oneN1[1:end] == zero.(tN1[1:end])
    for i in eachindex(tN1.coeffs)
        @test tN1.coeffs[i].order == zeroN1.coeffs[i].order == oneN1.coeffs[i].order
    end
    @test string(t1N*t1N) ==
        "  1.0 x₁² + 𝒪(‖x‖³) + ( 2.0 x₁ + 𝒪(‖x‖³)) t + ( 1.0 + 𝒪(‖x‖³)) t² + ( 2.0 x₂² + 𝒪(‖x‖³)) t³ + 𝒪(t⁴)"
    @test !isnan(tN1)
    @test !isinf(tN1)

    @test mod(tN1+1,1.0) == 0+tN1
    @test mod(tN1-1.125,2) == 0.875+tN1
    @test (rem(tN1+1.125,1.0))[0][1] == 0.125 + t
    @test (rem(tN1-1.125,2))[0][1] == -1.125 + t
    @test mod2pi(-3pi+tN1)[0][1][0] ≈ pi
    @test mod2pi(0.125+2pi+tN1)[0][1][0] ≈ 0.125
    @test mod(t1N+1.125,1.0) == 0.125+t1N
    @test mod(t1N-1.125,2) == 0.875+t1N
    @test (rem(t1N+1.125,1.0))[0] == 0.125 + t1N[0]
    @test (rem(t1N-1.125,2))[0] == -1.125 + t1N[0]
    @test mod2pi(-3pi+t1N)[0][0][1] ≈ pi
    @test mod2pi(0.125+2pi+t1N)[0][0][1] ≈ 0.125

    @test abs(tN1+1) == 1+tN1
    @test abs(tN1-1) == 1-tN1
    @test_throws DomainError abs(tN1)
    @test_throws DomainError abs(t1N)

    @test convert(Array{Taylor1{TaylorN{Float64}},1}, [tN1, tN1]) == [t1N, t1N]
    @test convert(Array{Taylor1{TaylorN{Float64}},2}, [tN1 tN1]) == [t1N t1N]
    @test convert(Array{TaylorN{Taylor1{Float64}},1}, [t1N, t1N]) == [tN1, tN1]
    @test convert(Array{TaylorN{Taylor1{Float64}},2}, [t1N t1N]) == [tN1 tN1]

    @test evaluate(t1N, 0.0) == TaylorN(xH, 2)
    @test t1N() == TaylorN(xH, 2)
    @test string(evaluate(t1N, 0.0)) == " 1.0 x₁ + 𝒪(‖x‖³)"
    @test string(evaluate(t1N^2, 1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)"
    @test string((t1N^2)(1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)"
    v = zeros(TaylorN{Float64},2)
    @test evaluate!([t1N, t1N^2], 0.0, v) == nothing
    @test v == [TaylorN(1), TaylorN(1)^2]

    # Tests for functions of mixtures
    t1N = Taylor1([zero(TaylorN(Float64,1)), one(TaylorN(Float64,1))], 6)
    t = Taylor1(3)
    xHt = HomogeneousPolynomial([one(t), zero(t)])
    yHt = HomogeneousPolynomial([zero(t), t])
    x = TaylorN(1, order=2)
    y = TaylorN(2, order=2)
    xN1 = TaylorN([HomogeneousPolynomial(zero(t), 0), xHt, zero(yHt)], 2)
    yN1 = TaylorN([HomogeneousPolynomial(zero(t), 0), zero(xHt), yHt], 2)
    for fn in (exp, log, sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh, atanh)
        if fn == asin || fn == acos || fn == atanh
            cc = 0.5
        elseif fn == asinh || fn == acosh
            cc = 1.5
        else
            cc = 1.0
        end
        @test x*fn(cc+t1N) == fn(cc+t)*xN1
        @test t*fn(cc+xN1) == fn(cc+x)*t1N
    end

    vt = zeros(Taylor1{Float64},2)
    @test evaluate!([tN1, tN1^2], [t, t], vt) == nothing
    @test vt == [2t, 4t^2]

    tint = Taylor1(Int, 10)
    t = Taylor1(10)
    x = TaylorN( [HomogeneousPolynomial(zero(t), 5), HomogeneousPolynomial([one(t),zero(t)])], 5)
    y = TaylorN(typeof(tint), 2, order=5)
    @test typeof(x) == TaylorN{Taylor1{Float64}}
    @test eltype(y) == Taylor1{Int}
    @test -x == 0 - x
    @test +y == y
    @test one(y)/(1+x) == 1 - x + x^2 - x^3 + x^4 - x^5
    @test one(y)/(1+y) == 1 - y + y^2 - y^3 + y^4 - y^5
    @test (1+y)/one(t) == 1 + y
    @test typeof(y+t) == TaylorN{Taylor1{Float64}}

    # See #92 and #94
    δx, δy = set_variables("δx δy")
    xx = 1+Taylor1(δx,5)
    @test typeof(xx) == Taylor1{TaylorN{Float64}}
    @test eltype(xx) == TaylorN{Float64}
    @test !isnan(xx)
    @test !isnan(δx)
    @test !isinf(xx)
    @test !isinf(δx)
    @test +xx == xx
    @test -xx == 0 - xx
    @test xx/1.0 == 1.0*xx
    @test xx + xx == xx*2
    @test xx - xx == zero(xx)
    @test xx*xx == xx^2
    @test xx/xx == one(xx)
    @test xx*δx + Taylor1(typeof(δx),5) == δx + δx^2 + Taylor1(typeof(δx),5)
    @test xx/(1+δx) == one(xx)
    @test typeof(xx+δx) == Taylor1{TaylorN{Float64}}

    #testing evaluate and function-like behavior of Taylor1, TaylorN for mixtures:
    t = Taylor1(25)
    p = cos(t)
    q = sin(t)
    a = [p,q]
    dx = set_variables("x", numvars=4, order=10)
    P = sin.(dx)
    v = [1.0,2,3,4]
    F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))]
    Q = F(v+dx)
    diff_evals = cos(sin(dx[1]))-p(P[1])
    @test norm( norm.(map(x->x.coeffs, diff_evals.coeffs),Inf) , Inf) < 1e-15
    #evaluate a Taylor1 at a TaylorN
    @test p(P) == evaluate(p, P)
    @test q(Q) == evaluate(q, Q)
    #evaluate an array of Taylor1s at a TaylorN
    aT1 = [p,q,p^2,log(1+q)] #an array of Taylor1s
    @test aT1(Q[4]) == evaluate(aT1, Q[4])
    @test (aT1.^2)(Q[3]) == evaluate(aT1.^2, Q[3])
    #evaluate a TaylorN at an array of Taylor1s
    @test P[1](aT1) == evaluate(P[1], aT1)
    @test P[1](aT1) == evaluate(P[1], (aT1...,))
    @test Q[2](aT1) == evaluate(Q[2], [aT1...])
    #evaluate an array of TaylorN{Float64} at an array of Taylor1{Float64}
    @test P(aT1) == evaluate(P, aT1)
    @test Q(aT1) == evaluate(Q, aT1)
    #test evaluation of an Array{TaylorN{Taylor1}} at an Array{Taylor1}
    aH1 = [
        HomogeneousPolynomial([Taylor1(rand(2))]),
        HomogeneousPolynomial([Taylor1(rand(2)),Taylor1(rand(2)),
            Taylor1(rand(2)),Taylor1(rand(2))])
        ]
    bH1 = [
        HomogeneousPolynomial([Taylor1(rand(2))]),
        HomogeneousPolynomial([Taylor1(rand(2)),Taylor1(rand(2)),
            Taylor1(rand(2)),Taylor1(rand(2))])
        ]
    aTN1 = TaylorN(aH1); bTN1 = TaylorN(bH1)
    x = [aTN1, bTN1]
    δx = [Taylor1(rand(3)) for i in 1:4]
    @test typeof(x) == Array{TaylorN{Taylor1{Float64}},1}
    @test typeof(δx) == Array{Taylor1{Float64},1}
    x0 = Array{Taylor1{Float64}}(undef, length(x))
    eval_x_δx = evaluate(x,δx)
    @test x(δx) == eval_x_δx
    evaluate!(x,δx,x0)
    @test x0 == eval_x_δx
    @test typeof(evaluate(x[1],δx)) == Taylor1{Float64}
    @test x() == map(y->y[0][1], x)
    for i in eachindex(x)
        @test evaluate(x[i],δx) == eval_x_δx[i]
        @test x[i](δx) == eval_x_δx[i]
    end
    p11 = Taylor1([sin(t),cos(t)])
    @test evaluate(p11,t) == sin(t)+t*cos(t)
    @test p11(t) == sin(t)+t*cos(t)
    a11 = Taylor1([t,t^2,exp(-t),sin(t),cos(t)])
    b11 = t+t*(t^2)+(t^2)*(exp(-t))+(t^3)*sin(t)+(t^4)*cos(t)
    diff_a11b11 = a11(t)-b11
    @test norm(diff_a11b11.coeffs,Inf) < 1E-19

    X, Y = set_variables(Taylor1{Float64}, "x y")
    @test typeof( norm(X) ) == Float64
    @test norm(X) > 0
    @test norm(X+Y) == sqrt(2)
    @test norm(-10X+4Y,Inf) == 10.

    X,Y = convert(Taylor1{TaylorN{Float64}},X), convert(Taylor1{TaylorN{Float64}},Y)
    @test typeof( norm(X) ) == Float64
    @test norm(X) > 0
    @test norm(X+Y) == sqrt(2)
    @test norm(-10X+4Y,Inf) == 10.

    @test TaylorSeries.rtoldefault(TaylorN{Taylor1{Int}}) == 0
    @test TaylorSeries.rtoldefault(Taylor1{TaylorN{Int}}) == 0
    for T in (Float64, BigFloat)
        @test TaylorSeries.rtoldefault(TaylorN{Taylor1{T}}) == sqrt(eps(T))
        @test TaylorSeries.rtoldefault(Taylor1{TaylorN{T}}) == sqrt(eps(T))
        @test TaylorSeries.real(TaylorN{Taylor1{T}}) == TaylorN{Taylor1{T}}
        @test TaylorSeries.real(Taylor1{TaylorN{T}}) == Taylor1{TaylorN{T}}
        @test TaylorSeries.real(TaylorN{Taylor1{Complex{T}}}) == TaylorN{Taylor1{T}}
        @test TaylorSeries.real(Taylor1{TaylorN{Complex{T}}}) == Taylor1{TaylorN{T}}
    end

    rndT1(ord1) = Taylor1(-1 .+ 2rand(ord1+1)) # generates a random Taylor1 with order `ord`
    nmonod(s, d) = binomial(d+s-1, d) #number of monomials in s variables with exact degree d
    #rndHP generates a random `ordHP`-th order homog. pol. of Taylor1s, each with order `ord1`
    rndHP(ordHP, ord1) = HomogeneousPolynomial( [rndT1(ord1) for i in 1:nmonod(get_numvars(), ordHP)] )
    #rndTN generates a random `ordHP`-th order TaylorN of of Taylor1s, each with order `ord1`
    rndTN(ordN, ord1) = TaylorN([rndHP(i, ord1) for i in 0:ordN])

    P = rndTN(get_order(), 3)
    @test P ≈ P
    Q = deepcopy(P)
    Q[2][2] = Taylor1([NaN, Inf])
    @test isnan(Q)
    @test isinf(Q)
    @test !isfinite(Q)
    Q[2][2] = P[2][2]+sqrt(eps())/2
    @test isapprox(P, Q, rtol=1.0)
    Q[2][2] = P[2][2]+10sqrt(eps())
    @test !isapprox(P, Q, atol=sqrt(eps()), rtol=0)
    @test P ≉ Q^2
    Q[2][2] = P[2][2]+eps()/2
    @test isapprox(Q, Q, atol=eps(), rtol=0)
    @test isapprox(Q, P, atol=eps(), rtol=0)
    Q[2][1] = P[2][1]-10eps()
    @test !isapprox(Q, P, atol=eps(), rtol=0)
    @test P ≉ Q^2

    X, Y = set_variables(BigFloat, "x y", numvars=2, order=6)
    p1N = Taylor1([X^2,X*Y,Y+X,Y^2])
    q1N = Taylor1([X^2,(1.0+sqrt(eps(BigFloat)))*X*Y,Y+X,Y^2])
    @test p1N ≈ p1N
    @test p1N ≈ q1N

    Pv = [rndTN(get_order(), 3), rndTN(get_order(), 3)]
    Qv = convert.(Taylor1{TaylorN{Float64}}, Pv)

    @test TaylorSeries.jacobian(Pv) == TaylorSeries.jacobian(Qv)

    @test_throws ArgumentError Taylor1(2) + TaylorN(1)
    @test_throws ArgumentError Taylor1(2) - TaylorN(1)
    @test_throws ArgumentError Taylor1(2) * TaylorN(1)
    @test_throws ArgumentError TaylorN(2) / Taylor1(1)
end

@testset "Tests with nested Taylor1s" begin
    ti = Taylor1(3)
    to = Taylor1([zero(ti), one(ti)], 9)
    @test string(to) == " ( 1.0 + 𝒪(t⁴)) t + 𝒪(t¹⁰)"
    @test string(to^2) == " ( 1.0 + 𝒪(t⁴)) t² + 𝒪(t¹⁰)"
    @test ti + to == Taylor1([ti, one(ti)], 9)
    tito = ti * to
    @test tito == Taylor1([zero(ti), ti], 9)
    @test tito / to == ti
    @test get_order(tito/to) == get_order(to)-1
    @test tito / ti == to
    @test get_order(tito/ti) == get_order(to)
    # @test get_order((tito/ti)[1]) == get_order(ti)-1
    @test ti^2-to^2 == (ti+to)*(ti-to)
    @test sin(to) ≈ Taylor1(one(ti) .* sin(Taylor1(10)).coeffs, 9)
    @test to(1 + ti) == 1 + ti
    @test to(1 + ti) isa Taylor1{Float64}
    @test ti(1 + to) == 1 + to
    @test constant_term(ti+to) == ti
    @test linear_polynomial(ti*to) == Taylor1([zero(ti), ti], 9)
    @test get_order(linear_polynomial(to)) == get_order(to)
    @test nonlinear_polynomial(to+ti*to^2) == Taylor1([zero(ti), zero(ti), ti], 9)
    @test ti(1 + to) isa Taylor1{Taylor1{Float64}}
    @test sqrt(tito^2) == tito
    @test get_order(sqrt(tito^2)) == get_order(to) >> 1
    @test (tito^3)^(1/3) == tito
    @test get_order(sqrt(tito^2)) == get_order(to) >> 1
    ti2to = ti^2 * to
    tti = (ti2to/to)/ti
    @test get_order(tti) == get_order(to)-1
    @test get_order(tti[0]) == get_order(ti)-1
end


text:
"""
    SimpleWeightedDiGraph{T, U}

A type representing a directed graph with weights of type `U`.

Note that adding or removing vertices or edges is not particularly performant;
see MetaGraphs.jl for possible alternatives.

The primary constructor takes a sparse adjacency matrix as input, of which
the transpose is stored.
To provide the transpose directly, the keyword argument `permute=true` can be used.
"""
mutable struct SimpleWeightedDiGraph{T<:Integer, U<:Real} <: AbstractSimpleWeightedGraph{T, U}
    weights::SparseMatrixCSC{U, T} # indexed by [dst, src]
    function SimpleWeightedDiGraph{T, U}(adjmx::SparseMatrixCSC{U,T}; permute=true) where {T <: Integer, U <: Real}
        dima, dimb = size(adjmx)
        isequal(dima, dimb) || error("Adjacency / distance matrices must be square")
        permute ? new{T, U}(permutedims(adjmx)) : new{T, U}(adjmx)
    end

end

SimpleWeightedDiGraph{T}(adjmx::SparseMatrixCSC{U, T}; permute=true) where {T<:Integer, U<:Real} =
    SimpleWeightedDiGraph{T, U}(adjmx; permute=permute)

SimpleWeightedDiGraph(adjmx::SparseMatrixCSC{U, T}; permute=true) where {T<:Integer, U<:Real} =
    SimpleWeightedDiGraph{T, U}(adjmx; permute=permute)

SimpleWeightedDiGraph(m::AbstractMatrix{U}) where {U <: Real} =
    SimpleWeightedDiGraph{Int, U}(SparseMatrixCSC{U, Int}(m))
SimpleWeightedDiGraph{T}(m::AbstractMatrix{U}) where {T <: Integer, U <: Real} =
    SimpleWeightedDiGraph{T, U}(SparseMatrixCSC{U, T}(m))
SimpleWeightedDiGraph{T, U}(m::AbstractMatrix) where {T <: Integer, U <: Real} =
    SimpleWeightedDiGraph{T, U}(SparseMatrixCSC{U, T}(m))

SimpleWeightedDiGraph(g::SimpleWeightedDiGraph) = SimpleWeightedDiGraph(copy(g.weights), permute=false)
SimpleWeightedDiGraph{T,U}(g::SimpleWeightedDiGraph) where {T <: Integer, U <: Real} =
    SimpleWeightedDiGraph(SparseMatrixCSC{U, T}(copy(g.weights)), permute=false)

ne(g::SimpleWeightedDiGraph) = nnz(g.weights)

function SimpleWeightedDiGraph{T,U}(n::Integer = 0) where {T<:Integer, U<:Real}
    weights = spzeros(U, T, T(n), T(n))
    return SimpleWeightedDiGraph{T, U}(weights)
end

# Graph()
SimpleWeightedDiGraph() = SimpleWeightedDiGraph{Int, Float64}()

# Graph(6), Graph(0x5)
SimpleWeightedDiGraph(n::T) where T<:Integer = SimpleWeightedDiGraph{T, Float64}(n)

# Graph(UInt8)
SimpleWeightedDiGraph(::Type{T}) where T<:Integer = SimpleWeightedDiGraph{T, Float64}(zero(T))

# Graph(UInt8, Float32)
SimpleWeightedDiGraph(::Type{T}, ::Type{U}) where {T <: Integer, U <: Real} = SimpleWeightedDiGraph{T, U}(zero(T))

# DiGraph(AbstractGraph, ::Type{U})
function SimpleWeightedDiGraph(g::Graphs.AbstractGraph{T}, ::Type{U}=Float64) where {U <: Real, T}
    return SimpleWeightedDiGraph{T}(adjacency_matrix(g, U))
end

"""
    DiGraph(g::AbstractGraph, x::Real)

Construct a weighted digraph from other graph `g` with initial weight `x`.
"""
function SimpleWeightedDiGraph(g::Graphs.AbstractGraph{T}, x::U) where {U <: Real, T}
    m = adjacency_matrix(g, U)'
    return SimpleWeightedDiGraph{T, U}(x .* m, permute=false)
end

# DiGraph(srcs, dsts, weights)
function SimpleWeightedDiGraph(i::AbstractVector{T}, j::AbstractVector{T}, v::AbstractVector{U}; combine = +) where T<:Integer where U<:Real
    m = max(maximum(j), maximum(i))
    SimpleWeightedDiGraph{T, U}(sparse(j, i, v, m, m, combine), permute=false)
end

Graphs.SimpleDiGraph(g::SimpleWeightedDiGraph) = SimpleDiGraph(g.weights')

edgetype(::SimpleWeightedDiGraph{T, U}) where T<:Integer where U<:Real = SimpleWeightedGraphEdge{T,U}

edges(g::SimpleWeightedDiGraph) = (SimpleWeightedEdge(x[2], x[1], x[3]) for x in zip(findnz(g.weights)...))
weights(g::SimpleWeightedDiGraph) = g.weights'

inneighbors(g::SimpleWeightedDiGraph, v::Integer) = g.weights[v,:].nzind

# add_edge! will overwrite weights.
function add_edge!(g::SimpleWeightedDiGraph, e::SimpleWeightedGraphEdge)
    T = eltype(g)
    U = weighttype(g)
    s_, d_, w = Tuple(e)

    if w == zero(U)
        @warn "Note: adding edges with a zero weight to this graph type has no effect." maxlog=1 _id=:swd_add_edge_zero
        return false
    end

    s = T(s_)
    d = T(d_)
    (s in vertices(g) && d in vertices(g)) || return false
    @inbounds g.weights[d, s] = w
    return true
end

function rem_edge!(g::SimpleWeightedDiGraph, e::SimpleWeightedGraphEdge)
    has_edge(g, e) || return false
    U = weighttype(g)
    @inbounds g.weights[dst(e), src(e)] = zero(U)
    return true
end

copy(g::SimpleWeightedDiGraph) =  SimpleWeightedDiGraph(copy(g.weights'))

==(g::SimpleWeightedDiGraph, h::SimpleWeightedDiGraph) = g.weights == h.weights

is_directed(::Type{<:SimpleWeightedDiGraph}) = true

"""
    g[e::SimpleWeightedGraph, Val{:weight}()]

Equivalent to g[src(e), dst(e)].
"""
function Base.getindex(g::SimpleWeightedDiGraph{T, U}, e::AbstractEdge, ::Val{:weight}) where {T, U, S}
    return g.weights[dst(e), src(e)]
end

"""
    g[e::SimpleWeightedGraph, i::Integer, j::Integer, Val{:weight}()]

Return the weight of edge (i, j).
"""
function Base.getindex(g::SimpleWeightedDiGraph{T, U}, i::Integer, j::Integer, ::Val{:weight}) where {T, U, S}
    return g.weights[j, i]
end


text:
module TestSquares
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.CategoricalAlgebra.Squares

t  = FinFunction([1,2,4], 5)
l  = FinFunction([1,2,3], 3)
b = FinFunction([2,3,1], 3)
f  = FinFunction([2,3,2,1,3], 3)

α = SquareDiagram(t, b, l, f)

t₁ = FinFunction([1,1,1,1,1], 1)
r  = FinFunction([1], 1)
b₁ = FinFunction([1,1,1], 1)

β = SquareDiagram(t₁, b₁, f, r)
γ = composeH(α, β)
@test top(γ) == t⋅t₁
@test bottom(γ) == b⋅b₁
@test left(γ) == left(α)
@test right(γ) == right(β)

t = FinFunction([1,2,4], 5)
b = FinFunction([2,1,3], 3)
l = FinFunction([2,1,3], 3)
r = FinFunction([1,2,2,3,3], 3)
α = SquareDiagram(t,b,l,r)
β = SquareDiagram(FinFunction([1,2,1,2,1], 2),
                  FinFunction([1,1,1],1),
                  FinFunction([1,2,2,3,3], 3),
                  FinFunction([1,1],1)
)
γ = composeH(α, β)
@test collect(top(γ)) == [1,2,2]
@test collect(bottom(γ)) == [1,1,1]

α = SquareDiagram(l,r, t,b)
β = SquareDiagram(FinFunction([1,2,2,3,3], 3),
                  FinFunction([1,1],1),
                  FinFunction([1,2,1,2,1], 2),
                  FinFunction([1,1,1],1)
)
γ = composeV(α, β)
@test_throws AssertionError composeH(α, β)
@test collect(left(γ)) == [1,2,2]
@test collect(right(γ)) == [1,1,1]

@test collect(idH(FinSet(3))) == collect(1:3)
@test collect(idV(FinSet(5))) == collect(1:5)
@test collect(idH(FinSet(0))) == []

@test collect(composeH(FinFunction([2,1], 3), FinFunction([2, 1, 3]))) == [1,2]
@test collect(composeV(FinFunction([2,1], 3), FinFunction([2, 1, 3]))) == [1,2]

end


text:
using Statistics
using DataFrames
using JLD2, FileIO

using Cairo, Gadfly
Gadfly.set_default_plot_size(40cm,20cm)

## plotting

sap_ = getFactorType(newfg, :x0x1f1)

##

Caesar.plotScatterAlign(sap_, bw=0.001)

##

# snt_ = deepcopy(snt)
snt = overlayScatterMutate(sap_; user_offset=[-20.0;+100; pi/4], sample_count=200, bw=1e-4, user_coords=[0;-50;0.]) #, user_coords=snt_.best_coords)
plotScatterAlign(snt)

## sensitivity sweep on bw

SNT_10 = [[overlayScatterMutate(sap_; user_offset=[-20.0;+20;pi/4], sample_count=100, bw) for i in 1:200] for bw in exp10.(-6:1:-1)];

mkdir("/tmp/caesar/rex/SNT_10")
@save "/tmp/caesar/rex/SNT_10/SNT_10.jld2" SNT_10

## histogram sensitivities

df = DataFrame(s=Int[], bw=Float64[], group=Symbol[], μ=Float64[], σ=Float64[], _μ=Float64[], μ_=Float64[] )

for (i,s) in enumerate(-6:-1)

  tl = "e$s"
  X_ = (x->x.best_coords[1]).(SNT[i])
  Y_ = (x->x.best_coords[2]).(SNT[i])
  R_ = (x->x.best_coords[3]).(SNT[i])

  μ_x = Statistics.mean(X_); σ_x = Statistics.std(X_)
  μ_y = Statistics.mean(Y_); σ_y = Statistics.std(Y_)
  μ_r = Statistics.mean(R_); σ_r = Statistics.std(R_)

  _μx, μx_ = μ_x - σ_x, μ_x + σ_x
  _μy, μy_ = μ_y - σ_y, μ_y + σ_y
  _μr, μr_ = μ_r - σ_r, μ_r + σ_r

  push!( df, ( s,exp10(s),:x, μ_x,σ_x, _μx,μx_ ) )
  push!( df, ( s,exp10(s),:y, μ_y,σ_y, _μy,μy_ ) )
  push!( df, ( s,exp10(s),:r, μ_r,σ_r, _μr,μr_ ) )

  # Hx = Gadfly.plot(x=X_, Geom.histogram, Guide.title("X, $tl"))
  # Hy = Gadfly.plot(x=Y_, Geom.histogram, Guide.title("Y, $tl"))
  # Ht = Gadfly.plot(x=R_, Geom.histogram, Guide.title("θ, $tl"))
  # hstack(Hx, Hy, Ht);
end

##

df_xy = filter( row -> row.group != :r, df )
df_r  = filter( row -> row.group == :r,  df )

p_xy = Gadfly.plot(df_xy, x=:bw, y=:μ, ymin=:_μ, ymax=:μ_, color=:group, Geom.line, Geom.errorbar, Scale.x_log10)
p_r  = Gadfly.plot(df_r,  x=:bw, y=:μ, ymin=:_μ, ymax=:μ_, color=:group, Geom.line, Geom.errorbar, Scale.x_log10)

hstack(p_xy, p_r)

##=============================================================================================
## dev testing on one scatter image with known transform

# sweeps_ = sweeps[5:6]
# sweeps_[2] = sweeps[5]

# domain = (range(-100,100,length=976),range(-100,100,length=976))

# newfg = buildGraphChain!( sweeps_,
#                           ScatterAlignPose2,
#                           (_, data) -> (data.currData,data.nextData,domain),
#                           stopAfter=1,
#                           doRef = false,
#                           inflation_fct=0.0,
#                           solverParams=SolverParams(graphinit=false, inflateCycles=1) )
# #

##

# fct = getFactorType(newfg, :x0x1f1)
# Xtup = sampleFactor(newfg, :x0x1f1)
# e0 = identity_element(SpecialEuclidean(2))
# δ = calcFactorResidualTemporary(fct, (Pose2,Pose2), Xtup, (e0, e0))

#


text:
# See https://en.wikipedia.org/wiki/Matrix_(mathematics)
# Examples from https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/

using LinearAlgebra

expressions = [
    :(inv([1 2 3; 4 1 6; 7 8 1]))
    :(det([1 2 3; 4 1 6; 7 8 1]))
    :(eigvals([-4. -17.; 2. 2.]))
    :(eigvecs([-4. -17.; 2. 2.]))
    :(factorize([1.5 2 -4; 3 -1 -6; -10 2.3 4]))
    :(factorize([1.5 2 -4; 2 -1 -3; -4 -3 5]))
    :(([1 0; 1 -2] \ [32; -4]))
    :(([1 0; 1 -2] * [32.0, 18.0]))
    :(dot([1; 1], [2; 3]))
    :(cross([0; 1; 0], [0; 0; 1]))
    :(UpperHessenberg([1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]))
    :(cholesky([4. 12. -16.; 12. 37. -43.; -16. -43. 98.]).U)
    :(cholesky([4. 12. -16.; 12. 37. -43.; -16. -43. 98.]).L)
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
module AMDGPU

### Imports ###

using CEnum
using Setfield
using Libdl
using LLVM, LLVM.Interop
using GPUArrays
using GPUCompiler
using Adapt
using Requires
import LinearAlgebra
using Printf
import Core: LLVMPtr

### Exports ###

export HSAAgent, HSAQueue, HSAExecutable, HSAKernelInstance, HSASignal
export get_agents, profile, get_first_isa_string, launch!
export get_default_agent, get_default_queue
export has_rocm_gpu

export ROCArray, ROCVector, ROCMatrix, ROCVecOrMat
export roc

### HSA Runtime ###

const libhsaruntime = "libhsa-runtime64.so.1"
include(joinpath(@__DIR__, "hsa", "HSA.jl"))
import .HSA: Agent, Queue, Executable, Status, Signal

# Load binary dependencies
include(joinpath(dirname(@__DIR__), "deps", "loaddeps.jl"))

struct Adaptor end

const RT_LOCK = Threads.ReentrantLock()

include("extras.jl")
include("error.jl")
include("agent.jl")
include("queue.jl")
include("signal.jl")
include("executable.jl")
include("kernel.jl")
include("memory.jl")

@enum DeviceRuntime HSA_rt OCL_rt
const RUNTIME = Ref{DeviceRuntime}(HSA_rt)
#=
if get(ENV, "AMDGPUNATIVE_OPENCL", "") != ""
    RUNTIME[] = OCL_rt
end
=#
include("runtime.jl")
include("statussignal.jl")
include("sync.jl")

const ci_cache = GPUCompiler.CodeCache()
Base.Experimental.@MethodTable(method_table)

# Device sources must load _before_ the compiler infrastructure
# because of generated functions.
include(joinpath("device", "tools.jl"))
include(joinpath("device", "pointer.jl"))
include(joinpath("device", "array.jl"))
include(joinpath("device", "gcn.jl"))
include(joinpath("device", "runtime.jl"))
include(joinpath("device", "llvm.jl"))
include(joinpath("device", "globals.jl"))

include("query.jl")
include("compiler.jl")
include("execution_utils.jl")
include("execution.jl")
include("exceptions.jl")
include("reflection.jl")

### ROCArray ###

include("array.jl")
#include("subarray.jl")
#include("indexing.jl")
include("broadcast.jl")
#include("matmul.jl")
include("mapreduce.jl")
#include("gpuarray_interface.jl")

allowscalar(x::Bool) = GPUArrays.allowscalar(x)

### Initialization and Shutdown ###

const HSA_REFCOUNT = Threads.Atomic{UInt}(0)
function hsaref!()
    #=
    if Threads.atomic_add!(HSA_REFCOUNT, UInt(1)) > typemax(UInt)-10
        Core.println("HSA_REFCOUNT OVERFLOW!")
        exit(1)
    end
    =#
end
function hsaunref!()
    #=
    if Threads.atomic_sub!(HSA_REFCOUNT, UInt(1)) == 1
        HSA.shut_down()
    end
    =#
end

# Utilities
include("utils.jl")

# Load HIP and ROCm external libraries
const libhip = "libamdhip64.so"
if functional(:hip)
include(joinpath(@__DIR__, "hip", "HIP.jl"))
functional(:rocblas)      && include(joinpath(@__DIR__, "blas", "rocBLAS.jl"))
#functional(:rocsparse)  && include("sparse/rocSPARSE.jl")
#functional(:rocsolver)   && include("solver/rocSOLVER.jl")
#functional(:rocalution) && include("solver/rocALUTION.jl")
if functional(:rocrand)
include(joinpath(@__DIR__, "rand", "rocRAND.jl"))
include(joinpath(@__DIR__, "random.jl"))
end
functional(:rocfft)      && include(joinpath(@__DIR__, "fft", "rocFFT.jl"))
#functional(:miopen)     && include("dnn/MIOpen.jl")

# Ensure external libraries are up to date
function check_library(name, path)
    path === nothing && return
    if !ispath(path)
        @warn "$name library has changed. Please run Pkg.build(\"AMDGPU\") and restart Julia."
    end
end
check_library("rocBLAS", librocblas)
check_library("rocSPARSE", librocsparse)
check_library("rocSOLVER", librocsolver)
check_library("rocALUTION", librocalution)
check_library("rocRAND", librocrand)
check_library("rocFFT", librocfft)
check_library("MIOpen", libmiopen)
end # functional(:hip)

function __init__()
    if !configured && build_reason != "unknown"
        if build_reason == "Build did not occur"
            @warn """
            AMDGPU.jl has not been built
            Please run Pkg.build("AMDGPU") and restart Julia
            """
            return
        end

        # Some other reason
        @warn """
        Build failed to start!
        Reason: $build_reason
        """
        return
    end
    if hsa_configured
        # Make sure we load the library found by the last `] build`
        push!(Libdl.DL_LOAD_PATH, dirname(libhsaruntime_path))
        # TODO: Do the same (if possible) for the debug library

        # Initialize the HSA runtime
        HSA_REFCOUNT[] = 0
        status = HSA.init()
        if status == HSA.STATUS_SUCCESS
            hsaref!()
            HSA_REFCOUNT[] = 1
            # Register shutdown hook
            atexit() do
                hsaunref!()
            end

            # Populate the default agent
            agents = get_agents(:gpu)
            if length(agents) > 0
                DEFAULT_AGENT[] = first(agents)
            end
        else
            @warn "HSA initialization failed with code $status"
        end
    else
        @warn """
        HSA runtime has not been built, runtime functionality will be unavailable.
        Please run Pkg.build("AMDGPU") and reload AMDGPU.
        Reason: $hsa_build_reason
        """

        if parse(Bool, get(ENV, "JULIA_AMDGPU_CORE_MUST_LOAD", "0"))
            print_build_diagnostics()
            error("Failed to load HSA runtime, but HSA must load, bailing out")
        end
    end

    # Check whether ld.lld was found
    if !functional(:lld)
        @warn """
        ld.lld was not found, compilation functionality will be unavailable.
        Please run Pkg.build("AMDGPU") and reload AMDGPU.
        Reason: $lld_build_reason
        """

        if parse(Bool, get(ENV, "JULIA_AMDGPU_CORE_MUST_LOAD", "0"))
            print_build_diagnostics()
            error("Failed to find ld.lld, but ld.lld must exist, bailing out")
        end
    end

    # Check whether device intrinsics are available
    if !functional(:device_libs)
        @warn """
        ROCm-Device-Libs were not found, device intrinsics will be unavailable.
        Please run Pkg.build("AMDGPU") and reload AMDGPU.
        Reason: $device_libs_build_reason
        """

        if parse(Bool, get(ENV, "JULIA_AMDGPU_CORE_MUST_LOAD", "0"))
            print_build_diagnostics()
            error("Failed to find Device Libs, but Device Libs must exist, bailing out")
        end
    end

    # Check whether HIP is available
    if functional(:hip)
        push!(Libdl.DL_LOAD_PATH, dirname(libhip_path))
    else
        @warn """
        HIP library has not been built, HIP integration will be unavailable.
        Please run Pkg.build("AMDGPU") and reload AMDGPU.
        Reason: $hip_build_reason
        """

        if parse(Bool, get(ENV, "JULIA_AMDGPU_HIP_MUST_LOAD", "0"))
            print_build_diagnostics()
            error("Failed to load HIP runtime, but HIP must load, bailing out")
        end
    end

    # Check whether external libraries are available
    if use_artifacts && !functional(:rocrand)
        @warn """
        rocRAND_jll failed to load, RNG functionality will be unavailable.
        Please run Pkg.build("AMDGPU") and reload AMDGPU.
        Reason: $rocrand_build_reason
        """
    end

    # Load optional OpenCL integrations
    @require OpenCL="08131aa3-fb12-5dee-8b74-c09406e224a2" include("opencl.jl")
end

end # module


text:
@doc raw"""
    Bbox3

An object `b` of the type [`Bbox3`](@ref) is a bounding box in the
three-dimensional Euclidean space ``Ε^3``.
"""
Bbox3

@doc raw"""
    Bbox3()

Introduces an *empty* bounding box with lower left corner point at ``(∞, ∞, ∞)``
and with upper right corner point at ``(-∞, -∞, -∞)``, ``∞`` being
[`Inf`](@ref).
"""
Bbox3()

"""
    Bbox3(xmin::Float64, ymin::Float64, zmin::Float64, xmax::Float64, ymax::Float64, zmax::Float64)
    Bbox3(xmin::Real, ymin::Real, zmin::Real, xmax::Real, ymax::Real, zmax::Real)

Introduces a bounding box `b` with lower left corner at `(xmin, ymin, zmin)` and
with upper right corner at `(xmax, ymax, zmax)`.
"""
Bbox3(xmin::Real, ymin::Real, zmin::Real, xmax::Real, ymax::Real, zmax::Real) =
    Bbox3(convert(Float64, xmin), convert(Float64, ymin), convert(Float64, zmin),
          convert(Float64, xmax), convert(Float64, ymax), convert(Float64, zmax))

"""
    dimension(b::Bbox3)

Returns 3.
"""
dimension(b::Bbox3)

"""
    min(b::Bbox3, i::Integer)

Returns [`xmin(b)`](@ref) if `i == 0` or [`ymin(b)`](@ref) if `i == 1` or
[`zmin(b)`](@ref) if `i == 2`.

!!! info "Precondition"

    `i >= 0` and `i <= 2`.
"""
min(b::Bbox3, i::Integer)

"""
    max(b::Bbox3, i::Integer)

Returns [`xmax(b)`](@ref) if `i == 0` or [`ymax(b)`](@ref) if `i == 1` or [`zmax(b)`](@ref) if `i == 2`.

!!! info "Precondition"

    `i >= 0` or `i <= 1`.
"""
max(b::Bbox3, i::Integer)

"""
    +(b::Bbox3, c::Bbox3)

Returns a bounding box of `b` and `c`.
"""
+(b::Bbox3, c::Bbox3)

"""
    dilate(b::Bbox3, dist::Integer)

Dilates the bounding box by a specified number of ULP.
"""
dilate(b::Bbox3, dist::Integer)


text:
using Test
using ConstraintSolver
using JSON
using Random
using MathOptInterface, JuMP, Cbc, GLPK, Combinatorics
using ReferenceTests
using LinearAlgebra

const MOI = MathOptInterface
const CS = ConstraintSolver
const MOIU = MOI.Utilities

function CSTestOptimizer(; branch_strategy = :Auto)
    CS.Optimizer(logging = [], seed = 1, branch_strategy = branch_strategy)
end
function CSJuMPTestOptimizer(; branch_strategy = :Auto)
    JuMP.optimizer_with_attributes(
        CS.Optimizer,
        "logging" => [],
        "seed" => 4,
        "branch_strategy" => branch_strategy,
    )
end
cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
function CSCbcJuMPTestOptimizer(; branch_strategy = :Auto)
    JuMP.optimizer_with_attributes(
        CS.Optimizer,
        "logging" => [],
        "lp_optimizer" => cbc_optimizer,
        "seed" => 2,
        "branch_strategy" => branch_strategy,
    )
end

macro test_macro_throws(errortype, m)
    # See https://discourse.julialang.org/t/test-throws-with-macros-after-pr-23533/5878
    :(@test_throws $(esc(errortype)) try
        @eval $m
    catch err
        throw(err.error)
    end)
end

"""
    test_string(arr::AbstractArray)

In Julia 1.0 string(arr) starts with Array{Int, 1} or something. In 1.5 it doesn't.
This function removes Array{...} and starts with `[`. Additionally all white spaces are removed.
"""
function test_string(arr::AbstractArray)
    s = string(arr)
    if startswith(s, "Array")
        p = first(findfirst("[", s))
        s = s[p:end]
    end
    s = replace(s, " " => "")
    return s
end

test_stime = time()

include("general.jl")
include("sudoku_fcts.jl")

include("docs.jl")
include("fcts.jl")
include("unit/index.jl")
include("options.jl")
include("moi.jl")
include("scheduling.jl")
include("constraints/table.jl")
include("constraints/indicator.jl")
include("constraints/reified.jl")
include("constraints/equal_to.jl")

include("lp_solver.jl")

include("steiner.jl")
include("stable_set.jl")
include("small_special.jl")
include("maximum_weight_matching.jl")
include("small_eq_sum_real.jl")
include("sudoku.jl")
include("str8ts.jl")
include("eternity.jl")
include("killer_sudoku.jl")
include("graph_color.jl")
println("Time for all tests $(time()-test_stime)")


text:
using Makie
using MakieLayout

include("sirmodel.jl")

R0=2.0
N=1000
γ=0.1
β=γ*R0/N

sol = solveSIR(N, β, γ)

rep = calcEffectiveReproduction(sol)

scene, layout = layoutscene(resolution=(1200,900))

xticks=0:10:160

ax = layout[1,1] = LAxis(scene, title="SIR Model,N=$N,β=$β,γ=$γ", xticks=xticks)

ls = [lines!(ax, sol.t, sol[1,:], color=:blue),
      lines!(ax, sol.t, sol[2,:], color=:red),
      lines!(ax, sol.t, sol[3,:], color=:green)]

layout[1,1] = LLegend(scene, ls, ["Susceptible", "Infectious", "Removed"],
                      width=Auto(), height=Auto(), halign=:left, valign=:center,
                      tellwidth=false, tellheight=false, margin=(10,10,10,10))

bx = layout[2,1] = LAxis(scene, title="Effective Reproduce. \0=$R0",xticks=xticks)
lines!(bx, rep[2], rep[1], color=:blue)

#display(scene)
#readline()

save("plot.png",scene)

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
# https://github.com/SciML/DiffEqFlux.jl/issues/432#issuecomment-708079051
function destructure(m; cache = IdDict())
  xs = Zygote.Buffer([])
  Flux.fmap(m) do x
    if x isa AbstractArray
      push!(xs, x)
    else
      cache[x] = x
    end
    return x
  end
  return vcat(vec.(copy(xs))...), p -> _restructure(m, p, cache = cache)
end
function _restructure(m, xs; cache = IdDict())
  i = 0
  Flux.fmap(m) do x
    x isa AbstractArray || return cache[x]
    x = reshape(xs[i.+(1:length(x))], size(x))
    i += length(x)
    return x
  end
end

paramlength(m::Flux.Dense) = length(m.weight) + length(m.bias)

struct Broadcaster{M,P}
  model::M
  p::P
  paramlength::Int
end
function Broadcaster(model)
  p = initial_params(model)
  paramlength = length(p)
  Broadcaster(model,p,paramlength)
end
(m::Broadcaster)(xs, p) = [m.model(x, p) for x in xs] # mapfoldl(x -> m.model(x, p), vcat, xs)

Base.show(io::IO, m::Broadcaster) = print(io, "Broadcaster(", m.model, ")")
initial_params(m::Broadcaster) = m.p
paramlength(m::Broadcaster) = m.paramlength

# Flux.@functor Broadcaster #(model,)
Flux.trainable(m::Broadcaster) = (m.model,)

struct Mapper{V}
  W::V
  b::V
  p::V
  paramlength::Int
end
function Mapper(in::Integer)
  W = ones(Float32,in)
  b = fill(0.00001f0,in)
  p = vcat(W,b)
  Mapper(W, b, p, length(p))
end
# function (m::Mapper{<:AbstractVector{T}})(x::Matrix{T}, p=m.p) where T
function (m::Mapper)(x, p=m.p)
  Wl = size(m.W,1)
  W = @view p[1 : Wl]
  b = @view p[Wl + 1 : end]
  W .* x .+ b
end
Base.show(io::IO, m::Mapper) = print(io, "Mapper(", length(m.W), ")")
initial_params(m::Mapper) = m.p
paramlength(m::Mapper) = m.paramlength

Flux.@functor Mapper (p,)
Flux.trainable(m::Mapper) = (m.p,)

struct FluxLayerWrapper{FL,P,RE}
  layer::FL
  p::P
  re::RE
  paramlength::Int
end
function FluxLayerWrapper(layer)
  p, re = Flux.destructure(layer)
  FluxLayerWrapper(layer, p, re, length(p))
end
(m::FluxLayerWrapper)(x, p) = m.re(p)(x)
Base.show(io::IO, m::FluxLayerWrapper) = print(io, "FluxLayerWrapper(", layer, ")")
initial_params(m::FluxLayerWrapper) = m.p
paramlength(m::FluxLayerWrapper) = m.paramlength

reset_state!(m,p) = nothing

function reset_state!(m::Union{Flux.Chain, FastChain}, p)
  start_idx = 1
  for l in m.layers
    pl = paramlength(l)
    p_layer = p[start_idx:start_idx+pl-1]
    reset_state!(l, p_layer)
    start_idx += pl
  end
end

get_bounds(m) = Float32[],Float32[]
get_bounds(m::Broadcaster) = get_bounds(m.model)
get_bounds(m::FluxLayerWrapper) = get_bounds(m.layer)

function get_bounds(m::Mapper)
  lb = vcat(fill(-10.1, length(m.W)),
            fill(-10.1, length(m.b))) |> f32
  ub = vcat(fill(-10.1, length(m.W)),
            fill(-10.1, length(m.b))) |> f32
  lb, ub
end

function get_bounds(m::Union{Flux.Chain, FastChain})
  lb = vcat([get_bounds(layer)[1] for layer in m.layers]...)
  ub = vcat([get_bounds(layer)[2] for layer in m.layers]...)
  lb, ub
end

function get_bounds(m::FastDense)
  lb = vcat(fill(-10.1, m.out*m.in),
            fill(-10.1, m.out)) |> f32
  ub = vcat(fill(10.1, m.out*m.in),
            fill(10.1, m.out)) |> f32
  lb, ub
end

function get_bounds(m::Flux.Dense)
  lb = vcat(fill(-10.1, length(m.weight)),
            fill(-10.1, length(m.bias))) |> f32
  ub = vcat(fill(10.1, length(m.weight)),
            fill(10.1, length(m.bias))) |> f32
  lb, ub
end


text:
struct Transvection <: GSymbol
    id::Symbol
    i::UInt16
    j::UInt16
    inv::Bool

    function Transvection(id::Symbol, i::Integer, j::Integer, inv = false)
        @assert id in (:ϱ, :λ)
        return new(id, i, j, inv)
    end
end

ϱ(i, j) = Transvection(:ϱ, i, j)
λ(i, j) = Transvection(:λ, i, j)

function Base.show(io::IO, t::Transvection)
    id = if t.id === :ϱ
        'ϱ'
    else # if t.id === :λ
        'λ'
    end
    print(io, id, subscriptify(t.i), '.', subscriptify(t.j))
    t.inv && print(io, "^-1")
end

Base.inv(t::Transvection) = Transvection(t.id, t.i, t.j, !t.inv)

Base.:(==)(t::Transvection, s::Transvection) =
    t.id === s.id && t.i == s.i && t.j == s.j && t.inv == s.inv

Base.hash(t::Transvection, h::UInt) = hash(hash(t.id, hash(t.i)), hash(t.j, hash(t.inv, h)))

Base.@propagate_inbounds @inline function evaluate!(
    v::NTuple{T,N},
    t::Transvection,
    tmp = one(first(v)),
) where {T,N}
    i, j = t.i, t.j
    @assert 1 ≤ i ≤ length(v) && 1 ≤ j ≤ length(v)

    A = alphabet(parent(first(v)))

    @inbounds begin
        if t.id === :ϱ
            if !t.inv
                append!(word(v[i]), word(v[j]))
            else
                # append!(word(v[i]), inv(A, word(v[j])))
                for l in Iterators.reverse(word(v[j]))
                    push!(word(v[i]), inv(A, l))
                end
            end
        else # if t.id === :λ
            if !t.inv
                # prepend!(word(v[i]), word(v[j]))
                for l in Iterators.reverse(word(v[j]))
                    pushfirst!(word(v[i]), l)
                end
            else
                # prepend!(word(v[i]), inv(A, word(v[j])))
                for l in word(v[j])
                    pushfirst!(word(v[i]), inv(A, l))
                end
            end
        end

        _setnormalform!(v[i], false)
        _setvalidhash!(v[i], false)
    end
    normalform!(tmp, v[i])
    copyto!(v[i], tmp)

    return v
end

struct PermRightAut <: GSymbol
    perm::Vector{UInt8}

    function PermRightAut(p::AbstractVector{<:Integer})
        @assert sort(p) == 1:length(p)
        return new(p)
    end
end

function Base.show(io::IO, p::PermRightAut)
    print(io, 'σ')
    join(io, (subscriptify(Int(i)) for i in p.perm))
end

Base.inv(p::PermRightAut) = PermRightAut(invperm(p.perm))

Base.:(==)(p::PermRightAut, q::PermRightAut) = p.perm == q.perm
Base.hash(p::PermRightAut, h::UInt) = hash(p.perm, hash(PermRightAut, h))

evaluate!(v::NTuple{T,N}, p::PermRightAut, tmp = nothing) where {T,N} = v[p.perm]


text:
using UserNSSandbox_jll

# Use User Namespaces to provide isolation on Linux hosts whose kernels support it
export UserNamespacesExecutor, UnprivilegedUserNamespacesExecutor, PrivilegedUserNamespacesExecutor

abstract type UserNamespacesExecutor <: SandboxExecutor; end

# A version of `chmod()` that hides all of its errors.
function chmod_recursive(root::String, perms, use_sudo::Bool)
    files = String[]
    try
        files = readdir(root)
    catch e
        if !isa(e, Base.IOError)
            rethrow(e)
        end
    end
    for f in files
        path = joinpath(root, f)
        try
            if use_sudo
                run(`$(sudo_cmd()) chmod $(string(perms, base=8)) $(path)`)
            else
                chmod(path, perms)
            end
        catch e
            if !isa(e, Base.IOError)
                rethrow(e)
            end
        end
        if isdir(path)
            chmod_recursive(path, perms, use_sudo)
        end
    end
end

function cleanup(exe::UserNamespacesExecutor)
    if exe.persistence_dir !== nothing && isdir(exe.persistence_dir)
        # Because a lot of these files are unreadable, we must `chmod +r` them before deleting
        chmod_recursive(exe.persistence_dir, 0o777, isa(exe, PrivilegedUserNamespacesExecutor))
        try
            rm(exe.persistence_dir; force=true, recursive=true)
        catch
        end
    end
end

# Because we can run in "privileged" or "unprivileged" mode, let's treat
# these as two separate, but very similar, executors.
mutable struct UnprivilegedUserNamespacesExecutor <: UserNamespacesExecutor
    persistence_dir::Union{String,Nothing}
    UnprivilegedUserNamespacesExecutor() = new(nothing)
end
mutable struct PrivilegedUserNamespacesExecutor <: UserNamespacesExecutor
    persistence_dir::Union{String,Nothing}
    PrivilegedUserNamespacesExecutor() = new(nothing)
end

Base.show(io::IO, exe::UnprivilegedUserNamespacesExecutor) = write(io, "Unprivileged User Namespaces Executor")
Base.show(io::IO, exe::PrivilegedUserNamespacesExecutor) = write(io, "Privileged User Namespaces Executor")

function executor_available(::Type{T}; verbose::Bool=false) where {T <: UserNamespacesExecutor}
    # If we're on a platform that doesn't even have `sandbox` available, return false
    if !UserNSSandbox_jll.is_available()
        return false
    end
    return with_executor(T) do exe
        return check_kernel_version(;verbose) &&
               check_overlayfs_loaded(;verbose) &&
               probe_executor(exe; test_read_only_map=true, test_read_write_map=true, verbose)
    end
end

function check_kernel_version(;verbose::Bool = false)
    # Don't bother to do anything on non-Linux
    if !Sys.islinux()
        return false
    end
    kernel_version = get_kernel_version()

    # If we were unable to parse any part of the version number, then warn and exit.
    if kernel_version === nothing
        @warn("Unable to check version number")
        return false
    end

    # Otherwise, we have a kernel version and if it's too old, we should freak out.
    if kernel_version < v"3.18"
        @warn("Kernel version too old: detected $(kernel_version), need at least 3.18!")
        return false
    end

    if verbose
        @info("Parsed kernel version \"$(kernel_version)\"")
    end
    return true
end

function check_overlayfs_loaded(;verbose::Bool = false)
    if !Sys.islinux()
        return false
    end

    # If the user has disabled this check, return `true`
    if parse(Bool, get(ENV, "SANDBOX_SKIP_OVERLAYFS_CHECK", "false"))
        return true
    end

    mods = get_loaded_modules()
    if verbose
        @info("Found $(length(mods)) loaded modules")
    end

    filter!(mods) do (name, size, count, deps, state, addr)
        return name == "overlay"
    end
    if isempty(mods)
        @warn("Could not find loaded `overlay` module, try `sudo modprobe overlay` or export SANDBOX_SKIP_OVERLAYFS_CHECK=true to disable this check!")
        return false
    end

    if verbose
        @info("Found loaded `overlay` module")
    end
    return true
end

function build_executor_command(exe::UserNamespacesExecutor, config::SandboxConfig, user_cmd::Cmd)
    # While we would usually prefer to use the `executable_product()` function to get a
    # `Cmd` object that has all of the `PATH` and `LD_LIBRARY_PATH` environment variables
    # set properly so that the executable product can be run, we are careful to ensure
    # that `sandbox` has no dependencies (as much as that is possible).
    cmd_string = String[UserNSSandbox_jll.sandbox_path]

    # Enable verbose mode on the sandbox wrapper itself
    if config.verbose
        push!(cmd_string, "--verbose")
    end

    # Extract the rootfs, as it's treated specially
    append!(cmd_string, ["--rootfs", config.read_only_maps["/"]])

    # Add our `--cd` command
    append!(cmd_string, ["--cd", config.pwd])

    # Add in read-only mappings (skipping the rootfs)
    for (dst, src) in config.read_only_maps
        if dst == "/"
            continue
        end
        append!(cmd_string, ["--map", "$(src):$(dst)"])
    end

    # Add in read-write mappings
    for (dst, src) in config.read_write_maps
        append!(cmd_string, ["--workspace", "$(src):$(dst)"])
    end

    # Add in entrypoint, if it is set
    if config.entrypoint !== nothing
        append!(cmd_string, ["--entrypoint", config.entrypoint])
    end

    # If we have a `--persist` argument, check to see if we already have a persistence_dir
    # setup, if we do not, create a temporary directory and set it into our executor
    if config.persist
        if exe.persistence_dir === nothing
            persist_parent_dir = get(ENV, "SANDBOX_PERSISTENCE_DIR", tempdir())
            mkpath(persist_parent_dir)
            exe.persistence_dir = mktempdir(persist_parent_dir)
        end
        append!(cmd_string, ["--persist", exe.persistence_dir])
    end

    # For each platform requested by `multiarch`, ensure its matching interpreter is registered.
    register_requested_formats!(config.multiarch_formats; verbose=config.verbose)

    # Set the user and group, if requested
    append!(cmd_string, ["--uid", string(config.uid), "--gid", string(config.gid)])

    # Set the custom tmpfs_size, if requested
    if config.tmpfs_size !== nothing
        append!(cmd_string, ["--tmpfs-size", config.tmpfs_size])
    end

    # If we're running in privileged mode, we need to add `sudo` (or `su`, if `sudo` doesn't exist)
    if isa(exe, PrivilegedUserNamespacesExecutor)
        # Next, prefer `sudo`, but allow fallback to `su`. Also, force-set our
        # environmental mappings with sudo, because many of these are often  lost
        # and forgotten due to `sudo` restrictions on setting `LD_LIBRARY_PATH`, etc...
        if get(sudo_cmd(), 1, "") == "sudo"
            sudo_envs = vcat([["-E", "$k=$(config.env[k])"] for k in keys(config.env)]...)
            if user_cmd.env !== nothing
                append!(sudo_envs, vcat([["-E", pair] for pair in user_cmd.env]...))
            end
            prepend!(cmd_string, String[sudo_cmd()..., sudo_envs...])
        else
            prepend!(cmd_string, sudo_cmd())
        end
    end

    # Finally, append the user-requested command string
    push!(cmd_string, "--")
    append!(cmd_string, user_cmd.exec)

    # Construct a `Cmd` object off of those, with the SandboxConfig's env (if this is an unprivileged runner):
    sandbox_cmd = setenv(Cmd(cmd_string), String[])
    if isa(exe, UnprivilegedUserNamespacesExecutor)
        sandbox_cmd = setenv(sandbox_cmd, config.env)

        # If the user has provided an environment with their command, merge that in as well
        if user_cmd.env !== nothing
            sandbox_cmd = addenv(sandbox_cmd, user_cmd.env)
        end
    end

    # If the user has asked that this command be allowed to fail silently, pass that on
    if user_cmd.ignorestatus
        sandbox_cmd = ignorestatus(sandbox_cmd)
    end

    return sandbox_cmd
end


text:
using Printf
using LinearAlgebra
using SparseArrays

function tensor_grid(grids)
    n = 1
    for g in grids
        n *= length(g)
    end
    X = zeros(length(grids), n)
    for (i,x) in enumerate(Base.product(grids...))
        X[:,i] = [x...]
    end
    return X
end

function show_array(arr, name=nothing)
    if name != nothing
        @printf "Array %s\n" name
    else
        @printf "Array\n"
    end
    show(IOContext(stdout, :compact => true), "text/plain", arr)
    @printf "\n"
end

# Test matrix generation
An = n -> spdiagm(-1 => -ones(n-1), 0 => 2*ones(n), 1 => -ones(n-1))

# order-2 Laplacian in 1, 2, or 3D
function Ad(n, d)
    I = sparse(UniformScaling{Float64}(1.0), n, n)
    if d == 1
        A = An(n)
    elseif d == 2
        A = kron(An(n), I) + kron(I, An(n))
    elseif d == 3
        A = kron(An(n), I, I) + kron(I, An(n), I) + kron(I, I, An(n))
    end
    return A
end

# order-2/4/6/8 Laplacian in 1, 2, or 3D
function Ad_o(n, d, o)
    order2 = ([1, -2, 1], [1, 0, 1])
    order4 = ([-1/12, 4/3, -5/2, 4/3, -1/12], [2, 1, 0, 1, 2])
    order6 = ([1/90, -3/20, 3/2, -49/18, 3/2, -3/20, 1/90], [3, 2, 1, 0, 1, 2, 3])
    order8 = ([-1/560, 8/315, -1/5, 8/5, -205/72, 8/5, -1/5, 8/315, -1/560], [4, 3, 2, 1, 0, 1, 2, 3, 4])
    e = ones(n)
    I = sparse(UniformScaling{Float64}(1.0), n, n)
    if o == 2
        Ao = spdiagm(([c*e[1:n-s] for (c,s) in zip(order2...)]...), (-1, 0, 1), n, n)
    elseif o == 4
        Ao = spdiagm(([c*e[1:n-s] for (c,s) in zip(order4...)]...), (-2, -1, 0, 1, 2), n, n)
    elseif o == 6
        Ao = spdiagm(([c*e[1:n-s] for (c,s) in zip(order6...)]...), (-3, -2, -1, 0, 1, 2, 3), n, n)
    elseif o == 8
        Ao = spdiagm(([c*e[1:n-s] for (c,s) in zip(order8...)]...), (-4, -3, -2, -1, 0, 1, 2, 3, 4), n, n)
    end
    if d == 1
        A = - Ao
    elseif d == 2
        A = kron(- Ao, I) + kron(I, - Ao)
    elseif d == 3
        A = kron(- Ao, I, I) + kron(I, - Ao, I) + kron(I, I, - Ao)
    end
    return A
end

function count_sorted(v::AbstractArray{V,1}) where {V}
    if length(v) == 0
        return 0
    end
    lastv = v[1]
    n = 1
    for i = 2:length(v)
        if v[i] != lastv
            n += 1
            lastv = v[i]
        end
    end
    return n
end

function sorted_to_ptr(v::AbstractArray{V}) where {V}
    l = count_sorted(v)
    r = Array{Int64,1}(undef, l+1)
    r[1] = 1
    r[l+1] = length(v)+1
    if length(v) == 0 || length(v) == 1
        return r
    end
    lastv = v[1]
    n = 1
    for i = 2:length(v)
        if v[i] != lastv
            r[n+1] = i
            n += 1
            lastv = v[i]
        end
    end
    return r
end

# [1, 2, 1] -> [1, 2, 4, 5]
function sizes_to_ptr(v::AbstractArray)
    if length(v) == 0
        return [1]
    else
        return vcat([1], 1 .+ cumsum(v))
    end
end

# Fill Adst with Asrc[n,s]
function fill_dense_block!(Asrc::SparseMatrixCSC{Tv,Ti}, Adst::StridedMatrix{Tv}, n::UnitRange{Ti}, s::UnitRange{Ti},) where {Tv, Ti <: Integer}
    for j = 1:length(s)
        k = A.colptr[s[j]]:(A.colptr[s[j]+1]-1)
        I = A.rowval[k]
        V = A.nzval[k]
        ids = searchsortedfirst(I, n[1]):searchsortedlast(I, n[end])
        for i in ids
            Adst[I[i] - n[1] + 1, j] = V[i]
        end
    end
end

# In a binary tree, returns true if (l1, s1) and (l2, s2) are on same path to the root
#			(3,1)
#	      /       \
#     (2,1)      (2,2)
#     /  \       /   \
# (1,1) (1,2) (1,3) (1,4)
function are_connected(l1, s1, l2, s2, maxLevel)
    @assert l1 > 0 && l1 <= maxLevel && s1 > 0 && s1 <= 2^(maxLevel-l1)
    @assert l2 > 0 && l2 <= maxLevel && s2 > 0 && s2 <= 2^(maxLevel-l2)
    if l1 == l2
        return s1 == s2
    end
	L = [l1, l2]
    S = [s1, s2]
    (lmax, imax) = findmax(L)
    (lmin, imin) = findmin(L)
    smax = S[imax]
    smin = S[imin]
    while lmin < lmax
        lmin += 1
        smin = div(smin+1,2)
    end
    smin == smax
end

# BinTree iterator
# Generate
# (1, 1), (1, 2), (1, 3), (1, 4)
#     (2, 1),         (2, 2)
#             (3, 1)
struct BinTreeIt
    maxLevel::Int64
end

function Base.iterate(iter::BinTreeIt, state=(0, 2^iter.maxLevel))
    (l, s) = state
    if l >= iter.maxLevel
        return nothing
    end
    if s >= 2^(iter.maxLevel-l)
        next = (l+1, 1)
    else
        next = (l, s+1)
    end
    return (next, next)
end

# Fills colptr and rowval (assumed large enough) with A[list, list], excluding the diagonal
# List is assumed to be sorted
# base-0
# Ignore the diagonal
function fill_nodiag_list_sorted!(A::SparseMatrixCSC{Tv,Ti}, colptr::Array{Tj}, rowval::Array{Tj}, list::AbstractArray{Tk}) where {Ti <: Integer, Tj <: Integer, Tk <: Integer, Tv}
    n = length(list)
    colptr[1] = 1
    if n == 0
        return
    end
    for j = 1:n
        colptr[j+1] = colptr[j] # Initially, nothing
        j_ = list[j]
        i_ = 1
        for k = A.colptr[j_]:(A.colptr[j_+1]-1)
            i = A.rowval[k]
            if i == j_ # Skip diagonal
                continue
            end
            i_ = searchsortedfirst(list, Tk(i), i_, length(list), ord) # Ok since sorted
            if i_ < n+1 && list[i_] == Tk(i) # Yes, at i_
                rowval[colptr[j+1]] = i_
                colptr[j+1] += 1
            end
        end
    end
    nz = colptr[n+1]-1
    colptr[1:n+1] .-= 1
    rowval[1:nz]  .-= 1
    return
end


text:
#EdgeExistence problem in hackerearth.com

path = "C:\\Users\\didier1708\\Documents\\GitHub\\Julia-playground\\EdgeExistence.txt";
r = readlines(path);
myParse = function(str::String)
    a = [parse(Int, word) for word in split(str)];
    return a;
end
rs = map(myParse, r);
k = rs[1][1];
n = rs[1][2];
Edges = rs[2:1+n];
q = rs[n+2][1];
Queries = rs[n+3:n+2+q];
spit = function(query)
    answer = ""
    (query in Edges) ? answer="YES" : answer ="NO";
    return answer
end
Answers = map(spit, Queries);
map(println, Answers);

max_stars_repo_name: dan-seol/Julia-playground

max_issues_repo_name: dan-seol/Julia-playground

max_forks_repo_name: dan-seol/Julia-playground


text:
export CP

mutable struct CP <: PopulationModel
    k::Int
    cp::Union{CANDECOMP{Float64,3},Float64}
	function CP(k, cp=NaN)
		new(k,cp)
	end
end

function fit!(model::CP, Z::Array{Float64,3})
    model.cp = candecomp(Z, model.k, ([randn(x,model.k) for x in size(Z)]...,),
    compute_error=true, method=:ALS, verbose=false, maxiter=1000)
end

weights(model::CP) = model.cp.factors[1]

factors(model::CP, dim=0) = dim == 0 ? model.cp.factors : model.cp.factors[dim]


text:
using FilePaths

JS_DIR = "js"

begin
    download("https://cdn.jsdelivr.net/npm/vega@5", joinpath(JS_DIR, "vega.min.js"))
    download("https://cdn.jsdelivr.net/npm/vega-lite@4", joinpath(JS_DIR, "vega-lite.min.js"))
    download("https://cdn.jsdelivr.net/npm/vega-embed@6", joinpath(JS_DIR, "vega-embed.min.js"))
end

max_stars_repo_name: benjavicente/covid-maps

max_issues_repo_name: benjavicente/covid-maps

max_forks_repo_name: benjavicente/covid-maps


text:
using Revise
using Random
using BenchmarkTools
using POMDPs
using POMDPModelTools
# using CuArrays
using Flux
using DeepQLearning
include("test/test_env.jl")

mdp = TestMDP((5,5), 4, 6)
# mdp = SimpleGridWorld()
rng = MersenneTwister(1)
mdp = TestMDP((5,5), 1, 6)
model = Chain(x-> flattenbatch(x), Dense(25, 32), Dense(32,32), Dense(32, 32, relu), Dense(32, length(actions(mdp))))

solver = DeepQLearningSolver(batch_size = 128, eval_freq = 10_000, save_freq=10_000, qnetwork = model, prioritized_replay=true,
                             max_steps=1000, learning_rate=0.001,train_start=500,log_freq=100000,
                             recurrence=false,trace_length=10, double_q=true, dueling=false, rng=rng, verbose=false)

@btime policy = solve($solver, $mdp)

policy = solve(solver, mdp)

env = MDPEnvironment(mdp)
o = reset!(env)

using RLInterface
using LinearAlgebra

env = MDPEnvironment(mdp)
replay = DeepQLearning.initialize_replay_buffer(solver, env)
active_q = solver.qnetwork
policy = NNPolicy(env.problem, active_q, ordered_actions(env.problem), length(obs_dimensions(env)))
active_q = getnetwork(policy)
target_q = deepcopy(active_q)

s_batch, a_batch, r_batch, sp_batch, done_batch, indices, importance_weights = DeepQLearning.sample(replay)

q_values = active_q(s_batch)
best_a = Flux.onecold(q_values)

argmax(q_values, dims=1)

best_a = [CartesianIndex(Flux.argmax(q_values[:, i]), i) for i=1:solver.batch_size]
q_values[best_a] == dropdims(maximum(q_values, dims=1), dims=1)

a_batch
q_values[a_batch]

@btime q_sa = diag(view($q_values, $a_batch, :));
@btime q_sa = q_values[a, :]
q_sa2 = q_values[:][a_batch]

@btime getindex.(Ref($q_values), $a_batch, 1:$solver.batch_size);
a_batch = CartesianIndex.(a_batch, 1:solver.batch_size)

A = rand(4, 32)
a = rand(1:4, 32) # indices of the first column of A
@btime diag(view($A, $a, :)); # first way
@btime getindex.(Ref($A), $a, 1:size($A, 2));

function test(active_q, target_q)
    p = Flux.params(active_q)
    loss_val = nothing
    td_err = nothing
    gs = Flux.gradient(p) do
        # compute q_vals and targets
        q_values = active_q(s_batch)
        q_sa = diag(view(q_values, a_batch, :))
        γ = discount(env.problem)
        q_sp_max = dropdims(maximum(target_q(sp_batch), dims=1), dims=1)
        q_targets = r_batch .+ (1f0 .- done_batch) .* γ .* q_sp_max
        td_err = q_sa .- q_targets
        loss_val = mean(huber_loss, importance_weights.*q_sa)
        loss_val
    end
    # for (i, g) in gs.grads
    #     @show g
    #     @show g.contents
    # end
    DeepQLearning.globalnorm(p, gs)
    return loss_val, td_err
end

gs.grads

@show [a for (i, a) in gs.grads]

# benchmark
# standard dqn 650ms 282MB
# double dqn   706ms 329MB
# zygote dqn   764.5 175MB
# ddqn zygote  832.5 184MB

using Profile
using ProfileView
@profile 1+1
Profile.clear()

@profile solve(solver, mdp)

ProfileView.view()

### Try on SubHunt

using Revise
using POMDPs
using SubHunt
using RLInterface
using DeepQLearning
using Flux

solver = DeepQLearningSolver(qnetwork= Chain(Dense(8, 32, relu), Dense(32,32,relu), Dense(32, 6)),
                             max_steps=100_000)
solve(solver, SubHuntPOMDP())

### get q_sa
na = 4
bs = 32
tl = 100

s_batch = [rand(2, bs) for i=1:solver.trace_length]
a_batch = [rand(1:na, bs) for i=1:solver.trace_length]
q_values = model.(s_batch) # vector of size trace_length nactions x batch_size
q_sa = [diag(view(q_values[i], a_batch[i], :)) for i=1:solver.trace_length]

using Flux: onehotbatch

onehotbatch(a_batch[1], 1:4)

a_batch[1]

sbatch = [rand(2, 32) for i=1:tl]
q_values = model.(sbatch)
qsa = getindex.(q_values, a_batch)

using POMDPModels
using RLInterface

env = POMDPEnvironment(TigerPOMDP())

function simulate(env::AbstractEnvironment, nsteps::Int = 10)
    done = false
    r_tot = 0.0
    step = 1
    o = reset!(env)
    while !done && step <= nsteps
        action = sample_action(env) # take random action
        obs, rew, done, info = step!(env, action)
        @show obs, rew, done, info
        r_tot += rew
        step += 1
    end
    return r_tot
end

@show simulate(env)


text:
module DensityEstimationML

using StatsBase
using Distributions
using TensorFlow
using StaticArrays

using Reexport

include("./prob_utils.jl")
include("estimator.jl")
include("network.jl")
include("./GenerateDatasets.jl")

@reexport using .GenerateDatasets

end # module


text:
module FilesLH

using ArgCheck, DocStringExtensions, Parameters
import Base.show
import CommonLH: ask_yes_no

# Directories
export clear_directory, is_dir_empty, files_in_dir, files_not_in_dir2, find_common_base_dir, dir_diff_report, make_dir, right_dirs
# File names
export add_extension, change_extension, has_extension, extensions_equal, filesep, paths_equal
# File transfer
export remote_copy, rsync_command, rsync_dir, git_upload_dir, make_remote_dir, deploy_docs
# Computers
export Computer, computerLH, computerLongleaf
export home_dir, run_local, is_remote, common_base_dir, remote_and_local_path
export ComputerList, defaultCompList
export add_computer, remove_computer, current_computer, get_computer
# Slurm
export sbatch_commands

include("computers.jl")
include("directories.jl")
include("file_names.jl")
include("file_transfer.jl")
include("sbatch.jl")

end # module