text:
se_kernel, se_n_hyper = include("se_kernel.jl")
scale_kernel, scale_n_hyper = include("scale_kernel.jl")

"""
Custom function created by combining kernels from kernel_coder().
hyperparameters == ["λ, shift, σ1, σ2"]
"""
function custom_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    outputs::Vector{<:Integer}=[1,1])

    @assert length(hyperparameters)==4 "hyperparameters is the wrong length"
    dorder_len = 6
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"
    @assert minimum(dorder) >= 0
    @assert maximum(outputs) < 3
    @assert minimum(outputs) > 0

    λ = hyperparameters[1]
    shift = hyperparameters[2]
    σ1 = hyperparameters[3]
    σ2 = hyperparameters[4]

    # dλ = dorder[3]
    # dshift = dorder[4]
    dσ1 = dorder[5]
    dσ2 = dorder[6]

    if outputs==[1,1]
        return se_kernel([λ], δ, dorder[1:3]) + scale_kernel([σ1], δ, append!(dorder[1:2],[dσ2])) * Int(δ == 0)
    end

    if outputs in [[1,2],[2,1]]
        return se_kernel([λ, -shift], δ, dorder[1:4]; shift_ind=2)
    end

    if outputs==[2,2]
        return se_kernel([λ], δ, dorder[1:3]) + scale_kernel([σ2], δ, append!(dorder[1:2],[dσ2])) * Int(δ == 0)
    end

end

return custom_kernel, 4  # the function handle and the number of kernel hyperparameters

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
# ---
# title: 55. Jump Game
# id: problem55
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Greedy
# link: <https://leetcode.com/problems/jump-game/description/>
# hidden: true
# ---
#
# Given an array of non-negative integers, you are initially positioned at the
# first index of the array.
#
# Each element in the array represents your maximum jump length at that
# position.
#
# Determine if you are able to reach the last index.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,3,1,1,4]
#     Output: true
#     Explanation: Jump 1 step from index 0 to 1, then 3 steps to the last index.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,2,1,0,4]
#     Output: false
#     Explanation: You will always arrive at index 3 no matter what. Its maximum jump length is 0, which makes it impossible to reach the last index.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 3 * 10^4`
#   * `0 <= nums[i][j] <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
# online methods
module olm
dir = @__DIR__

push!(LOAD_PATH, dir * "/../sim")
push!(LOAD_PATH, dir * "/../dis")
using sim
using dis

const olm_dt = 1.0 # time increment for planning step
const olm_h = 1e-1 # time increment for numerical integration
#const olm_gamma = 0.9999 # discount factor
const olm_gamma = 0.999 # discount factor

mutable struct Tree
  value::Any
  next::AbstractArray{Tree, 1}

  Tree() = new(nothing, Tree[])
  Tree(value) = new(value, Tree[])
end

mutable struct GroupTree
  value::Any
  next::AbstractArray{Array{GroupTree, 1}, 1}

  GroupTree() = new(nothing, Array{GroupTree, 1}[])
  GroupTree(value) = new(value, Array{GroupTree, 1}[])
end

include(dir * "/olm_mcts.jl")
include(dir * "/olm_fwds.jl")
include(dir * "/olm_pofs.jl")
include(dir * "/olm_pomc.jl")

end


text:
using IOLA
using IOLA.Utils
using HDF5

max_points = 100
res_size = 2048
codecs_num = length(instances(Codec.CodecType))

parameters = Dict{Codec.CodecType, Tuple{Function,Int64,Int64}}()

function to_float64(arr::AbstractArray{T}) where T<:Integer
    return convert.(Float64, arr)./typemax(T)
end

h5open("dataset.h5") do dataset
    data = dataset["data"]
    siglen = size(data,1)
    po = max_points + 1
    @label loop1
    po -= 1
    prev_rpo = 0
    for codec in instances(Codec.CodecType)
        par = Codec.getparams(codec)
        len = (siglen-par.hop)/po
        len = round(len/par.length)*par.length
        rpo = div(siglen-par.hop,len)
        prev_rpo = prev_rpo == 0 ? rpo : prev_rpo
        if  prev_rpo != rpo
            @goto loop1
        end
    end
    reslen = Int(prev_rpo)
    println("Number of points: $reslen")

    for codec in instances(Codec.CodecType)
        par = Codec.getparams(codec)
        len = (siglen-par.hop)/po
        len = round(len/par.length)*par.length
        push!(parameters, codec => (Codec.gettransform(par), len, par.hop))
    end
    h5open("results.h5", "w") do results
        C = d_create(results, "C", datatype(Float64),
                     dataspace(codecs_num,res_size),
                     "chunk", (codecs_num,1))
        KD = d_create(results, "KD", datatype(Float64),
                      dataspace(reslen,2,codecs_num,res_size),
                      "chunk", (reslen,2,codecs_num,1))
        p = d_create(results, "p", datatype(Int64),
                     dataspace(codecs_num,res_size),
                     "chunk", (codecs_num,1))

        C_tmp = zeros(Float64,codecs_num)
        KD_tmp = zeros(Float64,reslen,2,codecs_num)
        p_tmp = zeros(Int64,codecs_num)

        for data_i in 1:res_size
            print("Progress: $data_i/$res_size\r")
            fdata = to_float64(view(data[:,data_i],:,1))
            for (i,codec) in enumerate(instances(Codec.CodecType))
                KD_tmp[:,:,i] = analyze(fdata, parameters[codec]...)
                p_tmp[i] = findperiod(KD_tmp[:,2,i], 3:parameters[codec][3])
                C_tmp[i] = radiusofmean(KD_tmp[:,1,i], 2π*mod.(KD_tmp[:,2,i],
                                                      p_tmp[i])./p_tmp[i])
            end
            C[:,data_i] = C_tmp
            KD[:,:,:,data_i] = KD_tmp
            p[:,data_i] = p_tmp
        end
    end
end

max_stars_repo_name: Ashymad/IOLA

max_issues_repo_name: Ashymad/IOLA

max_forks_repo_name: Ashymad/IOLA


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

mutable struct P4estMPIInterfaceContainer{NDIMS, uEltype<:Real, NDIMSP2} <: AbstractContainer
  u                ::Array{uEltype, NDIMSP2}       # [primary/secondary, variable, i, j, interface]
  local_element_ids::Vector{Int}                   # [interface]
  node_indices     ::Vector{NTuple{NDIMS, Symbol}} # [interface]
  local_sides      ::Vector{Int}                   # [interface]

  # internal `resize!`able storage
  _u           ::Vector{uEltype}
end

@inline nmpiinterfaces(interfaces::P4estMPIInterfaceContainer) = length(interfaces.local_sides)
@inline Base.ndims(::P4estMPIInterfaceContainer{NDIMS}) where NDIMS = NDIMS

function Base.resize!(mpi_interfaces::P4estMPIInterfaceContainer, capacity)
  @unpack _u, local_element_ids, node_indices, local_sides = mpi_interfaces

  n_dims = ndims(mpi_interfaces)
  n_nodes = size(interfaces.u, 3)
  n_variables = size(interfaces.u, 2)

  resize!(_u, 2 * n_variables * n_nodes^(n_dims-1) * capacity)
  interfaces.u = unsafe_wrap(Array, pointer(_u),
    (2, n_variables, ntuple(_ -> n_nodes, n_dims-1)..., capacity))

  resize!(local_element_ids, capacity)

  resize!(node_indices, capacity)

  resize!(local_sides, capacity)

  return nothing
end

# Create MPI interface container and initialize interface data
function init_mpi_interfaces(mesh::ParallelP4estMesh, equations, basis, elements)
  NDIMS = ndims(elements)
  uEltype = eltype(elements)

  # Initialize container
  n_mpi_interfaces = count_required_surfaces(mesh).mpi_interfaces

  _u = Vector{uEltype}(undef, 2 * nvariables(equations) * nnodes(basis)^(NDIMS-1) * n_mpi_interfaces)
  u = unsafe_wrap(Array, pointer(_u),
    (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS-1)..., n_mpi_interfaces))

  local_element_ids = Vector{Int}(undef, n_mpi_interfaces)

  node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, n_mpi_interfaces)

  local_sides = Vector{Int}(undef, n_mpi_interfaces)

  mpi_interfaces = P4estMPIInterfaceContainer{NDIMS, uEltype, NDIMS+2}(
    u, local_element_ids, node_indices, local_sides, _u)

  init_mpi_interfaces!(mpi_interfaces, mesh)

  return mpi_interfaces
end

function init_mpi_interfaces!(mpi_interfaces, mesh::ParallelP4estMesh)
  init_surfaces!(nothing, nothing, nothing, mpi_interfaces, mesh)

  return mpi_interfaces
end

# Overload init! function for boundaries and regular interfaces since it must call the appropriate
# init_surfaces! function for parallel p4est meshes
function init_interfaces!(interfaces, mesh::ParallelP4estMesh)
  init_surfaces!(interfaces, nothing, nothing, nothing, mesh)

  return interfaces
end

function init_boundaries!(boundaries, mesh::ParallelP4estMesh)
  init_surfaces!(nothing, nothing, boundaries, nothing, mesh)

  return boundaries
end

function reinitialize_containers!(mesh::ParallelP4estMesh, equations, dg::DGSEM, cache)
  # Re-initialize elements container
  @unpack elements = cache
  resize!(elements, ncells(mesh))
  init_elements!(elements, mesh, dg.basis)

  required = count_required_surfaces(mesh)

  # resize interfaces container
  @unpack interfaces = cache
  resize!(interfaces, required.interfaces)

  # resize boundaries container
  @unpack boundaries = cache
  resize!(boundaries, required.boundaries)

  # resize mortars container
  @unpack mortars = cache
  resize!(mortars, required.mortars)

  # resize mpi_interfaces container
  @unpack mpi_interfaces = cache
  resize!(mpi_interfaces, required.mpi_interfaces)

  # re-initialize containers together to reduce
  # the number of iterations over the mesh in p4est
  init_surfaces!(interfaces, mortars, boundaries, mpi_interfaces, mesh)
end

# A helper struct used in initialization methods below
mutable struct ParallelInitSurfacesIterFaceUserData{Interfaces, Mortars, Boundaries, MPIInterfaces, Mesh}
  interfaces      ::Interfaces
  interface_id    ::Int
  mortars         ::Mortars
  mortar_id       ::Int
  boundaries      ::Boundaries
  boundary_id     ::Int
  mpi_interfaces  ::MPIInterfaces
  mpi_interface_id::Int
  mesh            ::Mesh
end

function ParallelInitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mpi_interfaces, mesh)
  return ParallelInitSurfacesIterFaceUserData{
    typeof(interfaces), typeof(mortars), typeof(boundaries), typeof(mpi_interfaces), typeof(mesh)}(
      interfaces, 1, mortars, 1, boundaries, 1, mpi_interfaces, 1, mesh)
end

function init_surfaces_iter_face_parallel(info, user_data)
  # Unpack user_data
  data = unsafe_pointer_to_objref(Ptr{ParallelInitSurfacesIterFaceUserData}(user_data))

  # Function barrier because the unpacked user_data above is type-unstable
  init_surfaces_iter_face_inner(info, data)
end

# 2D
cfunction(::typeof(init_surfaces_iter_face_parallel), ::Val{2}) = @cfunction(init_surfaces_iter_face_parallel, Cvoid, (Ptr{p4est_iter_face_info_t}, Ptr{Cvoid}))

# Function barrier for type stability, overload for parallel P4estMesh
function init_surfaces_iter_face_inner(info, user_data::ParallelInitSurfacesIterFaceUserData)
  @unpack interfaces, boundaries, mpi_interfaces = user_data

  if info.sides.elem_count == 2
    # Two neighboring elements => Interface or mortar

    # Extract surface data
    sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))

    if sides[1].is_hanging == 0 && sides[2].is_hanging == 0
      # No hanging nodes => normal interface or MPI interface
      if sides[1].is.full.is_ghost == 1 || sides[2].is.full.is_ghost == 1 # remote side => MPI interface
        if mpi_interfaces !== nothing
          init_mpi_interfaces_iter_face_inner(info, sides, user_data)
        end
      else
        if interfaces !== nothing
          init_interfaces_iter_face_inner(info, sides, user_data)
        end
      end
    else
      # Hanging nodes => mortar
      error("ParallelP4estMesh does not support non-conforming meshes.")
    end
  elseif info.sides.elem_count == 1
    # One neighboring elements => boundary
    if boundaries !== nothing
      init_boundaries_iter_face_inner(info, user_data)
    end
  end

  return nothing
end

function init_surfaces!(interfaces, mortars, boundaries, mpi_interfaces, mesh::ParallelP4estMesh)
  # Let p4est iterate over all interfaces and call init_surfaces_iter_face
  iter_face_c = cfunction(init_surfaces_iter_face_parallel, Val(ndims(mesh)))
  user_data = ParallelInitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mpi_interfaces,
                                                   mesh)

  iterate_p4est(mesh.p4est, user_data; iter_face_c=iter_face_c)

  return nothing
end

# Initialization of MPI interfaces after the function barrier
function init_mpi_interfaces_iter_face_inner(info, sides, user_data)
  @unpack mpi_interfaces, mpi_interface_id, mesh = user_data
  user_data.mpi_interface_id += 1

  if sides[1].is.full.is_ghost == 1
    local_side = 2
  elseif sides[2].is.full.is_ghost == 1
    local_side = 1
  else
    error("should not happen")
  end

  # Get local tree, one-based indexing
  tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1)
  # Quadrant numbering offset of the local quadrant at this interface
  offset = tree.quadrants_offset
  tree_quad_id = sides[local_side].is.full.quadid # quadid in the local tree
  # ID of the local neighboring quad, cumulative over local trees
  local_quad_id = offset + tree_quad_id

  # p4est uses zero-based indexing, convert to one-based indexing
  mpi_interfaces.local_element_ids[mpi_interface_id] = local_quad_id + 1
  mpi_interfaces.local_sides[mpi_interface_id] = local_side

  # Local face at which the interface lies
  local_face = sides[local_side].face

  # Save mpi_interfaces.node_indices dimension specific in containers_[23]d_parallel.jl
  init_mpi_interface_node_indices!(mpi_interfaces, local_face, local_side, info.orientation,
                                   mpi_interface_id)

  return nothing
end

# Iterate over all interfaces and count
# - (inner) interfaces
# - mortars
# - boundaries
# - (MPI) interfaces at subdomain boundaries
# and collect the numbers in `user_data` in this order.
function count_surfaces_iter_face_parallel(info, user_data)
  if info.sides.elem_count == 2
    # Two neighboring elements => Interface or mortar

    # Extract surface data
    sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))

    if sides[1].is_hanging == 0 && sides[2].is_hanging == 0
      # No hanging nodes => normal interface or MPI interface
      if sides[1].is.full.is_ghost == 1 || sides[2].is.full.is_ghost == 1 # remote side => MPI interface
        # Unpack user_data = [mpi_interface_count] and increment mpi_interface_count
        ptr = Ptr{Int}(user_data)
        id = unsafe_load(ptr, 4)
        unsafe_store!(ptr, id + 1, 4)
      else
        # Unpack user_data = [interface_count] and increment interface_count
        ptr = Ptr{Int}(user_data)
        id = unsafe_load(ptr, 1)
        unsafe_store!(ptr, id + 1, 1)
      end
    else
      # Hanging nodes => mortar
      error("ParallelP4estMesh does not support non-conforming meshes.")
    end
  elseif info.sides.elem_count == 1
    # One neighboring elements => boundary

    # Unpack user_data = [boundary_count] and increment boundary_count
    ptr = Ptr{Int}(user_data)
    id = unsafe_load(ptr, 3)
    unsafe_store!(ptr, id + 1, 3)
  end

  return nothing
end

# 2D
cfunction(::typeof(count_surfaces_iter_face_parallel), ::Val{2}) = @cfunction(count_surfaces_iter_face_parallel, Cvoid, (Ptr{p4est_iter_face_info_t}, Ptr{Cvoid}))

function count_required_surfaces(mesh::ParallelP4estMesh)
  # Let p4est iterate over all interfaces and call count_surfaces_iter_face_parallel
  iter_face_c = cfunction(count_surfaces_iter_face_parallel, Val(ndims(mesh)))

  # interfaces, mortars, boundaries, mpi_interfaces
  user_data = [0, 0, 0, 0]

  iterate_p4est(mesh.p4est, user_data; iter_face_c=iter_face_c)

  # Return counters
  return (interfaces     = user_data[1],
          mortars        = user_data[2],
          boundaries     = user_data[3],
          mpi_interfaces = user_data[4])
end

end # @muladd


text:
module evolution

export init_population, mutate, eval_chromo_similarity

using Random: rand
using StatsBase: sample, Weights

include("constants.jl")

function init_population(cols_number::Int, tabu_list::Set)::Population
    population = Population()
    while length(population) < POPULATION_SIZE
        chromo_size = rand(MIN_CHROMO_SIZE:INIT_MAX_CHROMO_SIZE)
        random_chromo = sample(1:cols_number, chromo_size, replace = false)

        chromo_signature = hash(random_chromo)
        if !(chromo_signature in tabu_list)
            push!(population, random_chromo)
            push!(tabu_list, chromo_signature)
        end
    end
    return population
end

function eval_chromo_similarity(chromo1::Chromo, chromo2::Chromo)::Float64
    length(intersect(chromo1, chromo2)) / min(length(chromo1), length(chromo2))
end

function mutation_swap(chromo::Chromo)::Chromo
    random_idx1 = rand(1:length(chromo))
    random_idx2 = rand(1:length(chromo))

    while random_idx1 == random_idx2
        random_idx2 = rand(1:length(chromo))
    end

    chromo = copy(chromo)
    chromo[random_idx1], chromo[random_idx2] = chromo[random_idx2], chromo[random_idx1]
    return chromo
end

function mutation_substitution(chromo::Chromo, cols_number::Int)::Chromo
    random_col = rand(1:cols_number)
    while random_col in chromo
        random_col = rand(1:cols_number)
    end
    random_substitution_point = rand(1:length(chromo))
    chromo = copy(chromo)
    chromo[random_substitution_point] = random_col
    return chromo
end

function mutation_insertion(chromo::Chromo, cols_number::Int)::Chromo
    random_col = rand(1:cols_number)
    while random_col in chromo
        random_col = rand(1:cols_number)
    end
    random_insertion_point = rand(1:(length(chromo) + 1))
    return insert!(copy(chromo), random_insertion_point, random_col)
end

function mutation_deletion(chromo::Chromo)::Chromo
    deletion_point = rand(1:length(chromo))
    return deleteat!(copy(chromo), deletion_point)
end

function crossover(chromo1::Chromo, chromo2::Chromo)::Chromo
    cut1_idx = rand(1:length(chromo1))
    cut2_idx = rand(1:length(chromo2))
    new_chromo = chromo1[1:cut1_idx]
    for col_number in chromo2[cut2_idx:end]
        if !(col_number in new_chromo)
            push!(new_chromo, col_number)
        end
    end
    return new_chromo
end

function tournament_selection(
    scored_population::ScoredPopulation,
    penalties::Vector{Int},
)::Chromo
    best_chromo = nothing
    best_fitness = -Inf

    for _ = 1:TOURNAMENT_SIZE
        random_chromo, fitness = rand(scored_population)

        penalty = sum(map(col -> penalties[col], random_chromo))
        penalty /= length(random_chromo)
        penalty = OVERLAP_PENALTY^penalty

        if fitness - penalty >= best_fitness
            best_chromo = random_chromo
            best_fitness = fitness
        end
    end
    return best_chromo
end

function mutate(
    population::Population,
    old_scored_population::ScoredPopulation,
    tabu_list::Set,
    penalties::Vector{Int},
    cols_number::Int,
)::Tuple{Population,Int}
    tabu_hits = 0

    while length(population) < POPULATION_SIZE
        chromo1 = tournament_selection(old_scored_population, penalties)

        mutation = if rand() < RATE_CROSSOVER
            CROSSOVER
        else
            sample(collect(keys(ALL_MUTATIONS)), Weights(collect(values(ALL_MUTATIONS))))
        end

        new_chromo = if mutation == CROSSOVER
            chromo2 = tournament_selection(old_scored_population, penalties)
            crossover(chromo1, chromo2)
        elseif mutation == SWAP
            mutation_swap(chromo1)
        elseif mutation == SUBSTITUTION
            mutation_substitution(chromo1, cols_number)
        elseif mutation == INSERTION
            mutation_insertion(chromo1, cols_number)
        elseif mutation == DELETION
            mutation_deletion(chromo1)
        else
            error("Unsupported mutation $(mutation)")
        end

        @debug "$mutation: '$chromo1' -> '$new_chromo'"

        length(new_chromo) < MIN_CHROMO_SIZE && continue

        chromo_signature = hash(new_chromo)
        if !(chromo_signature in tabu_list)
            push!(population, new_chromo)
            push!(tabu_list, chromo_signature)
            for col in new_chromo
                penalties[col] += 1
            end
        else
            tabu_hits += 1
        end
    end

    return population, tabu_hits
end

end


text:
############################################################################################
#Include
include("stepnumber.jl")
include("kernel.jl")
include("config.jl")
include("diagnostics.jl")
include("postprocessing.jl")

############################################################################################
#export


text:
using RiskAdjustedLinearizations, JLD2, Test
include(joinpath(dirname(@__FILE__), "..", "..", "examples", "crw", "crw.jl"))

# Solve model
m_crw = CoeurdacierReyWinant()
m = crw(m_crw)
try
    solve!(m, m.z, m.y, m.Ψ; algorithm = :homotopy, step = .5, verbose = :none)
catch e
    local sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference/crw_sss.jld2"), "r")
    update!(ral, sssout["z_rss"], sssout["y_rss"], sssout["Psi_rss"])
end

# Calculate consumption at state zₜ
crw_cₜ(m, zₜ) = exp(m.y[1] + (m.Ψ * (zₜ - m.z))[1])

# Evaluates m_{t + 1} + r_{t + 1}
function crw_logSDFxR(m, zₜ, εₜ₊₁, Cₜ)
    zₜ₊₁, yₜ₊₁ = simulate(m, εₜ₊₁, zₜ)

    return log(m_crw.β) - m_crw.γ * (yₜ₊₁[1] - log(Cₜ)) + zₜ₊₁[2]
end

# Calculate 𝔼ₜ[exp(mₜ₊₁ + rₜ₊₁)] via quadrature
std_norm_mean = zeros(2)
std_norm_sig  = ones(2)
crw_𝔼_quadrature(f::Function) = gausshermite_expectation(f, std_norm_mean, std_norm_sig, 10)

# Calculate implied state variable(s)
function crw_endo_states(m, zₜ, zₜ₋₁, c_impl)
    # rₜ, yₜ are exogenous while Nₜ = exp(rₜ) * Aₜ₋₁ + Yₜ is entirely pre-determined.
    # Thus, our implied state variable will be foreign asset Aₜ = Nₜ - Cₜ.

    # zₜ₋₁ may be the previous period's implied state, so we start from there
    # to calculate Aₜ₋₁.
    yₜ₋₁ = m.y + m.Ψ * (zₜ₋₁ - m.z) # Calculate implied jump variables last period
    Cₜ₋₁ = exp(yₜ₋₁[1])             # to get the implied consumption last period.
    Aₜ₋₁ = zₜ₋₁[1] - Cₜ₋₁           # Given that consumption, we compute implied foreign assets yesterday.
    Nₜ   = exp(zₜ[2]) * Aₜ₋₁ + exp(zₜ[3]) # Now we can get implied resources available today.

    return vcat(zₜ, Nₜ - exp(c_impl)) # This gives us implied foreign assets today, along with other state variables
end

# Load draws from bivariate standard normal
shocks = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "crw_shocks.jld2"), "r")["shocks"]

@testset "Calculate Euler Equation Errors using Gauss-Hermite quadrature" begin

    # Calculate Euler Equation errors
    out1 = out2 = out3 = out4 = out5 = NaN
    for i in 1:100
        out1, out2, out3, out4, out5 = try
            abs.(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature; c_init = m.y[1] * 1.1)),
            abs.(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, m.z * 1.1; c_init = m.y[1] * 1.1)),
            abs.(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, m.z * 1.1;
                                      c_init = m.y[1] * 1.1, method = :newton)),
            abs(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, shocks,
                                     summary_statistic = x -> norm(x, Inf))),
            abs(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, shocks,
                                     summary_statistic = x -> norm(x, 2)))
        catch e
            NaN, NaN, NaN, NaN, NaN
        end

        if !isnan(out1)
            break
        end
        if i == 100
            out1, out2, out3, out4, out5 = abs.(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature;
                                                                     c_init = m.y[1] * 1.1)),
            abs.(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, m.z * 1.1; c_init = m.y[1] * 1.1)),
            abs.(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, m.z * 1.1;
                                      c_init = m.y[1] * 1.1, method = :newton)),
            abs(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, shocks,
                                     summary_statistic = x -> norm(x, Inf))),
            abs(euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, shocks,
                                     summary_statistic = x -> norm(x, 2)))
        end
    end

    @test out1 < 1e-10
    @test out2 < 1e-3
    @test out3 < 1e-3
    @test out4 < .6
    @test out5 < .6

    out6 = euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, shocks,
                                return_soln = true)
    _states, _ = simulate(m, size(shocks, 2), shocks, m.z)
    _impl_c = Vector{eltype(_states)}(undef, size(_states, 2))
    for i in 1:length(_impl_c)
        _impl_c[i] = crw_cₜ(m, (@view _states[:, i]))
    end
    @test abs.(maximum((out6 - _impl_c) ./ out6)) < .6

    c_ral, c_impl, endo_states_ral, endo_states_impl =
        dynamic_euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature, crw_endo_states, 1, shocks;
                                     return_soln = true)
    c_err, endo_states_err = dynamic_euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature,
                                                          crw_endo_states, 1, shocks; return_soln = false)
    @test_throws DimensionMismatch dynamic_euler_equation_error(m, crw_cₜ, crw_logSDFxR, crw_𝔼_quadrature,
                                                                crw_endo_states, 0, shocks; return_soln = false)
    @test c_err < 2e-5
    @test endo_states_err < 1e-3
    @test c_err == norm((c_ral - c_impl) ./ c_ral, Inf)
    @test endo_states_err == norm(vec(endo_states_ral - endo_states_impl) ./ vec(endo_states_ral), Inf)
end

nothing


text:
using Base.Iterators, PolygonOps

abstract type Edge end
struct EdgeX <: Edge
    range :: Tuple{Number, Number} # [a, b]
    y :: Function # y(x)
    a :: Number # ∂y
end

struct EdgeY <: Edge
    range :: Tuple{Number, Number} # [a, b]
    x :: Function # x(y)
    a :: Number # ∂x
end

function Edge(A, B)::Edge
    xRange, yRange = minmax.(A, B)
    (xa, ya), (xb, yb) = A, B
    if xa != xb
        a = (ya - yb) // (xa - xb)
        b = ya - a * xa
        EdgeX(xRange, x -> a * x + b, a)
    else
        a = (xa - xb) // (ya - yb)
        b = xa - a * ya
        EdgeY(yRange, y -> a * y + b, a)
    end
end

wrap1(xs) = flatten((rest(xs, 2), [xs[1]]))

onEdge((x, y), edge :: EdgeX) = (edge.range[1] <= x <= edge.range[2]) && y ≈ edge.y(x)
onEdge((x, y), edge :: EdgeY) = (edge.range[1] <= y <= edge.range[2]) && x ≈ edge.x(y)
onEdge(point, edges :: Vector{Edge}) = any([onEdge(point, edge) for edge in edges])

function prepareEdge(shape, Γ) :: Tuple{Any, Vector{Edge}}
    Γn, Γd = [], []
    for (edge, f) ∈ zip(Edge.(shape, wrap1(shape)), Γ)
        if f != 0 push!(Γn, (f, edge))
        else push!(Γd, edge) end
    end
    Γn, Γd
end

max_stars_repo_name: Pan-Maciek/FEM

max_issues_repo_name: Pan-Maciek/FEM

max_forks_repo_name: Pan-Maciek/FEM


text:
function define_indices_d(::Type{Gridded{Linear}}, d, pad)
    symix, symixp, symx = Symbol("ix_",d), Symbol("ixp_",d), Symbol("x_",d)
    quote
        $symix = clamp($symix, 1, size(itp, $d)-1)
        $symixp = $symix + 1
    end
end

function coefficients(::Type{Gridded{Linear}}, N, d)
    symix, symixp, symx = Symbol("ix_",d), Symbol("ixp_",d), Symbol("x_",d)
    sym, symp, symfx = Symbol("c_",d), Symbol("cp_",d), Symbol("fx_",d)
    symk, symkix = Symbol("k_",d), Symbol("kix_",d)
    quote
        $symkix = $symk[$symix]
        $symfx = ($symx - $symkix)/($symk[$symixp] - $symkix)
        $sym = 1 - $symfx
        $symp = $symfx
    end
end

function gradient_coefficients(::Type{Gridded{Linear}}, d)
    sym, symp = Symbol("c_",d), Symbol("cp_",d)
    symk, symix = Symbol("k_",d), Symbol("ix_",d)
    symixp = Symbol("ixp_",d)
    quote
        $symp = 1/($symk[$symixp] - $symk[$symix])
        $sym = - $symp
    end
end

# This assumes fractional values 0 <= fx_d <= 1, integral values ix_d and ixp_d (typically ixp_d = ix_d+1,
#except at boundaries), and an array itp.coefs
function index_gen{IT<:DimSpec{Gridded}}(::Type{Gridded{Linear}}, ::Type{IT}, N::Integer, offsets...)
    if length(offsets) < N
        d = length(offsets)+1
        sym = Symbol("c_", d)
        symp = Symbol("cp_", d)
        return :($sym * $(index_gen(IT, N, offsets..., 0)) + $symp * $(index_gen(IT, N, offsets..., 1)))
    else
        indices = [offsetsym(offsets[d], d) for d = 1:N]
        return :(itp.coefs[$(indices...)])
    end
end


text:
#=
tester:
- Julia version: 1.3.1
- Author: erathorn
- Date: 2020-10-07
OPENBLAS_NUM_THREADS=5 JULIA_NUM_THREADS=5 /home/jo/Julia14/julia-1.4.2/bin/julia -O3
=#
#using Distributed
#addprocs(2)
using Revise
using Pkg
Pkg.activate(".")
using MCPhylo
using Random
#using LinearAlgebra
#LinearAlgebra.BLAS.set_num_threads(5)

#@everywhere Random.seed!(42)
#
#@everywhere using LinearAlgebra
#@everywhere BLAS.set_num_threads(10)
#include("./src/MCPhylo.jl")
#using .MCPhylo
#using Random
Random.seed!(42)

#mt, df = make_tree_with_data("Example.nex", binary=true); # load your own nexus file
mt, df = make_tree_with_data("untracked_files/simulation_PNUTS_Paper/out_JC_20-600.nex", binary=false); # load your own nexus file

#mt, df = make_tree_with_data("untracked_files/simulation_PNUTS_Paper/out_Res_20-600.nex", binary=true); # load your own nexus file
#mt, df = make_tree_with_data("untracked_files/Dravidian.cc.phy.nex"); # load your own nexus file

mt2 = deepcopy(mt)
randomize!(mt2)
my_data = Dict{Symbol, Any}(
  :mtree => mt,
  :df => df,
  :nnodes => size(df)[3],
  :nbase => size(df)[1],
  :nsites => size(df)[2],
);

# model setup
model =  Model(
    df = Stochastic(3, (mtree, mypi) ->  PhyloDist(mtree, mypi, [1.0], [1.0], JC), false),
    mypi = Stochastic(1, () -> Dirichlet(4,1)),
    mtree = Stochastic(Node(), () -> TreeDistribution(CompoundDirichlet(1.0,1.0,0.100,1.0)), true)
     )
# intial model values
inits = [ Dict{Symbol, Union{Any, Real}}(
    :mtree => mt,
    :mypi=> rand(Dirichlet(4,1)),
    :df => my_data[:df],
    :nnodes => my_data[:nnodes],
    :nbase => my_data[:nbase],
    :nsites => my_data[:nsites],
    :a => rand(),
    ),
    Dict{Symbol, Union{Any, Real}}(
        :mtree => mt2,
        :mypi=> rand(Dirichlet(4,1)),
        :df => my_data[:df],
        :nnodes => my_data[:nnodes],
        :nbase => my_data[:nbase],
        :nsites => my_data[:nsites],
        :a => rand()
        )
    ]

scheme = [PNUTS(:mtree, tree_depth=5),#, target=0.6, targetNNI=4, tree_depth=5),
          #PPHMC(:mtree, 0.001, 10, 0.003),
          SliceSimplex(:mypi),
          ]

setsamplers!(model, scheme);

# do the mcmc simmulation. if trees=true the trees are stored and can later be
# flushed ot a file output.
sim = mcmc(model, my_data, inits, 100, burnin=50,thin=1, chains=2, trees=true)

# # request more runs
# sim = mcmc(sim, 1000)

# # write the output to a path specified as the second argument
# to_file(sim, "example_run")

# using BenchmarkTools
pd = PhyloDist(mt, [0.25,0.25,0.25,0.25], [1.0], [1.0], JC)

gradlogpdf(pd, df)
# @timev logpdf(pd, df)

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
#
# Date created: 2022-03-10
# Author: aradclif
#
#
############################################################################################

function aminusbbody(N::Int, D)
    body = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> Symbol(:i_, d), N)...)
    e = Expr(:(=), c, Expr(:call, :-, a, b))
    push!(body.args, e)
    body
end

function aminusb_quote(N::Int, D)
    ls = loopgen(N)
    body = aminusbbody(N, D)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return C
    end
end

@generated function aminusb!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, B::AbstractArray{T, N},
                             dims::D) where {Tₒ, T, N, D}
    aminusb_quote(N, D)
end

function exp!(A::AbstractArray{T, N}) where {T, N}
    @turbo for i ∈ eachindex(A)
        A[i] = exp(A[i])
    end
    A
end

function lvlogsoftmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = lvlogsoftmax1(A)
    else
        B = lvlse(A, dims=dims)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(B, d), N)
        C = similar(A, promote_type(T, eltype(B)))
        aminusb!(C, A, B, Dᴮ′)
    end
    return C
end
lvlogsoftmax(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvlogsoftmax(A, (dims,))
lvlogsoftmax(A::AbstractArray{T, N}; dims=:) where {T, N} = lvlogsoftmax(A, dims)
lvlogsoftmax(A::AbstractArray{T, N}) where {T, N} = lvlogsoftmax1(A)
lvlogsoftmax(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvlogsoftmax1(A)

function lvlogsoftmax1(A::AbstractArray{T, N}) where {T, N}
    b = lvlse1(A)
    C = similar(A, Base.promote_op(exp, T))
    @turbo for i ∈ eachindex(A)
        C[i] = A[i] - b
    end
    C
end

################ threaded version

function taminusb_quote(N::Int, D)
    ls = loopgen(N)
    body = aminusbbody(N, D)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return C
    end
end
@generated function taminusb!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    taminusb_quote(N, D)
end

function lvtlogsoftmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = lvtlogsoftmax1(A)
    else
        B = lvtlse(A, dims=dims)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(B, d), N)
        C = similar(A, promote_type(T, eltype(B)))
        taminusb!(C, A, B, Dᴮ′)
    end
    return C
end
lvtlogsoftmax(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvtlogsoftmax(A, (dims,))
lvtlogsoftmax(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtlogsoftmax(A, dims)
lvtlogsoftmax(A::AbstractArray{T, N}) where {T, N} = lvtlogsoftmax1(A)
lvtlogsoftmax(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvtlogsoftmax1(A)

function lvtlogsoftmax1(A::AbstractArray{T, N}) where {T, N}
    b = lvtlse1(A)
    C = similar(A, Base.promote_op(exp, T))
    @tturbo for i ∈ eachindex(A)
        C[i] = A[i] - b
    end
    C
end


text:
struct Rayleigh end

scattering(m, x, ::Rayleigh) = 8 * x^4 / 3 * abs((m^2 - 1) / (m^2 + 2))^2
absorption(m, x, ::Rayleigh) = 4 * x * imag((m^2 - 1) / (m^2 + 2))
extinction(m, x, s::Rayleigh) = scattering(m, x, s) + absorption(m, x, s)
backscattering(m, x, s::Rayleigh) = 3 * scattering(m, x, s) / 2

# s = Scatterer(1.33 + 0.01im, 25., 870.)
# scattering(s, Rayleigh())

@doc raw"""
    rayleigh_scattering(m, x)

Computes Mie efficencies of a spherical particle in the Rayleigh regime (``x = \pi \, d_p / \lambda \ll 1``) given refractive index `m`, `wavelength`, and `diameter`. Uses Rayleigh-regime approximations:
```math
Q_{sca} = \frac{8x^4}{3} \left|{\frac{m^2-1}{m^2+2}}\right|^2
```
```math
Q_{abs} = 4x \: \text{Im}\left\{\frac{m^2-1}{m^2+2}\right\}
```
```math
Q_{ext} = Q_{sca}+ Q_{abs}
```
```math
Q_{back} = \frac{3Q_{sca}}{2}
```
```math
Q_{ratio} = 1.5
```
```math
Q_{pr} = Q_{ext}
```

# Arguments
- `m`: The complex refractive index, with the convention ``m = n+ik``.
- `wavelength`: The wavelength of incident light, in nanometers.
- `diameter`: The diameter of the particle, in nanometers.
- `n_medium`: The refractive index of the surrounding medium. This must be positive, nonzero, and real. Any imaginary part will be discarded.

# Returns
- `qext`, `qsca`, `qabs`, `qback`: The Mie efficencies described above.

# Examples
```julia-repl
julia> rayleigh_scattering(mx(1.33 + 0.01im, 25., 870.)...)

(0.0041753430994240295, 0.00011805645915412197, 0.004057286640269908, 0.00017708468873118297)
```
# Notes
-  different formula on [wiki](https://en.m.wikipedia.org/wiki/Rayleigh_scattering)
"""
function rayleigh_scattering(m::ComplexF64, x::Float64)

    if x ≈ 0
        return 0, 0, 0, 0
    end

    LL = (m^2 - 1) / (m^2 + 2) # Lorentz-Lorenz term

    qsca = 8 * x^4 / 3 * abs(LL)^2  # B&H eq 5.8
    qabs = 4 * x * imag(LL) # B&H eq. 5.11
    qext = qsca + qabs
    qback = 3 * qsca / 2 # B&H eq. 5.9

    # qratio = 1.5
    # qpr = qext

    return qext, qsca, qabs, qback
end

# rayleigh_scattering(mx(1.33 + 0.01im, 0., 870.)...)
# rayleigh_scattering(mx(1.33 + 0.01im, 25., 870.)...)
# #
# lambda = LinRange(600, 1000, 1001)
# qext =  [rayleigh_scattering(mx(1.33 + 0.01im, 25., l)...)[1] for l in lambda]
# plot(lambda, qext)
#

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


text:
tu = (1, 2, "Hello")
# 1 based
println(tu[1])
println(tu[3])

named = (first = 100, second = 10)
println(named[1])
println(named.first)

map(x -> x * 10, named) |> println

(a, b, c) = 2:4
println(a, b, c)

println(length(tu))
println(lastindex(tu))

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenBLAS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenBLAS")
JLLWrappers.@generate_main_file("OpenBLAS", UUID("4536629a-c528-5b80-bd46-f80d51c5b363"))
end  # module OpenBLAS_jll


text:
using Optim, Random, Distributions, Suppressor, Statistics, LinearAlgebra
using FixedEffectModels, DataFrames
includet("estimator_binary.jl")

struct UnobservedType
        v::Array{Float64}
        c::Array{Float64}
end

function UnobservedType(n)
        v = rand(Uniform(7, 12), n)
        c = rand(Uniform(5, 10), n)
        return UnobservedType(v, c)
end

function good_demand(p, v)
        return (p.<v).*1.0
end

function good_supply(p, c, tech, subsidy::Bool)
        if subsidy
            return (p .- tech .> c)
        else
            return tech.*(p .> c)
        end

end

function eqm_price(U, types::UnobservedType, tech, subsidy)
        function eqm_condition(p)
            zp = mean(good_demand(p .+ U, types.v)) - mean(good_supply(p .+ U, types.c, tech, subsidy))
            return abs(zp)
        end
        r = optimize(eqm_condition, 0.1, 50.0).minimizer
        return r
end

#technology that increases production by 20% at the same cost
function simulate_eqm_data(types, n; subsidy= false, p = 0.5)
    #types = UnobservedType(n)
    W = rand(Bernoulli(p), n)
    #W = append!(ones(5000), zeros(5000))
    U = 0.18 .* rand([-1,1], n)
    #binary subsidy
    if subsidy
        tech = W.*0.5 + (1 .- W).*0.0
    #production technology intervention
    else
        tech = W*1.2 + (1 .- W).*1.0
    end
    pstar = eqm_price(U, types, tech, subsidy)
    #println("Ground truth direct: ", direct_ground(types, pstar))
    D = good_demand(pstar .+ U, types.v)
    S = good_supply(pstar .+ U, types.c, tech, subsidy)
    Y = S.*(pstar .+ tech.*subsidy .+ U .- types.c)
    return EqmData(Y, D, S, W, U)
end

function run_simulation(n, S)
        results = zeros((S, 4))
        for i in 1:S
                types = UnobservedType(n)
                data = simulate_eqm_data(types, n)
                estimates = estimate_τ(data)
                results[i, 1] = estimates.τD
                results[i, 2] = estimates.τI
                results[i, 3] = estimates.σD
                results[i, 4] = estimates.σI
        end

        gt = [0.22222222222, 0.7333333333333*35/10.5^2*(-1.0)]
        for j in 1:2
            println("Effect Size: ", mean(results[:, j]))
            println("Bias: ", mean(results[:, j].- gt[j]))
            println("SD: ", std(results[:, j]))
            ciu = results[:, j] .+ results[:, j+2]*1.96
            cil = results[:, j] .- results[:, j+2]*1.96
            covg = (gt[j] .< ciu) .& (gt[j] .> cil)
            println("Coverage: ", mean(covg))
            println("CI Width: ", mean(ciu .- cil))
        end
        return results
end

max_stars_repo_name: evanmunro/market-equilibrium

max_issues_repo_name: evanmunro/market-equilibrium

max_forks_repo_name: evanmunro/market-equilibrium


text:
# syntax: proto2
# COV_EXCL_START
using ProtoBuf
import ProtoBuf.meta

# service methods for ContainerManagementProtocolService
const _ContainerManagementProtocolService_methods = MethodDescriptor[
        MethodDescriptor("startContainers", 1, hadoop.yarn.StartContainersRequestProto, hadoop.yarn.StartContainersResponseProto),
        MethodDescriptor("stopContainers", 2, hadoop.yarn.StopContainersRequestProto, hadoop.yarn.StopContainersResponseProto),
        MethodDescriptor("getContainerStatuses", 3, hadoop.yarn.GetContainerStatusesRequestProto, hadoop.yarn.GetContainerStatusesResponseProto),
        MethodDescriptor("increaseContainersResource", 4, hadoop.yarn.IncreaseContainersResourceRequestProto, hadoop.yarn.IncreaseContainersResourceResponseProto),
        MethodDescriptor("updateContainer", 5, hadoop.yarn.ContainerUpdateRequestProto, hadoop.yarn.ContainerUpdateResponseProto),
        MethodDescriptor("signalToContainer", 6, hadoop.yarn.SignalContainerRequestProto, hadoop.yarn.SignalContainerResponseProto),
        MethodDescriptor("localize", 7, hadoop.yarn.ResourceLocalizationRequestProto, hadoop.yarn.ResourceLocalizationResponseProto),
        MethodDescriptor("reInitializeContainer", 8, hadoop.yarn.ReInitializeContainerRequestProto, hadoop.yarn.ReInitializeContainerResponseProto),
        MethodDescriptor("restartContainer", 9, hadoop.yarn.ContainerIdProto, hadoop.yarn.RestartContainerResponseProto),
        MethodDescriptor("rollbackLastReInitialization", 10, hadoop.yarn.ContainerIdProto, hadoop.yarn.RollbackResponseProto),
        MethodDescriptor("commitLastReInitialization", 11, hadoop.yarn.ContainerIdProto, hadoop.yarn.CommitResponseProto)
    ] # const _ContainerManagementProtocolService_methods
const _ContainerManagementProtocolService_desc = ServiceDescriptor("hadoop.yarn.ContainerManagementProtocolService", 1, _ContainerManagementProtocolService_methods)

ContainerManagementProtocolService(impl::Module) = ProtoService(_ContainerManagementProtocolService_desc, impl)

mutable struct ContainerManagementProtocolServiceStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    ContainerManagementProtocolServiceStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_ContainerManagementProtocolService_desc, channel))
end # mutable struct ContainerManagementProtocolServiceStub

mutable struct ContainerManagementProtocolServiceBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    ContainerManagementProtocolServiceBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_ContainerManagementProtocolService_desc, channel))
end # mutable struct ContainerManagementProtocolServiceBlockingStub

startContainers(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.StartContainersRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[1], controller, inp, done)
startContainers(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.StartContainersRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[1], controller, inp)

stopContainers(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.StopContainersRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[2], controller, inp, done)
stopContainers(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.StopContainersRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[2], controller, inp)

getContainerStatuses(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainerStatusesRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[3], controller, inp, done)
getContainerStatuses(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.GetContainerStatusesRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[3], controller, inp)

increaseContainersResource(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.IncreaseContainersResourceRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[4], controller, inp, done)
increaseContainersResource(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.IncreaseContainersResourceRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[4], controller, inp)

updateContainer(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerUpdateRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[5], controller, inp, done)
updateContainer(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerUpdateRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[5], controller, inp)

signalToContainer(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.SignalContainerRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[6], controller, inp, done)
signalToContainer(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.SignalContainerRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[6], controller, inp)

localize(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ResourceLocalizationRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[7], controller, inp, done)
localize(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ResourceLocalizationRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[7], controller, inp)

reInitializeContainer(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ReInitializeContainerRequestProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[8], controller, inp, done)
reInitializeContainer(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ReInitializeContainerRequestProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[8], controller, inp)

restartContainer(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerIdProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[9], controller, inp, done)
restartContainer(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerIdProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[9], controller, inp)

rollbackLastReInitialization(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerIdProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[10], controller, inp, done)
rollbackLastReInitialization(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerIdProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[10], controller, inp)

commitLastReInitialization(stub::ContainerManagementProtocolServiceStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerIdProto, done::Function) = call_method(stub.impl, _ContainerManagementProtocolService_methods[11], controller, inp, done)
commitLastReInitialization(stub::ContainerManagementProtocolServiceBlockingStub, controller::ProtoRpcController, inp::hadoop.yarn.ContainerIdProto) = call_method(stub.impl, _ContainerManagementProtocolService_methods[11], controller, inp)

export ContainerManagementProtocolService, ContainerManagementProtocolServiceStub, ContainerManagementProtocolServiceBlockingStub, startContainers, stopContainers, getContainerStatuses, increaseContainersResource, updateContainer, signalToContainer, localize, reInitializeContainer, restartContainer, rollbackLastReInitialization, commitLastReInitialization
# COV_EXCL_STOP


text:
#########
#File: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\scripts\simulationsAndPlotsPaper\reciprocity_p_star\simulate_AR_and_estimate_variance_edges_GWESP.jl
#Project: c:\Users\digan\Dropbox\Dynamic_Networks\repos\ScoreDrivenExponentialRandomGraphs\scripts\simulationsAndPlotsPaper\reciprocity_p_star
#Created Date: Wednesday April 28th 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Sunday May 30th 2021 12:04:11 pm
#Modified By:  Domenico Di Gangi
#-----
#Description: test sequence of simulations in the current folder for a single realization
#-----
########

using DrWatson
@quickactivate "ScoreDrivenExponentialRandomGraphs"

using ProjUtilities
using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets, DynNets, Utilities

model = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges +  gwesp(decay = 0.25,fixed = TRUE)", true), indTvPar = [true, true], scoreScalingType="FISH_D")

N = 100
T = 110

import ScoreDrivenERGM.DynNets: reference_model, seq_of_obs_from_seq_of_mats, estimate_and_filter
import ScoreDrivenERGM.Utilities.AReg

dgpSettings = (type = "AR", opt = (B =[0.95], sigma = [0, 0.2],  minVals = [-3.0, 0], maxVals = [-0.1, 2] ))

parDgpT = DynNets.sample_time_var_par_from_dgp(reference_model(model), dgpSettings.type, N, T;  dgpSettings.opt..., maxAttempts=100000)

@elapsed A_T_dgp = StaticNets.sample_ergm_sequence(reference_model(model).staticModel, N, parDgpT, 1)[:,:,:,1]

obsT = seq_of_obs_from_seq_of_mats(model, A_T_dgp)

filtModel = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges +  gwesp(decay = 0.25,fixed = TRUE)", true), indTvPar = [false, true], scoreScalingType="FISH_D")

@elapsed ~, vEstSdResPar, fVecT_filt, ~, ~, conv_flag, ftot_0 = estimate_and_filter(filtModel, N, obsT;  show_trace = true)

using PyPlot

plot(fVecT_filt[2,:])
plot(parDgpT[2,:])

mapslices( x-> AReg.fitARp(Float64.(x), 1)[2], fVecT_filt, dims=2 )

fVecT_filt
#endregion

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
@testset "testing size" begin
for i in 1:10
    testcube = rand(1:10, 10, 10, 10)
    test_sparsecube = SatelliteImages.sparse_cube(testcube)
    @test size(test_sparsecube) == (10, 10, 10)
end
end


text:
function create_grids()

    function xgrid(theta,xval)
        N      = length(xval)
        xub    = vcat(xval[2:N],Inf)
        xtran1  = zeros(N,N)
        xtran1c = zeros(N,N)
        lcdf   = zeros(N)
        for i=1:length(xval)
            xtran1[:,i]   = (xub[i] .>= xval).*(1 .- exp.(-theta*(xub[i] .- xval)) .- lcdf)
            lcdf        .+= xtran1[:,i]
            xtran1c[:,i] .+= lcdf
        end
        return xtran1, xtran1c
    end

    zval   = collect([.25:.01:1.25]...)
    zbin   = length(zval)
    xval   = collect([0:.125:25]...)
    xbin   = length(xval)
    tbin   = xbin*zbin
    xtran  = zeros(tbin,xbin)
    xtranc = zeros(xbin,xbin,xbin)
    for z=1:zbin
        xtran[(z-1)*xbin+1:z*xbin,:],xtranc[:,:,z] = xgrid(zval[z],xval)
    end

    return zval,zbin,xval,xbin,xtran
end


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

mutable struct HopLinkProperties <: SwaggerModel
    roundTripTimeMin::Any # spec type: Union{ Nothing, Int32 } # spec name: roundTripTimeMin
    roundTripTimeAvg::Any # spec type: Union{ Nothing, Int32 } # spec name: roundTripTimeAvg
    roundTripTimeMax::Any # spec type: Union{ Nothing, Int32 } # spec name: roundTripTimeMax

    function HopLinkProperties(;roundTripTimeMin=nothing, roundTripTimeAvg=nothing, roundTripTimeMax=nothing)
        o = new()
        validate_property(HopLinkProperties, Symbol("roundTripTimeMin"), roundTripTimeMin)
        setfield!(o, Symbol("roundTripTimeMin"), roundTripTimeMin)
        validate_property(HopLinkProperties, Symbol("roundTripTimeAvg"), roundTripTimeAvg)
        setfield!(o, Symbol("roundTripTimeAvg"), roundTripTimeAvg)
        validate_property(HopLinkProperties, Symbol("roundTripTimeMax"), roundTripTimeMax)
        setfield!(o, Symbol("roundTripTimeMax"), roundTripTimeMax)
        o
    end
end # type HopLinkProperties

const _property_map_HopLinkProperties = Dict{Symbol,Symbol}(Symbol("roundTripTimeMin")=>Symbol("roundTripTimeMin"), Symbol("roundTripTimeAvg")=>Symbol("roundTripTimeAvg"), Symbol("roundTripTimeMax")=>Symbol("roundTripTimeMax"))
const _property_types_HopLinkProperties = Dict{Symbol,String}(Symbol("roundTripTimeMin")=>"Int32", Symbol("roundTripTimeAvg")=>"Int32", Symbol("roundTripTimeMax")=>"Int32")
Base.propertynames(::Type{ HopLinkProperties }) = collect(keys(_property_map_HopLinkProperties))
Swagger.property_type(::Type{ HopLinkProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_HopLinkProperties[name]))}
Swagger.field_name(::Type{ HopLinkProperties }, property_name::Symbol) =  _property_map_HopLinkProperties[property_name]

function check_required(o::HopLinkProperties)
    true
end

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


text:
@testset "test `dfgadflysegplot.jl`" begin
	using DataFrames, Gadfly
	x = randn(50);
	y = rand(1:9999,50);
	df = DataFrame(:col1 => x, :col2 => y);
	(df2, colind) = forsegplot(df, [:col1, :col2]);
	@test isequal(df2[!,:x0], x[1:end-1])
	@test isequal(df2[!,:x1], x[2:end])
	@test isequal(df2[!,:y0], y[1:end-1])
	@test isequal(df2[!,:y1], y[2:end])

	plotwithouterror = true;
	try
		Gadfly.plot(df2, Geom.segment; colind...);
	catch
		plotwithouterror = false;
	end
	@test plotwithouterror

end


text:
struct FileStorage <: DataStorage
    path::Any
    verbose::Any
    force::Any
    preferences::Any
end

FileStorage(path; verbose = true, force = false, kwargs...) =
    FileStorage(path, verbose, force, Dict(kwargs...))

gethandle(storage::FileStorage) = joinpath(storage.path, "data")

function savemetadata(storage::FileStorage, metadata)
    file = string(joinpath(storage.path, "meta.jld2"))
    JLD2.save(file, metadata)
end

function loadmetadata(storage::FileStorage)
    JLD2.load(string(joinpath(storage.path, "meta.jld2")))
end

function savedataset!(path, storage::FileStorage, blocks::Tuple, datas::Tuple)
    ensuredir(storage, path)
    for (i, (block, data)) in enumerate(zip(blocks, datas))
        subpath = joinpath(path, "block_$(lpad(i, 4, '0'))")
        savedataset!(subpath, storage, block, data)
    end
end

function savedataset!(path, storage::FileStorage, blocks::Tuple, datas)
    ensuredir(storage, path)
    for (i, block) in enumerate(blocks)
        data = mapobs(obs -> obs[i], datas)
        subpath = joinpath(path, "block_$(lpad(i, 4, '0'))")
        savedataset!(subpath, storage, block, data)
    end
end

function ensuredir(storage::FileStorage, path)
    if isdir(path) || isfile(path)
        if storage.force
            rm(path, recursive = true)
        else
            error("Directory $path already exists, aborting.")
        end
    end
    mkdir(path)
end

function loaddataset(path, storage::FileStorage, blocks::Tuple)
    datas = []
    for (i, block) in enumerate(blocks)
        subpath = joinpath(path, "block_$(lpad(i, 4, '0'))")
        data = loaddataset(subpath, storage, block)
        push!(datas, data)
    end
    return Tuple(datas)
end


text:
function load(path::String, sample::Union{EzXML.Node,Nothing}; transform::Function=x->asinh(x/250), kwargs...)

	#######################################
    params, data = loadFCS(path; kwargs...)
	params = sample === nothing ? getMetaData(params) : findchannels(sample)

	###################################### data with human readable channel names
	channels = dropmissing(params,"E")
	transform!(channels, ["N","S"] => ByRow((N,S)-> ~ismissing(S) ? S != "" ? S : N : N) => "name")
	data = DataFrame(data,channels.name)

	###################################### biexponential transformation
	@. data = transform(data)

    ################################ load metadata from workspace
	groups = loadGroups(data,sample)
	gating = gatingGraph(sample; transform=transform,
		channelMap=Dict( param.N => ismissing(param.S) ? param.N : param.S for param ∈ eachrow(params) )
	)

	labels = gate(data,gating)
	return data,labels,groups,gating
end

function load(path::String; files::String=joinpath(dirname(path),"*.fcs"), transform::Function=x->asinh(x/250), cols::Symbol=:setequal, kwargs...)

	@assert( isfile(path), "no such file: $path")
	@assert( length(glob(files)) ≠ 0, "no FCS files found using pattern: $files")

	try
		workspace = root(readxml(path))
	catch error
		throw(LoadError(path,error.line,error))
	end

	workspace = root(readxml(path))
	datasets = map( dataset -> basename(dataset["uri"]), findall("//DataSet",workspace) )
	@assert( length(datasets) == length(unique(datasets)), "FCS files under a workspace must have unique names" )

	data,labels,groups = DataFrame(),DataFrame(),DataFrame()
	gatings = Dict()

	for path ∈ glob(files)
		@info "Loading $path"
		sample = findsample(path,workspace)

		fcs,label,group,gating = load(path,sample;transform=transform,kwargs...)

		append!(data,fcs,cols=cols)
		append!(labels,label,cols=:union)
		append!(groups,group,cols=:union)
		gatings[path] = gating
	end

	map( name->replace!(labels[!,name],missing=>false), names(labels) )
	map( name->replace!(groups[!,name],missing=>false), names(groups) )

	disallowmissing!(labels)
	disallowmissing!(groups)

	transform!( labels, AsTable(filter(x->x≠"Unlabelled",names(labels))) => ByRow(x->~any(x)) => "Unlabelled" )
	transform!( groups, AsTable(filter(x->x≠"Ungrouped",names(groups))) =>  ByRow(x->~any(x)) => "Ungrouped" )

	return data,labels,groups,gatings
end

function findsample(uri::String,workspace::EzXML.Node)
	uri = escapeuri(uri)

	for code ∈ [("%2F"=>"/"), ("%5C"=>"\\"), ("%3A"=>":"), ("%2C"=>",")]
		uri = replace(uri, code)
	end

	sample = findfirst("//DataSet[contains(@uri,'$(basename(uri))')]",workspace)
	sample === nothing && @warn("""no metadata found""")
	return sample
end

function findchannels(sample::EzXML.Node)

	params = DataFrame( "keyword" => findall("..//Keywords/Keyword[starts-with(@name,'\$P') and contains(@name,'N')]", sample))
	keys = ["G","R","V","S","B","N","E","AR"]

	for key ∈ keys
		transform!(params, "keyword" => ByRow( x->getparam(chop(x["name"],head=0,tail=1)*key, sample)) => key)
	end

	transform!(params, "keyword" => ByRow( x->parse(Int,chop(x["name"],head=2,tail=1))) => "index")
	sort!(params,"index")

	transform!(params,["N","S"]=>ByRow((N,S)-> ~ismissing(S) ? S : stripname(N,params.N))=>"S")
	return select!(params,keys)
end

function stripname(name::String,names::AbstractArray{<:String})
    names = filter(x->occursin(x,name)&(x≠name), names)
    return isempty(names) ? missing : first(names)
end

function getparam(param::String,sample::EzXML.Node)
	keyword = findfirst("""..//Keywords/Keyword[@name='$param']""", sample)
	return keyword !== nothing ? keyword["value"] != "" ? keyword["value"] : missing : missing
end


text:
"""
# SequentialSamplingModels.jl

A collection of sequential sampling models based on the Distributions.jl API.

## Currently Supported models

- `LBA`: Linear Ballistic Accumulator
- `LNR`: Lognormal Race Model
- `Wald`: a shifted Wald represented a single boundary diffusion process
- `WaldMixture`: a shifted Wald represented a single boundary diffusion process with across-trial
    variability in the drift rate
- `AttentionalDiffusion`: a drift diffusion model in which the accumulation process is determined by the
utility of a visually attended option
"""
module SequentialSamplingModels
    using Distributions, Parameters, ConcreteStructs, PrettyTables
    import Distributions: pdf, logpdf, rand, loglikelihood, mean, std, cdf
    import Distributions: logccdf
    export SequentialSamplingModel, Wald, WaldMixture, LNR, LBA, DiffusionRace, AttentionalDiffusion
    export pdf, cdf, logpdf, rand, loglikelihood, mean, std

    include("utilities.jl")
    include("LogNormalRace.jl")
    include("Wald.jl")
    include("LBA.jl")
    include("DiffusionRace.jl")
    include("AttentionalDiffusion.jl")
end


text:
using FormulationLattice
using Base.Test

let
    @Literals(A, B, C, D)
    cl = A ∨ ((B ∨ C) ∧ D)
    formtrack = FormulaState[]
    cl2 = dnf(cl, formtrack)
    @test cl2 == A ∨ (B ∧ D) ∨ (C ∧ D)
end

let
    @Literals(A, B, C, D)
    cl = (A ∨ B) ∧ (C ∨ D)
    formtrack = FormulaState[]
    cl2 = dnf(cl, formtrack)
    @test cl2 == (A ∧ C) ∨ (A ∧ D) ∨ (B ∧ C) ∨ (B ∧ D)
end

let
    @Literals(A, B, C, D)
    cl = (A ∨ B) ∧ C ∧ D
    formtrack = FormulaState[]
    cl2 = dnf(cl, formtrack)
    @test cl2 == (A ∧ C ∧ D) ∨ (B ∧ C ∧ D)
end

let
    @Literals(A, B, C, D, E, F)
    cl = (A ∨ B) ∧ (C ∨ D) ∧ (E ∨ F)
    formtrack = FormulaState[]
    cl2 = dnf(cl, formtrack)
    @test cl2 == (A ∧ C ∧ E) ∨ (A ∧ C ∧ F) ∨ (A ∧ D ∧ E) ∨ (A ∧ D ∧ F) ∨
                 (B ∧ C ∧ E) ∨ (B ∧ C ∧ F) ∨ (B ∧ D ∧ E) ∨ (B ∧ D ∧ F)
end

let
    @Literals(A12, A13)
    @Literals(B12, B13)
    @Literals(C12, C13)
    cl = (A12 ∨ B12 ∨ C12) ∧ (A13 ∨ B13 ∨ C13)
    formtrack = FormulaState[]
    cl2 = dnf(cl, formtrack)
    @test cl2 == (A12 ∧ A13) ∨ (A12 ∧ B13) ∨ (A12 ∧ C13) ∨
                 (B12 ∧ A13) ∨ (B12 ∧ B13) ∨ (B12 ∧ C13) ∨
                 (C12 ∧ A13) ∨ (C12 ∧ B13) ∨ (C12 ∧ C13)
end

let
    @Literals(A12, A13, A23)
    @Literals(B12, B13, B23)
    @Literals(C12, C13, C23)
    cl = (A12 ∨ B12 ∨ C12) ∧ (A13 ∨ B13 ∨ C13) ∧ (A23 ∨ B23 ∨ C23)
    @Literals(A12, A13)
    @Literals(B12, B13)
    @Literals(C12, C13)
    cl = (A12 ∨ B12 ∨ C12) ∧ (A13 ∨ B13 ∨ C13)
    formtrack = FormulaState[]
    cl2 = dnf(cl, formtrack)
    @test cl2 == (A12 ∧ A13 ∧ A23) ∨ (A12 ∧ A13 ∧ B23) ∨ (A12 ∧ A13 ∧ C23) ∨ (A12 ∧ B13 ∧ A23) ∨ (A12 ∧ B13 ∧ B23) ∨ (A12 ∧ B13 ∧ C23) ∨ (A12 ∧ C13 ∧ A23) ∨ (A12 ∧ C13 ∧ B23) ∨ (A12 ∧ C13 ∧ C23) ∨
                 (B12 ∧ A13 ∧ A23) ∨ (B12 ∧ A13 ∧ B23) ∨ (B12 ∧ A13 ∧ C23) ∨ (B12 ∧ B13 ∧ A23) ∨ (B12 ∧ B13 ∧ B23) ∨ (B12 ∧ B13 ∧ C23) ∨ (B12 ∧ C13 ∧ A23) ∨ (B12 ∧ C13 ∧ B23) ∨ (B12 ∧ C13 ∧ C23) ∨
                 (C12 ∧ A13 ∧ A23) ∨ (C12 ∧ A13 ∧ B23) ∨ (C12 ∧ A13 ∧ C23) ∨ (C12 ∧ B13 ∧ A23) ∨ (C12 ∧ B13 ∧ B23) ∨ (C12 ∧ B13 ∧ C23) ∨ (C12 ∧ C13 ∧ A23) ∨ (C12 ∧ C13 ∧ B23) ∨ (C12 ∧ C13 ∧ C23)
end

let
    @Literals(A12, A13, A14, A23, A24, A34)
    @Literals(B12, B13, B14, B23, B24, B34)
    @Literals(C12, C13, C14, C23, C24, C34)
    @Literals(D12, D13, D14, D23, D24, D34)
    cl = (A12 ∨ B12 ∨ C12 ∨ D12) ∧ (A13 ∨ B13 ∨ C13 ∨ D13) ∧
         (A14 ∨ B14 ∨ C14 ∨ D14) ∧ (A23 ∨ B23 ∨ C23 ∨ D23) ∧
         (A24 ∨ B24 ∨ C24 ∨ D24) ∧ (A34 ∨ B34 ∨ C34 ∨ D34)
    formtrack = FormulaState[]
    cl2 = dnf(cl, formtrack)
    @test length(cl2.clauses) == 4096
    @test all(cl2.clauses) do x
        length(x.clauses) == 6 &&
        isa(x,And) &&
        all(y -> isa(y,Literal), x.clauses)
    end
end


text:
using Documenter, GraphicalModelLearning

makedocs(
    modules = [GraphicalModelLearning],
    format = Documenter.HTML(),
    sitename = "GraphicalModelLearning",
    authors = "Carleton Coffrin, Andrey Y. Lokhov, Sidhant Misra, Marc Vuffray, and contributors.",
    pages = [
        "Home" => "index.md",
        "Library" => "library.md"
    ]
)

deploydocs(
    repo = "github.com/lanl-ansi/GraphicalModelLearning.jl.git",
)


text:
using LULESH
using MPI
using Enzyme
using LLVM

LLVM.clopts("-memdep-block-scan-limit=70000")
LLVM.clopts("-dse-memoryssa-walklimit=10000")
LLVM.clopts("-attributor-max-iterations=128")
LLVM.clopts("-capture-tracking-max-uses-to-explore=256")

# Enzyme.API.printperf!(true)
# Enzyme.API.printall!(true)
# Enzyme.API.instname!(true)

Enzyme.API.inlineall!(true)
# Size of Domain is around 1024
Enzyme.API.maxtypeoffset!(32)
isdefined(Enzyme.API, :strictAliasing!) && Enzyme.API.strictAliasing!(true)
isdefined(Enzyme.API, :typeWarning!) &&  Enzyme.API.typeWarning!(false)
Enzyme.API.looseTypeAnalysis!(true)

function main(nx, structured, num_iters, mpi, enzyme)
    # TODO: change default nr to 11
    nr = 1
    balance = 1
    cost = 1
    floattype = Float64
    devicetype = Vector

    if mpi
        !MPI.Initialized() && MPI.Init()
        comm = MPI.COMM_WORLD
    else
        comm = nothing
    end

    if getMyRank(comm) == 0
        @info "Constructing LuleshProblem" num_iters structured nx nr balance cost ranks=getNumRanks(comm)
    end

    prob = LuleshProblem(num_iters, structured, nx, nr, balance, cost, devicetype, floattype, comm)

    # Set up the mesh and decompose. Assumes regular cubes for now
    # TODO: modify this constructor to account for new fields
    # TODO: setup communication buffers

    domain = Domain(prob)
    if enzyme
        shadowDomain = Domain(prob)
        @time Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
        domain = Domain(prob)
    end

    # getnodalMass = nodalMass(domain)

    # Initial domain boundary communication
    commRecv(domain, MSG_COMM_SBN, 1,
                domain.sizeX + 1, domain.sizeY + 1, domain.sizeZ + 1,
                true, false)
    fields = (domain.nodalMass,)
    commSend(domain, MSG_COMM_SBN, fields,
             domain.sizeX + 1, domain.sizeY + 1, domain.sizeZ + 1,
             true, false)
    commSBN(domain, fields)

    # End initialization
    if mpi
        MPI.Barrier(prob.comm::MPI.Comm)
    end

    if getMyRank(prob.comm) == 0
        if (structured)
            @info "Running" until=domain.stoptime domain=(nx,nx,nx)
        else
            @info "Running" until=domain.stoptime domain=domain.numElem
            @warn "Unstructured setup not supported"
        end
    end

    # timestep to solution
    start = getWtime(prob.comm)
    while domain.time < domain.stoptime
        # this has been moved after computation of volume forces to hide launch latencies
        timeIncrement!(domain)
        if enzyme
            Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
        else
            lagrangeLeapFrog(domain)
        end

        # checkErrors(domain, its, myRank)
        if getMyRank(prob.comm) == 0
            @info "Completed" cycle=domain.cycle time=domain.time dt=domain.deltatime
        end
        if domain.cycle >= num_iters
            break
        end
    end

    # Use reduced max elapsed time
    elapsed_time = getWtime(prob.comm) - start
    elapsed_timeG = comm_max(elapsed_time, prob.comm)

    if getMyRank(prob.comm) == 0
        verifyAndWriteFinalOutput(elapsed_timeG, domain, nx, getNumRanks(prob.comm))
    end

    if mpi
        MPI.Finalize()
    end
end

if !isinteractive()
    args = LULESH.parse_cmd()

    # assume cube subdomain geometry for now (nx)
    nx = args["N"]
    structured = args["s"]
    num_iters = args["num_iters"]
    mpi = args["mpi"]
    enzyme = args["enzyme"]
    main(nx, structured, num_iters, mpi, enzyme)
end


text:
# a basic create robot type node example

using Caesar, IncrementalInference
using CloudGraphs

using KernelDensityEstimate

# switch IncrementalInference to use CloudGraphs (Neo4j) data layer
# dbaddress = "localhost"
# dbusr = ""
# dbpwd = ""
# mongoaddress = "localhost"
session = "SESSTEST"
include(joinpath(dirname(@__FILE__) ,"blandauth.jl"))
Nparticles = 200
println("Attempting to solve session $(session)...")

configuration = CloudGraphs.CloudGraphConfiguration(dbaddress, 7474, dbusr, dbpwd, mongoaddress, 27017, false, "", "");
cloudGraph = connect(configuration, encodePackedType, getpackedtype, decodePackedType);
# register types of interest in CloudGraphs
registerGeneralVariableTypes!(cloudGraph)
Caesar.usecloudgraphsdatalayer!()

# Uncomment out for command line operation
# cloudGraph, addrdict = standardcloudgraphsetup(drawdepth=true)
# session = addrdict["session"]
# Nparticles = addrdict["num particles"]

# this is being replaced by cloudGraph, added here for development period
fg = Caesar.initfg(sessionname=session, cloudgraph=cloudGraph)

# Robot navigation and inference type stuff
N=Nparticles
doors = [-100.0;0.0;100.0;300.0]'
cov = [3.0]

# robot style, add first pose vertex
v1 = addNode!(fg,:x1,doors,N=N,labels=["POSE"])

# add a prior for the initial position of the robot
f0  = addFactor!(fg,[v1], Obsv2(doors, cov', [1.0]))

# add second pose vertex
tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg, :x2, tem, N=N,labels=["POSE"])

# now add the odometry factor between them
f1 = addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

v3=addNode!(fg,:x3,4.0*randn(1,N)+getVal(v2)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))
f2 = addFactor!(fg,[v3], Obsv2(doors, cov', [1.0]))

v4=addNode!(fg,:x4,2.0*randn(1,N)+getVal(v3)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

# if true
l1=addNode!(fg, :l1, 0.5*randn(1,N)+getVal(v3)+64.0, N=N,labels=["LANDMARK"])
addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))
addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0]))
# end

v5=addNode!(fg,:x5,2.0*randn(1,N)+getVal(v4)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

v6=addNode!(fg,:x6,1.25*randn(1,N)+getVal(v5)+40.0, N=N,labels=["POSE"])
addFactor!(fg,[v5;v6],Odo([40.0]',[1.25]',[1.0]))

v7=addNode!(fg,:x7,2.0*randn(1,N)+getVal(v6) +60.0, N=N,labels=["POSE"])
addFactor!(fg,[v6;v7],Odo([60.0]',[2.0]',[1.0]))

f3 = addFactor!(fg,[v7], Obsv2(doors, cov', [1.0]))

# Now operate with the data in the DB

tree = wipeBuildNewTree!(fg,drawpdf=true)
# run(`evince bt.pdf`)

# gc()
# recursive solving (single process, easy stack trace for debugging)
println("Starting inference over tree")
inferOverTreeR!(fg, tree)

# Get neighbors
# neighs = CloudGraphs.get_neighbors(fg.cg, cv1r)

gt = Dict{Symbol, Array{Float64,2}}()
gt[:x1]=([0.0;1.97304 ]')' # -0.0342366
gt[:x2]=([50.0; 2.83153 ]')' # 49.8797
gt[:x3]=([100.0; 1.65557 ]')' # 99.8351
gt[:x4]=([150.0; 1.64945 ]')' # 148.637
gt[:x5]=([200.0; 1.77992 ]')' # 198.62
gt[:x6]=([240.0; 2.20466 ]')' # 238.492
gt[:x7]=([300.0; 2.14353 ]')' # 298.467
gt[:l1]=([165.0; 1.17284 ]')' # 164.102

# get vertex back from DB
x1neoID = fg.cgIDs[fg.IDs[:x1]]
cv1r = CloudGraphs.get_vertex(fg.cg, x1neoID, false)
v1r = cloudVertex2ExVertex(cv1r)
plotKDE(getKDE(v1r))

plotKDE( getVertKDE(fg, :x2) )

# draw all beliefs
using Gadfly

DOYTICKS = false
xx,ll = ls(fg)
msgPlots = drawHorBeliefsList(fg, xx, gt=gt,nhor=2);
evalstr = ""
for i in 1:length(msgPlots)
    evalstr = string(evalstr, ",msgPlots[$(i)]")
end
pl = eval(parse(string("vstack(",evalstr[2:end],")")));
Gadfly.draw(PDF("cg4doors.pdf",15cm,20cm),pl) # can also do PNG


text:
using RealSense

include("rs-utils.jl")
include("gl-utils.jl")

# initialize GLFW window
glfwWidth = 1280
glfwHeight = 720
window = create_window_context(glfwWidth, glfwHeight)

# prepare drawing via OpenGL
window_frame = GLfloat[-1.0,-1.0,   1.0,-1.0,   1.0, 1.0,
                        1.0, 1.0,  -1.0, 1.0,  -1.0,-1.0]
vertex_texcoords = GLfloat[0.0,1.0,  1.0,1.0,  1.0,0.0,
                           1.0,0.0,  0.0,0.0,  0.0,1.0]
# create Vertex Buffer Objects
vbo = Ref{GLuint}(0)
glGenBuffers(1, vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo[])
glBufferData(GL_ARRAY_BUFFER, sizeof(window_frame), window_frame, GL_STATIC_DRAW)
texcoordsVBO = Ref{GLuint}(0)
glGenBuffers(1, texcoordsVBO)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_texcoords), vertex_texcoords, GL_STATIC_DRAW)
# create Vertex Array Objects
vao = Ref{GLuint}(0)
glGenVertexArrays(1, vao)
glBindVertexArray(vao[])
glBindBuffer(GL_ARRAY_BUFFER, vao[])
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)
# create shader program
vertSource= joinpath(@__DIR__, "rp.vert")
fragSource = joinpath(@__DIR__, "rp.frag")
shaderProgramID = createprogram(vertSource, fragSource)
# generate texture handle
tex = Ref{GLuint}(0)
glGenTextures(1, tex)

# initialize RealSense camera
err = Ref{Ptr{rs2_error}}(0)
ctx = rs2_create_context(RS2_API_VERSION, err)
checkerror(err)
device_list = rs2_query_devices(ctx, err)
checkerror(err)
dev_count = rs2_get_device_count(device_list, err)
checkerror(err)
@info "There are $dev_count connected RealSense devices."
_pipeline = rs2_create_pipeline(ctx, err)
checkerror(err)
# start the pipeline streaming
@static if is_apple()
    # we sleep 1s on MacOS, otherwise it will hit issue #1586:
    # https://github.com/IntelRealSense/librealsense/issues/1586
    sleep(1)
end
pipeline_profile = rs2_pipeline_start(_pipeline, err)
@assert err[] == C_NULL "Failed to start pipeline!"

# get the current device
dev = rs2_pipeline_profile_get_device(pipeline_profile, err)
checkerror(err)

# create Depth-Colorizer processing block that can be used to quickly visualize the depth data
color_map = rs2_create_colorizer(err)
checkerror(err)
# create target output queue for the results
frame_queue = rs2_create_frame_queue(1, err)
checkerror(err)
# set the output of the processing block to be the queue
rs2_start_processing_queue(color_map, frame_queue, err)
checkerror(err)

# OpenGL rendering loop
recording = false
pause_flag = false
while !GLFW.WindowShouldClose(window)
    updatefps!(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, glfwWidth, glfwHeight)
    glUseProgram(shaderProgramID)

    # wait next frame set and extract depth frame
    frames = rs2_pipeline_wait_for_frames(_pipeline, 5000, err)
    checkerror(err)
    depth_frame = get_depth_frame(frames)
    # select stream
    if depth_frame != C_NULL
        # frame reference is being "moved" into processing, so unless you
        # don't need frame anymore you must add an extra rs2_frame_add_ref
        rs2_frame_add_ref(depth_frame, err)
        checkerror(err)
        # invoke processing
        rs2_process_frame(color_map, depth_frame, err)
        checkerror(err)
        colorized_frame = rs2_wait_for_frame(frame_queue, 5000, err)
        checkerror(err)
        # get data
        depth_data = Ptr{UInt8}(rs2_get_frame_data(colorized_frame, err))
        checkerror(err)
        width = rs2_get_frame_width(depth_frame, err)
        height = rs2_get_frame_height(depth_frame, err)
        glBindTexture(GL_TEXTURE_2D, tex[])
        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, depth_data)
        mipmap()
        glBindVertexArray(vao[])
        glDrawArrays(GL_TRIANGLES, 0, 6)
        # release resource
        rs2_release_frame(colorized_frame)
    end
    # release resource
    rs2_release_frame(depth_frame)
    rs2_release_frame(frames)
    # check and call events
    GLFW.PollEvents()
    # press R to start recording
    if GLFW.GetKey(window, GLFW.KEY_R) && !recording
        # can current device be used as a record device?
        is_extendable = rs2_is_device_extendable_to(dev, RS2_EXTENSION_RECORD, err)
        checkerror(err)
        if is_extendable == 0
            @info "Start recording..."
            record_dev = rs2_create_record_device(dev, joinpath(@__DIR__, "a.bag"), err)
            checkerror(err)
            recording = true
        end
    end
    # press P to pause recording
    if GLFW.GetKey(window, GLFW.KEY_P) && recording && !pause_flag
        @info "Paused"
        rs2_record_device_pause(record_dev, err)
        pause_flag = true
    end
    # press Q to resume recording
    if GLFW.GetKey(window, GLFW.KEY_Q) && recording && pause_flag
        @info "Resumed"
        rs2_record_device_resume(record_dev, err)
        pause_flag = false
    end
    # press F to stop recording
    if GLFW.GetKey(window, GLFW.KEY_F) && recording
        rs2_delete_device(record_dev)
        @info "Recording stopped."
        recording = false
    end
    # swap the buffers
    GLFW.SwapBuffers(window)
end

# release rs2 resources
rs2_delete_frame_queue(frame_queue)
rs2_delete_processing_block(color_map)
rs2_delete_pipeline_profile(pipeline_profile)
rs2_delete_pipeline(_pipeline)
rs2_delete_device(dev)
rs2_delete_device_list(device_list)
rs2_delete_context(ctx)

# release OpenGL resources
GLFW.DestroyWindow(window)


text:
n = rand(1:10)
@test matrixdepot("forsythe", n) ==
                matrixdepot("forsythe", n, sqrt(eps(Float64)), zero(Float64))
@test matrixdepot("forsythe", n , 1, 2) == matrixdepot("forsythe", Float64, n, 1, 2)
T = Float32
a = convert(T, 3)
b = convert(T, 4)
@test matrixdepot("forsythe", T, n, 3, 4) == matrixdepot("forsythe", n, a, b)

if n != 1 # this test does not apply for n = 1
    A = zeros(T, n, n)
    for i = 1:n, j = 1:n
        A[i,j] = j == i ? b :
                  j == i + 1 ? 1.0 :
                  (i == n && j == 1)? a : 0.0
    end

    @test_approx_eq matrixdepot("forsythe", n, a, b) A
end
println("'forsythe' passed test...")


text:
# a belief that just stores the previous observation
# maintained by @zsunberg

# policies based on the previous observation only are often pretty good
# e.g. for the crying baby problem
"""
Updater that stores the most recent observation as the belief, missing if none is available.
"""
struct PreviousObservationUpdater <: Updater end

initialize_belief(u::PreviousObservationUpdater, d::Any) = d

update(bu::PreviousObservationUpdater, old_b, action, obs) = obs


text:
#=
Tools for solving the standard optimal savings / income fluctuation
problem for an infinitely lived consumer facing an exogenous income
process that evolves according to a Markov chain.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-18

References
----------

http://quant-econ.net/jl/ifp.html

=#
using Interpolations
using Optim

# utility and marginal utility functions
u(x) = log(x)
du(x) = 1 ./ x

"""
Income fluctuation problem

##### Fields

- `r::Float64` : Strictly positive interest rate
- `R::Float64` : The interest rate plus 1 (strictly greater than 1)
- `bet::Float64` : Discount rate in (0, 1)
- `b::Float64` :  The borrowing constraint
- `Pi::Matrix{Floa64}` : Transition matrix for `z`
- `z_vals::Vector{Float64}` : Levels of productivity
- `asset_grid::LinSpace{Float64}` : Grid of asset values

"""
type ConsumerProblem
    r::Float64
    R::Float64
    bet::Float64
    b::Float64
    Pi::Matrix{Float64}
    z_vals::Vector{Float64}
    asset_grid::LinSpace{Float64}
end

function ConsumerProblem(;r=0.01, bet=0.96, Pi=[0.6 0.4; 0.05 0.95],
                         z_vals=[0.5, 1.0], b=0.0, grid_max=16, grid_size=50)
    R = 1 + r
    asset_grid = linspace(-b, grid_max, grid_size)

    ConsumerProblem(r, R, bet, b, Pi, z_vals, asset_grid)
end

"""
Given a matrix of size `(length(cp.asset_grid), length(cp.z_vals))`, construct
an interpolation object that does linear interpolation in the asset dimension
and has a lookup table in the z dimension
"""
function Interpolations.interpolate(cp::ConsumerProblem, x::AbstractMatrix)
    sz = (length(cp.asset_grid), length(cp.z_vals))
    if size(x) != sz
        msg = "x must have dimensions $(sz)"
        throw(DimensionMismatch(msg))
    end

    itp = interpolate(x, (BSpline(Linear()), NoInterp()), OnGrid())
    scale(itp, cp.asset_grid, 1:sz[2])
end

"""
Apply the Bellman operator for a given model and initial value.

##### Arguments

- `cp::ConsumerProblem` : Instance of `ConsumerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output
- `;ret_policy::Bool(false)`: Toggles return of value or policy functions

##### Returns

None, `out` is updated in place. If `ret_policy == true` out is filled with the
policy function, otherwise the value function is stored in `out`.

"""
function bellman_operator!(cp::ConsumerProblem, V::Matrix, out::Matrix;
                           ret_policy::Bool=false)
    # simplify names, set up arrays
    R, Pi, bet, b = cp.R, cp.Pi, cp.bet, cp.b
    asset_grid, z_vals = cp.asset_grid, cp.z_vals

    z_idx = 1:length(z_vals)
    vf = interpolate(cp, V)

    # compute lower_bound for optimization
    opt_lb = 1e-8

    # solve for RHS of Bellman equation
    for (i_z, z) in enumerate(z_vals)
        for (i_a, a) in enumerate(asset_grid)

            function obj(c)
                y = 0.0
                for j in z_idx
                    y += vf[R*a+z-c, j] * Pi[i_z, j]
                end
                return -u(c)  - bet * y
            end

            res = optimize(obj, opt_lb, R.*a.+z.+b)
            c_star = res.minimum

            if ret_policy
                out[i_a, i_z] = c_star
            else
               out[i_a, i_z] = - obj(c_star)
            end

        end
    end
    out
end

bellman_operator(cp::ConsumerProblem, V::Matrix; ret_policy=false) =
    bellman_operator!(cp, V, similar(V); ret_policy=ret_policy)

"""
Extract the greedy policy (policy function) of the model.

##### Arguments

- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output

##### Returns

None, `out` is updated in place to hold the policy function

"""
get_greedy!(cp::ConsumerProblem, V::Matrix, out::Matrix) =
    bellman_operator!(cp, V, out, ret_policy=true)

get_greedy(cp::ConsumerProblem, V::Matrix) =
    bellman_operator(cp, V, ret_policy=true)

"""
The approximate Coleman operator.

Iteration with this operator corresponds to policy function
iteration. Computes and returns the updated consumption policy
c.  The array c is replaced with a function cf that implements
univariate linear interpolation over the asset grid for each
possible value of z.

##### Arguments

- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `c::Matrix`: Current guess for the policy function
- `out::Matrix` : Storage for output

##### Returns

None, `out` is updated in place to hold the policy function

"""
function coleman_operator!(cp::ConsumerProblem, c::Matrix, out::Matrix)
    # simplify names, set up arrays
    R, Pi, bet, b = cp.R, cp.Pi, cp.bet, cp.b
    asset_grid, z_vals = cp.asset_grid, cp.z_vals
    z_size = length(z_vals)
    gam = R * bet
    vals = Array(Float64, z_size)

    cf = interpolate(cp, c)

    # linear interpolation to get consumption function. Updates vals inplace
    cf!(a, vals) = map!(i->cf[a, i], vals, 1:z_size)

    # compute lower_bound for optimization
    opt_lb = 1e-8

    for (i_z, z) in enumerate(z_vals)
        for (i_a, a) in enumerate(asset_grid)
            function h(t)
                cf!(R*a+z-t, vals)  # update vals
                expectation = dot(du(vals), vec(Pi[i_z, :]))
                return abs(du(t) - max(gam * expectation, du(R*a+z+b)))
            end
            opt_ub = R*a + z + b  # addresses issue #8 on github
            res = optimize(h, min(opt_lb, opt_ub - 1e-2), opt_ub,
                           method=Optim.Brent())
            out[i_a, i_z] = res.minimum
        end
    end
    out
end

"""
Apply the Coleman operator for a given model and initial value

See the specific methods of the mutating version of this function for more
details on arguments
"""
coleman_operator(cp::ConsumerProblem, c::Matrix) =
    coleman_operator!(cp, c, similar(c))

function init_values(cp::ConsumerProblem)
    # simplify names, set up arrays
    R, bet, b = cp.R, cp.bet, cp.b
    asset_grid, z_vals = cp.asset_grid, cp.z_vals
    shape = length(asset_grid), length(z_vals)
    V, c = Array(Float64, shape...), Array(Float64, shape...)

    # Populate V and c
    for (i_z, z) in enumerate(z_vals)
        for (i_a, a) in enumerate(asset_grid)
            c_max = R*a + z + b
            c[i_a, i_z] = c_max
            V[i_a, i_z] = u(c_max) ./ (1 - bet)
        end
    end

    return V, c
end


text:
using DataFrames
using StatsBase
using CSV
using ProgressMeter
using OnlineStats
using HDF5

"""
    run(model, agglog, bactlog, foodlog; kwargs...)

Runs the given `model`. Aggregate data is logged to file at `agglog` (csv). Data of individual
bacteria is logged to `bactlog` (csv), and the food grid is written to `foodlog` (hdf5).

Keyword arguments
-----------------

- `nsteps = 2000` : Number of steps to run the simulation for
- `log_period = 1` : Periodicity of logging bacterium and food data. Use negative value
  to disable
- `log_bact = true` : Choose whether to log bacteria data
- `log_food = true` : Choose whether to log food map
- `chunk_size = 1000` : How many rows should be written to csv files at a time
- `float_precision = 2` : How many decimals of precision should there be in the csv files
"""
function run!(
    model::ABM,
    agglog::String,
    bactlog::String,
    foodlog::String;
    nsteps::Int = 2000,
    log_period::Int = 1,
    log_bact::Bool = true,
    log_food::Bool = true,
    chunk_size::Int = 1000,
    float_precision::Int = 2,
)
    agg_data = DataFrame(
        step = Int[],
        species = Int[],
        nbact = Int[],
        μ_energy = Float64[],
        σ_energy = Float64[],
        μ_sens = Float64[],
        σ_sens = Float64[],
        μ_repr = Float64[],
        σ_repr = Float64[],
        μ_speed = Float64[],
        σ_speed = Float64[],
    )
    bactlogdf = DataFrame(
        step = Int[],
        age = Int[],
        species = Int[],
        energy = Float64[],
        sens = Float64[],
        repr = Float64[],
        speed = Float64[],
    )

    agglog = open(agglog, "w+")
    log_bact && (bactlog = open(bactlog, "w+"))

    if log_food
        isfile(foodlog) && rm(foodlog)
        foodf = h5open(foodlog, "w")
        foodf["log"] = zeros((floor(Int, nsteps / log_period) + 1, size(model.food)...))
        fdata = HDF5.readmmap(foodf["log"])
    end
    appendbact = false
    appendagg = false

    p = Progress(nsteps; barglyphs = BarGlyphs("[=> ]"), color = :blue)
    for i = 0:nsteps
        nagents(model) == 0 && break
        count = [0 for _ = 1:model.nspecies]
        energy_means = [Mean() for _ = 1:model.nspecies]
        energy_vars = [Variance() for _ = 1:model.nspecies]
        sens_means = [Mean() for _ = 1:model.nspecies]
        sens_vars = [Variance() for _ = 1:model.nspecies]
        repr_means = [Mean() for _ = 1:model.nspecies]
        repr_vars = [Variance() for _ = 1:model.nspecies]
        speed_means = [Mean() for _ = 1:model.nspecies]
        speed_vars = [Variance() for _ = 1:model.nspecies]

        for a in allagents(model)
            count[a.species] += 1
            fit!(energy_means[a.species], a.energy)
            fit!(energy_vars[a.species], a.energy)
            fit!(sens_means[a.species], a.sensory_radius)
            fit!(sens_vars[a.species], a.sensory_radius)
            fit!(repr_means[a.species], a.reproduction_threshold)
            fit!(repr_vars[a.species], a.reproduction_threshold)
            fit!(speed_means[a.species], a.speed)
            fit!(speed_vars[a.species], a.speed)
        end
        for sp = 1:model.nspecies
            push!(
                agg_data,
                (
                    i,
                    sp,
                    count[sp],
                    round(value(energy_means[sp]); digits = float_precision),
                    round(√value(energy_vars[sp]); digits = float_precision),
                    round(value(sens_means[sp]); digits = float_precision),
                    round(√value(sens_vars[sp]); digits = float_precision),
                    round(value(repr_means[sp]); digits = float_precision),
                    round(√value(repr_vars[sp]); digits = float_precision),
                    round(value(speed_means[sp]); digits = float_precision),
                    round(√value(speed_vars[sp]); digits = float_precision),
                ),
            )
        end
        if size(agg_data, 1) >= chunk_size
            CSV.write(agglog, agg_data; append = appendagg)
            appendagg || (appendagg = true)
            empty!(agg_data)
        end

        if log_period > 0 && i % log_period == 0
            if log_bact
                for a in allagents(model)
                    push!(
                        bactlogdf,
                        (
                            i,
                            a.age,
                            a.species,
                            round(a.energy; digits = float_precision),
                            round(a.sensory_radius; digits = float_precision),
                            round(a.reproduction_threshold; digits = float_precision),
                            round(a.speed; digits = float_precision),
                        ),
                    )
                end
                if size(bactlogdf, 1) >= chunk_size
                    CSV.write(bactlog, bactlogdf; append = appendbact)
                    appendbact || (appendbact = true)
                    empty!(bactlogdf)
                end
            end

            log_food && (fdata[floor(Int, i / log_period)+1, :, :] = model.food)
        end
        Agents.step!(model, agent_step!, food_step!)
        ProgressMeter.next!(p; showvalues = [(:steps, i), (:population, nagents(model))])
    end

    size(agg_data, 1) > 0 && CSV.write(agglog, agg_data; append = appendagg)
    size(bactlogdf, 1) > 0 && CSV.write(bactlog, bactlogdf; append = appendbact)

    if log_food
        foodf["cap"] = model.food_data.food_cap
        close(foodf)
    end
    nothing
end

"""
    runconfig(config, dirname, foodpath; kwargs...)

Create a model with given `config`, run it using `foodmap` as the initial food distribution,
and log data to `dirname`.

Keyword arguments
-----------------
- See `run!`
- `saveconfig = false` whether the given config file should also be saved to `dirname` as BSON
"""
function runconfig(
    config,
    dirname::String,
    foodpath::String;
    nsteps::Int = 2000,
    log_period::Int = 10,
    log_bact::Bool = true,
    log_food::Bool = true,
    chunk_size::Int = 1000,
    float_precision::Int = 2,
    saveconfig::Bool = false,
)
    haskey(config, :initial_bacterium) &&
        config[:initial_bacterium] isa Vector &&
        (config[:initial_bacterium] = Vector{Bacterium}(config[:initial_bacterium]))
    mkpath(dirname)
    close(open(joinpath(dirname, "$(configname(config))"), "w"))

    food = image_to_foodmap(foodpath)
    model = initialize_model(food .* config[:food_data].food_cap; config...)

    run!(
        model,
        joinpath(dirname, "agg.csv"),
        joinpath(dirname, "bact.csv"),
        joinpath(dirname, "food.hdf5");
        nsteps,
        log_period,
        log_bact,
        log_food,
        chunk_size,
        float_precision,
    )
    saveconfig && bson(joinpath(dirname, "config.bson"), config)
    return model
end

"""
    configname(config)

Generate a highly abbreviated name for a config
"""
function configname(config::Dict)
    name = ""
    if haskey(config, :n_bacteria)
        if config[:n_bacteria] isa Vector
            name *= "nbx$(length(config[:n_bacteria]))"
        else
            name *= "nb$(config[:n_bacteria])"
        end
    end

    if haskey(config, :initial_bacterium)
        name *= "ib("
        if config[:initial_bacterium] isa Vector
            name *= "e"
            for b in config[:initial_bacterium]
                name *= "$(b.energy)_"
            end
            name *= "s"
            for b in config[:initial_bacterium]
                name *= "$(b.speed)_"
            end
            name *= "sr"
            for b in config[:initial_bacterium]
                name *= "$(b.sensory_radius)_"
            end
            name *= "r"
            for b in config[:initial_bacterium]
                name *= "$(b.reproduction_threshold)_"
            end
        else
            b = config[:initial_bacterium]
            name *= "e$(b.energy)s$(b.speed)sr$(b.sensory_radius)r$(b.reproduction_threshold)"
        end
        name *= ")"
    end

    if haskey(config, :food_data)
        f = config[:food_data]
        name *= "fd(c$(f.food_cap)r$(f.regen_rate))"
    end

    haskey(config, :eat_rate_factor) && (name *= "erf$(config[:eat_rate_factor])")
    haskey(config, :lifetime) && (name *= "l$(config[:lifetime])")
end

max_stars_repo_name: jiaangou/SimpleEvolution

max_issues_repo_name: jiaangou/SimpleEvolution

max_forks_repo_name: jiaangou/SimpleEvolution


text:
struct NexGaugePointer <: AbstractNexObject
    _nid::NexID

    widthable::IWidthable
end

function setproperty!(obj::NexGaugePointer, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

"""
    NexGauge(nexSerial, name; pid=pid, cid=cid)

A Nextion Gauge UI control.

On Nextion side, it have the following properties:
- bco
- val
- wid
- pco
"""
struct NexGauge <: AbstractNexObject
    _nid::NexID

    viewable::IViewable
    numericalvalued::INumericalValued
    colourable::IColourable
    touchable::ITouchable

    function NexGauge(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid, IViewable(nid), INumericalValued(nid, UInt16, RangeNumber{UInt16, 0:360}), IColourable(nid), ITouchable(nid))
    end
end

function setproperty!(obj::NexGauge, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexGauge, property::Symbol)
    if property == :pointer
        nid = NexID(obj)
        NexGaugePointer(nid, IWidthable(nid))
    else
        _getcommonproperty(obj, property)
    end
end


text:
######################################################################
# General purpose functionized patterns
# -----
# Licensed under MIT License.
import Base.Threads: @threads

export isinstanceof
isinstanceof(x::T2, ::Type{T1}) where {T1, T2<:T1} = true
isinstanceof(_, _) = false

export hassignature
function hassignature(fn, argtypes::Type...)
    sig = Tuple{typeof(fn), argtypes...}
    for method ∈ methods(fn)
        if sig <: method.sig
            return true
        end
    end
    return false
end

export isiterable
@generated function isiterable(x)
    if hassignature(iterate, x)
        :(true)
    else
        :(false)
    end
end

export curry
curry(fn, curryargs...; kwcurryargs...) = (moreargs...; kwargs...) -> fn(curryargs..., moreargs...; kwcurryargs..., kwargs...)

export decamelcase
"""`decamelcase(str; uppercase::Bool = false)`
Remove camel-case and transform into lowercase or uppercase with underscore separators (e.g. "fooBarBaz" to "foo_bar_baz").
If `uppercase` is true, the result will be uppercased, otherwise lowercased."""
function decamelcase(str::AbstractString; uppercase::Bool = false)
    casetransform = uppercase ? Base.uppercase : lowercase
    casetransform(replace(str, r"(\p{Ll})(\p{Lu})" => s"\1_\2"))
end

export indexof
indexof(arr, elem; by = identity, offset::Integer = 1, strict::Bool = false) = findnext(curr->strict ? by(curr) === elem : by(curr) == elem, arr, offset)

export shift!
function shift!(arr)
    item = arr[1]
    deleteat!(arr, 1)
    item
end

export unshift!
function unshift!(arr, value)
    insert!(arr, 1, value)
    arr
end

export hflip!
function hflip!(arr::AbstractArray)
    rows, cols = size(arr)
    if cols < 2 return arr end

    @threads for row ∈ 1:rows
        @threads for col ∈ 1:(cols÷2)
            tmp = arr[row, col]
            arr[row, col] = arr[row, cols-(col-1)]
            arr[row, cols-(col-1)] = tmp
        end
    end
    arr
end

export vflip!
function vflip!(arr::AbstractArray)
    rows, cols = size(arr)
    if rows < 2 return arr end

    @threads for col ∈ 1:cols
        @threads for row ∈ 1:(rows÷2)
            tmp = arr[row, col]
            arr[row, col] = arr[rows-(row-1), col]
            arr[rows-(row-1), col] = tmp
        end
    end
    arr
end

is_call_expr(_) = false
is_call_expr(expr::Expr) = expr.head === :call
is_call_assignment(expr::Expr) = expr.head === :(=) && is_call_expr(expr.args[2])

function Base.wait(test, cond::Condition; timeout::Number = Inf)
    value = nothing

    if timeout === Inf
        while !test()
            value = wait(cond)
        end
    else
        istimeout = false
        timer = Timer(timeout)
        @async begin
            wait(timer)
            istimeout = true
            notify(cond)
        end

        while !test() && !istimeout
            value = @await wait(cond)
        end

        close(timer)
        if istimeout
            throw(TimeoutError())
        end
    end

    return value
end


text:
include("../examples/smiley_examples.jl")
using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55

using BenchmarkTools
using ForwardDiff
using IntervalArithmetic
using NumberIntervals
using Turmeric
using StaticArrays
import Random

const SUITE = BenchmarkGroup()

Random.seed!(0)  # Seed the RNG to get consistent results
const tol = 1e-6
const methods = (
        Newton(),
        Krawczyk(),
    )

S = SUITE["Smiley"] = BenchmarkGroup()
for example in (SmileyExample22, SmileyExample52, SmileyExample54, SmileyExample55)
    s = S[example.title] = BenchmarkGroup()
    for method in methods
        s[string(method)] = @benchmarkable roots(
            $(example.f),
            $(NumberInterval.(example.region)),
            $method,
            $tol
        )
    end
end

S = SUITE["Rastrigin stationary points"] = BenchmarkGroup()
include("dietmar_ratz_functions.jl")

# Rastrigin function:
const A = 10

f(x, y) = 2A + x^2 - A * cos(2π * x) + y^2 - A * cos(2π * y)
f(X) = f(X[1], X[2])

ForwardDiff.gradient(f, X::IntervalBox) = ForwardDiff.gradient(f, X.v)
∇(f) = X -> ForwardDiff.gradient(f, X)

L = 5.0
X = NumberInterval.(SVector(-L..(L+1), -L..(L+1)))

for method in methods
    S[string(method)] = @benchmarkable roots($(∇(f)), $X, $method, 1e-5)
end

S = SUITE["Dietmar-Ratz"] = BenchmarkGroup()
X = NumberInterval(0.75, 1.75)

for (k, dr) in enumerate(dr_functions)
    s = S["Dietmar-Ratz $k"] = BenchmarkGroup()

    if k != 8  # dr8 is excluded as it has too many roots
        for method in methods
            s[string(method)] = @benchmarkable roots($dr, $X, $method, $tol)
        end
    end

    s["Automatic differentiation"] = @benchmarkable ForwardDiff.derivative($dr, $X)
end