text:
function bootstrap(data::AbstractArray,d::UnivariateDistribution;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  n = length(data)

  seed == 0 ? srand() : srand(seed)
  statistic = Array(Tuple{typeof(d),Float64},no_of_sims)
  for i=1:no_of_sims
    sim_data = sample(data, n, replace=true)
    statistic[i] =estimate_xmin(sim_data,typeof(d),xmins = xmins,xmax = xmax)
  end
  return statistic
end

function bootstrap(data::AbstractArray,distribution::Type{con_powerlaw};no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap(data,d,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end

function bootstrap(data::AbstractArray,distribution::Type{dis_powerlaw};no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap(data,d,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end

function bootstrap_p(data::AbstractArray,d::UnivariateDistribution;no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  sort_data = sort(data)
  α,θ = params(d)
  n = length(sort_data)
  tail_indx = findfirst(sort_data,θ)
  tail_p = length(sort_data[tail_indx:end])/n
  KS_stat = Kolmogorov_smirnov_test(sort_data[tail_indx:end],d)

  P = 0
  statistic = Array(Tuple{typeof(d),Float64},no_of_sims)
  seed == 0 ? srand() : srand(seed)
  for i=1:no_of_sims
    n1 = sum(map(x-> x>tail_p,rand(n)))
    n2 = n - n1
    sim_data = Array(Float64,0)
    append!(sim_data,sample(sort_data[1:tail_indx-1],n1,replace = true))
    append!(sim_data,rand(d,n2))
    statistic[i] = estimate_xmin(sim_data,typeof(d),xmins = xmins,xmax = xmax)
    if (KS_stat <= statistic[i][2])
      P +=1
    end
  end
  return statistic,(P/no_of_sims)
end

function bootstrap_p(data::AbstractArray,distribution::Type{con_powerlaw};no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap_p(data,d,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end

function bootstrap_p(data::AbstractArray,distribution::Type{dis_powerlaw};no_of_sims::Int64 = 10,xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5),seed::Int64 = 0)
  d,ks = estimate_xmin(data,distribution,xmins = xmins,xmax =xmax)
  bootstrap_p(data,d,no_of_sims = no_of_sims,xmins = xmins,xmax = xmax,seed =seed)
end


text:
export RegionGraph, random_region_graph, region_graph_2_pc
import Random: shuffle

"Root of region graph node hierarchy"
abstract type RegionGraph <: Tree end

######
## Partition
#####
mutable struct Partition
    partition::AbstractVector{T} where T <: RegionGraph
    Partition(n::AbstractVector{T}) where T <: RegionGraph = begin
        new(n)
    end
end

import Base: getindex, size, length
@inline Base.size(n::Partition) = size(n.partition)
@inline Base.length(n::Partition) = length(n.partition)
@inline Base.getindex(n::Partition, i::Int) = n.partition[i]

mutable struct RegionGraphInnerNode <: RegionGraph
    partitions::AbstractVector{Partition}
    variables::BitSet
    parent::Union{Nothing, RegionGraphInnerNode}

    RegionGraphInnerNode(partitions::AbstractVector{Partition}) = begin
        for partition in partitions
            for i = 1:length(partition), j = i+1 : length(partition)
                @assert isdisjoint(variables(partition[i]), variables(partition[j]))
            end
        end
        scope = variables(partitions[1][1])
        for ind = 2:size(partitions[1])[1]
            scope = scope ∪ variables(partitions[1][ind])
        end

        this = new(partitions, scope, nothing)
        for partition in partitions
            for i = 1:size(partition)[1]
                @assert isnothing(partition[i].parent)
                partition[i].parent = this
            end
        end
        this
    end
end

mutable struct RegionGraphLeafNode <: RegionGraph
    variables::BitSet
    parent::Union{Nothing, RegionGraphInnerNode}
    RegionGraphLeafNode(v) = new(v, nothing)
end

#####################
# Constructors
#####################

RegionGraph(v::AbstractVector{Var}) = RegionGraphLeafNode(BitSet(v))
RegionGraph(partitions::AbstractVector{Partition}) = RegionGraphInnerNode(partitions)

#####################
# Traits
#####################

import LogicCircuits.NodeType
@inline LogicCircuits.NodeType(::RegionGraphInnerNode) = Inner()
@inline LogicCircuits.NodeType(::RegionGraphLeafNode) = Leaf()

###################
## Methods
#################

import DirectedAcyclicGraphs: children

@inline children(n::RegionGraphInnerNode) = n.partitions

import LogicCircuits.Utils: variables

@inline variables(n::RegionGraphInnerNode) = n.variables
@inline variables(n::RegionGraphLeafNode) = n.variables

##################################################################

"""
    random_region_graph(X::AbstractVector{Var}, depth::Int = 5, replicas::Int = 2, num_splits::Int = 2)

- `X`: Vector of all variables to include; for the root region
- `depth`: how many layers to do splits
- `replicas`: number of replicas or paritions (replicas only used for the root region; for other regions only 1 parition (inner nodes), or 0 parition for leaves)
- `num_splits`: number of splits for each parition; split variables into random equaly sized regions
"""
function random_region_graph(X::AbstractVector{Var};
    depth::Int = 5, replicas::Int = 2, num_splits::Int = 2)::RegionGraph

    if length(X) < 2 || depth == 0
        # Cannot/should not split anymore
        return RegionGraph(X)
    end

    partitions = Vector{Partition}()
    for repeat = 1 : replicas
        cur_rg = split_rg(X, depth; num_splits=num_splits)
        push!(partitions, cur_rg)
    end

    # Validation: Each Partition should include the same set of variables
    prev_scope = nothing
    for cur_partition in partitions
        cur_scope = variables(cur_partition[1])
        for i = 2:length(cur_partition)
            cur_scope = cur_scope ∪ variables(cur_partition[i])
        end

        if !isnothing(prev_scope)
            @assert prev_scope == cur_scope "All partitions should include the same set of variables."
        else
            prev_scope = cur_scope
        end
    end

    RegionGraphInnerNode(partitions)
end

function split_rg(variables::AbstractVector{Var}, depth::Int; num_splits::Int = 2)::Partition
    partition(x, n) = begin
        ## TODO; might not work if num_splits > 2
        sz = ceil(Int, (length(x) / n))
        [ x[i:min(i+sz-1, length(x))] for i = 1:sz:length(x) ]
    end

    shuffle_variables = shuffle(variables)
    splits = partition(shuffle_variables, num_splits)
    cur_partition_regions = Vector{RegionGraph}()
    for split in splits
        # only 1 replicas for non-root
        rg_node = random_region_graph(split; depth=depth-1, replicas=1, num_splits)
        push!(cur_partition_regions, rg_node)
    end
    Partition(cur_partition_regions)
end

########################################
### Rat-SPNs in Juice PC data structure
########################################

"""
Makes sure the sum nodes does not have too many children. Makes balanced sums of sums to reduce children count.
"""
function balance_sum(children::Vector{ProbCircuit}, balance_childs_parents)::PlainSumNode
    if balance_childs_parents
        if length(children) <= 4
            PlainSumNode(children)
        else
            ls = 1:floor(Int32, length(children) / 2)
            rs = ceil(Int32, length(children) / 2):length(children)
            PlainSumNode([balance_sum(children[ls], balance_childs_parents), balance_sum(children[rs], balance_childs_parents)])
        end
    else
        PlainSumNode(children)
    end
end

"""
Makes sure input nodes don't have too many parents.
Makes a dummy sum node for each input per partition. Then nodes corresponding to the partition use
the dummy node as their children instead of the input node.
This way instead of num_nodes_root * num_nodes_leaf, we would have num_nodes_root parents nodes.
"""
function balanced_fully_factorized(vtrees::Vector{<:Vtree})::Vector{ProbCircuit}
    vars = variables(vtrees[1]) # assuming all vtrees same set of variables
    var_2_dummy_inputs = Dict(var => PlainSumNode([PlainProbLiteralNode(var), PlainProbLiteralNode(-var)]) for var in vars)

    balanced_recurse(nodes::Vector{<:Vtree})::Vector{ProbCircuit} = begin
        if isleaf(nodes[1])
            [PlainSumNode([var_2_dummy_inputs[variable(node)]]) for node in nodes]
        else
            lefts = balanced_recurse([node.left for node in nodes])
            rights = balanced_recurse([node.right for node in nodes])
            [PlainSumNode([PlainMulNode([left, right])]) for (left, right) in zip(lefts, rights)]
        end
    end

    balanced_recurse(vtrees)
end

"""
    region_graph_2_pc(node::RegionGraph; num_nodes_root, num_nodes_region, num_nodes_leaf, balance_childs_parents)

- `num_nodes_root`: number of sum nodes in the root region
- `num_nodes_leaf`: number of sum nodes per leaf region
- `num_nodes_region`: number of in each region except root and leaves
"""
function region_graph_2_pc(node::RegionGraph;
    num_nodes_root::Int = 4,  num_nodes_region::Int = 3, num_nodes_leaf::Int = 2, balance_childs_parents=true)::Vector{ProbCircuit}
    sum_nodes = Vector{ProbCircuit}()
    if isleaf(node)
        if false # balance_childs_parents # TODO fix, turning this on makes param learning faster but log likelihood improves much slower
            vtrees = [Vtree([PlainVtreeLeafNode(x) for x ∈ variables(node)], :balanced) for i=1:num_nodes_leaf]
            sum_nodes = balanced_fully_factorized(vtrees)
        else
            for i = 1:num_nodes_leaf
                vtree = Vtree([PlainVtreeLeafNode(x) for x ∈ variables(node)], :balanced)
                leaf_circuit = PlainProbCircuit(fully_factorized_circuit(ProbCircuit, vtree))
                push!(sum_nodes, leaf_circuit)
            end
        end
    else
        root_children = Vector{ProbCircuit}()
        # Foreach replication; usually only > 1 at root
        for partition in node.partitions
            partition_mul_nodes = Vector{ProbCircuit}()

            @assert length(partition) == 2 "Only supporting partitions of size 2 at the moment"
            lefts = region_graph_2_pc(partition[1]; num_nodes_root, num_nodes_region, num_nodes_leaf, balance_childs_parents)
            rights = region_graph_2_pc(partition[2]; num_nodes_root, num_nodes_region, num_nodes_leaf, balance_childs_parents)
            @assert all([issum(l) for l in lefts])
            @assert all([issum(r) for r in rights])

            for l in lefts, r in rights
                mul_node = PlainMulNode([l, r])
                push!(partition_mul_nodes, mul_node)
            end

            dummy_sum_node = balance_sum(partition_mul_nodes, balance_childs_parents)
            push!(root_children, dummy_sum_node)
        end

        # Repeat Sum nodes nodes based on where in Region Graph
        if isnothing(node.parent)
            # Root region
            for i = 1:num_nodes_root
                sum_node = balance_sum(root_children, balance_childs_parents)
                push!(sum_nodes, sum_node)
            end
        else
            # Inner region
            for i = 1:num_nodes_region
                sum_node = balance_sum(root_children, balance_childs_parents)
                push!(sum_nodes, sum_node)
            end
        end
    end

    sum_nodes
end

ext: jl


text:
module Assertions

export activate_assertions,
       deactivate_assertions,
       are_assertions_enabled

# source: https://discourse.julialang.org/t/defensive-programming-assert/8383/11

# functions to (de)activate assertions
function activate_assertions(m::Module)
    m.eval(:(are_assertions_enabled() = true))
    nothing
end
function deactivate_assertions(m::Module)
    m.eval(:(are_assertions_enabled() = false))
    nothing
end

# override Base.@assert
macro assert(exs...)
    quote
        if $__module__.are_assertions_enabled()
            Base.@assert $(map(esc, exs)...)
        end
    end
end

end  # module


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

@doc raw"""PodDisruptionBudgetSpec is a description of a PodDisruptionBudget.

    IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec(;
        maxUnavailable=nothing,
        minAvailable=nothing,
        selector=nothing,
    )

    - maxUnavailable::IoK8sApimachineryPkgUtilIntstrIntOrString : An eviction is allowed if at most \"maxUnavailable\" pods selected by \"selector\" are unavailable after the eviction, i.e. even in absence of the evicted pod. For example, one can prevent all voluntary evictions by specifying 0. This is a mutually exclusive setting with \"minAvailable\".
    - minAvailable::IoK8sApimachineryPkgUtilIntstrIntOrString : An eviction is allowed if at least \"minAvailable\" pods selected by \"selector\" will still be available after the eviction, i.e. even in the absence of the evicted pod.  So for example you can prevent all voluntary evictions by specifying \"100%\".
    - selector::IoK8sApimachineryPkgApisMetaV1LabelSelector : Label query over pods whose evictions are managed by the disruption budget.
"""
mutable struct IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec <: SwaggerModel
    maxUnavailable::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: maxUnavailable
    minAvailable::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: minAvailable
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector

    function IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec(;maxUnavailable=nothing, minAvailable=nothing, selector=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec, Symbol("maxUnavailable"), maxUnavailable)
        setfield!(o, Symbol("maxUnavailable"), maxUnavailable)
        validate_property(IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec, Symbol("minAvailable"), minAvailable)
        setfield!(o, Symbol("minAvailable"), minAvailable)
        validate_property(IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        o
    end
end # type IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec

const _property_map_IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec = Dict{Symbol,Symbol}(Symbol("maxUnavailable")=>Symbol("maxUnavailable"), Symbol("minAvailable")=>Symbol("minAvailable"), Symbol("selector")=>Symbol("selector"))
const _property_types_IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec = Dict{Symbol,String}(Symbol("maxUnavailable")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("minAvailable")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec[property_name]

function check_required(o::IoK8sApiPolicyV1beta1PodDisruptionBudgetSpec)
    true
end

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


text:
cipher_ascii = [
  79,59,12,2,79,35,8,28,20,2,3,68,8,9,68,45,0,12,9,67,68,4,7,5,23,27,1,21,79,85,78,79,85,71,38,10,71,27,12,2,79,6,2,8,13,9,1,13,9,8,68,19,7,1,71,
  56,11,21,11,68,6,3,22,2,14,0,30,79,1,31,6,23,19,10,0,73,79,44,2,79,19,6,28,68,16,6,16,15,79,35,8,11,72,71,14,10,3,79,12,2,79,19,6,28,68,32,0,0,
  73,79,86,71,39,1,71,24,5,20,79,13,9,79,16,15,10,68,5,10,3,14,1,10,14,1,3,71,24,13,19,7,68,32,0,0,73,79,87,71,39,1,71,12,22,2,14,16,2,11,68,2,25,
  1,21,22,16,15,6,10,0,79,16,15,10,22,2,79,13,20,65,68,41,0,16,15,6,10,0,79,1,31,6,23,19,28,68,19,7,5,19,79,12,2,79,0,14,11,10,64,27,68,10,14,15,2,
  65,68,83,79,40,14,9,1,71,6,16,20,10,8,1,79,19,6,28,68,14,1,68,15,6,9,75,79,5,9,11,68,19,7,13,20,79,8,14,9,1,71,8,13,17,10,23,71,3,13,0,7,16,71,
  27,11,71,10,18,2,29,29,8,1,1,73,79,81,71,59,12,2,79,8,14,8,12,19,79,23,15,6,10,2,28,68,19,7,22,8,26,3,15,79,16,15,10,68,3,14,22,12,1,1,20,28,72,
  71,14,10,3,79,16,15,10,68,3,14,22,12,1,1,20,28,68,4,14,10,71,1,1,17,10,22,71,10,28,19,6,10,0,26,13,20,7,68,14,27,74,71,89,68,32,0,0,71,28,1,9,27,
  68,45,0,12,9,79,16,15,10,68,37,14,20,19,6,23,19,79,83,71,27,11,71,27,1,11,3,68,2,25,1,21,22,11,9,10,68,6,13,11,18,27,68,19,7,1,71,3,13,0,7,16,71,
  28,11,71,27,12,6,27,68,2,25,1,21,22,11,9,10,68,10,6,3,15,27,68,5,10,8,14,10,18,2,79,6,2,12,5,18,28,1,71,0,2,71,7,13,20,79,16,2,28,16,14,2,11,9,
  22,74,71,87,68,45,0,12,9,79,12,14,2,23,2,3,2,71,24,5,20,79,10,8,27,68,19,7,1,71,3,13,0,7,16,92,79,12,2,79,19,6,28,68,8,1,8,30,79,5,71,24,13,19,1,
  1,20,28,68,19,0,68,19,7,1,71,3,13,0,7,16,73,79,93,71,59,12,2,79,11,9,10,68,16,7,11,71,6,23,71,27,12,2,79,16,21,26,1,71,3,13,0,7,16,75,79,19,15,0,
  68,0,6,18,2,28,68,11,6,3,15,27,68,19,0,68,2,25,1,21,22,11,9,10,72,71,24,5,20,79,3,8,6,10,0,79,16,8,79,7,8,2,1,71,6,10,19,0,68,19,7,1,71,24,11,21,
  3,0,73,79,85,87,79,38,18,27,68,6,3,16,15,0,17,0,7,68,19,7,1,71,24,11,21,3,0,71,24,5,20,79,9,6,11,1,71,27,12,21,0,17,0,7,68,15,6,9,75,79,16,15,10,
  68,16,0,22,11,11,68,3,6,0,9,72,16,71,29,1,4,0,3,9,6,30,2,79,12,14,2,68,16,7,1,9,79,12,2,79,7,6,2,1,73,79,85,86,79,33,17,10,10,71,6,10,71,7,13,20,
  79,11,16,1,68,11,14,10,3,79,5,9,11,68,6,2,11,9,8,68,15,6,23,71,0,19,9,79,20,2,0,20,11,10,72,71,7,1,71,24,5,20,79,10,8,27,68,6,12,7,2,31,16,2,11,
  74,71,94,86,71,45,17,19,79,16,8,79,5,11,3,68,16,7,11,71,13,1,11,6,1,17,10,0,71,7,13,10,79,5,9,11,68,6,12,7,2,31,16,2,11,68,15,6,9,75,79,12,2,79,
  3,6,25,1,71,27,12,2,79,22,14,8,12,19,79,16,8,79,6,2,12,11,10,10,68,4,7,13,11,11,22,2,1,68,8,9,68,32,0,0,73,79,85,84,79,48,15,10,29,71,14,22,2,79,
  22,2,13,11,21,1,69,71,59,12,14,28,68,14,28,68,9,0,16,71,14,68,23,7,29,20,6,7,6,3,68,5,6,22,19,7,68,21,10,23,18,3,16,14,1,3,71,9,22,8,2,68,15,26,
  9,6,1,68,23,14,23,20,6,11,9,79,11,21,79,20,11,14,10,75,79,16,15,6,23,71,29,1,5,6,22,19,7,68,4,0,9,2,28,68,1,29,11,10,79,35,8,11,74,86,91,68,52,0,
  68,19,7,1,71,56,11,21,11,68,5,10,7,6,2,1,71,7,17,10,14,10,71,14,10,3,79,8,14,25,1,3,79,12,2,29,1,71,0,10,71,10,5,21,27,12,71,14,9,8,1,3,71,26,23,
  73,79,44,2,79,19,6,28,68,1,26,8,11,79,11,1,79,17,9,9,5,14,3,13,9,8,68,11,0,18,2,79,5,9,11,68,1,14,13,19,7,2,18,3,10,2,28,23,73,79,37,9,11,68,16,
  10,68,15,14,18,2,79,23,2,10,10,71,7,13,20,79,3,11,0,22,30,67,68,19,7,1,71,8,8,8,29,29,71,0,2,71,27,12,2,79,11,9,3,29,71,60,11,9,79,11,1,79,16,15,
  10,68,33,14,16,15,10,22,73
]

str2int(s) = Int(s[1])

plain_text = ""
for i='a':'z', j='a':'z', k='a':'z'
  global cipher_ascii, plain_text
  println("Trying ", i, j, k)
  key = [Int(i), Int(j), Int(k)]
  plain_text = ""
  for n = 1:length(cipher_ascii)
    plain_text = string(plain_text, Char(cipher_ascii[n] ⊻ key[(n % 3) + 1]))
  end
  if occursin(" the ", plain_text)
    println(plain_text)
    println(reduce(+, str2int.(split(plain_text, ""))))
    break
  end
end

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff


text:
using Test
using QuantumOptics, CollectiveSpins

@testset "system" begin

# Make sure that the interface works as expected
spin1 = CollectiveSpins.Spin([0,0,0], delta=2)
spin2 = CollectiveSpins.Spin([1.2,0,0], delta=-1.)

S = SpinCollection([spin1, spin2], [0, 0, 1]; gammas=0.1)

S = SpinCollection(CollectiveSpins.geometry.chain(0.2, 4), [0., 1., 1.]; deltas=[1., 1., 1., 1.], gammas=2.)

cavitymode = CavityMode(1, delta=2, eta=3, kappa=4.)

CavitySpinCollection(cavitymode, S, 1)
CavitySpinCollection(cavitymode, S, [1, 2, 3, 4])

end # testset


text:
#=
Popis
  Funkcia vypocita pre kazdy vrchol hodnotu jeho closeness centrality.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozena vlastnost vrchola :: String
  normalize : parameter pre aplikovanie normalizacie :: Bool

Navratova hodnota
  closeness : slovnik obsahujuci mapovanie vrcholov a k nim prisluchajucich hodnot
           closeness centrality (index_vrchola => hodnota_closeness_centrality) :: Dict{Int64, Float64}()

Poznamky
  Hodnota centrality moze byt normalizovana nastavenim parametra normalize na True/False.
  Vypocet centrality je upraveny aby velkost grafu neovplyvnovala hodnotu centrality.
=#

function centrality_closeness(g::GenericGraph; dist_key::AbstractString="weight", normalize::Bool=true)

  epi = AttributeEdgePropertyInspector{Float64}(dist_key)
  closeness = Dict{Int64, Float64}()

  for v in Graphs.vertices(g)

      distances = filter(x -> x != Inf, dijkstra_shortest_paths(g, epi, [v]).dists)
      total_dist = sum(distances)

      if total_dist > 0 && num_vertices(g) > 1

        closeness[v.index] = (length(distances) - 1.0) / total_dist

        if normalize

          norm_coefficient = (length(distances)-1.0) / ( num_vertices(g) - 1 )
          closeness[v.index] *= norm_coefficient

        end

      end

  end

  return closeness

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
## test data available for use in PopGen
export dataset, nancycats, @nancycats, gulfsharks, @gulfsharks

"""
    dataset(::String)
Load an example dataset from either `"gulfsharks"` (SNP) or `"nancycats"` (microsatellite). Can also use `"sharks"` and `"cats"`
as shorthands. Use `?nancycats` and `?gulfsharks` to learn more about
these datasets.

### Example
```
ncats = dataset("nancycats")
gsharks = dataset("sharks")
```
"""
function dataset(name::String)
    lowercase(name) in ["gulfsharks", "sharks"] && return PopGen._gulfsharks()
    lowercase(name) in ["nancycats", "cats"] && return PopGen._nancycats()
end

"""
    @nancycats
Returns a `PopObj` of corresponding "nancycats" dataset as featured in
the R package `adegenet`. This is microsatellite data of 9 loci in 237
individuals across 17 populations.

Example:
```
ncats = @nancycats
```
"""
function nancycats()
    @warn "This method is deprecated, please use @nancycats instead"  * "\nnancycats() will be removed in the next release"
    filename = normpath(joinpath(@__DIR__,"../..","data", "datasets.jld2"))
    load(filename, "nancycats")
end

function _nancycats()
    filename = normpath(joinpath(@__DIR__,"../..","data", "datasets.jld2"))
    load(filename, "nancycats")
end

macro nancycats()
    return esc(quote
       PopGen._nancycats()
    end)
end

"""
    @gulfsharks
Returns a `PopObj` corresponding the Blacknose shark dataset as used in
Dimens et al. 2019. This is a mid-sized SNV dataset of 2209 haplotypes
in 212 individuals, across 7 populations.

Example:

```
sharks = @gulfsharks
```
"""
function gulfsharks()
    @warn "This method is deprecated, please use @gulfsharks instead" * "\ngulfsharks() will be removed in the next release"
    filename = normpath(joinpath(@__DIR__,"../..","data", "datasets.jld2"))
    load(filename, "gulfsharks")
end

function _gulfsharks()
    filename = normpath(joinpath(@__DIR__,"../..","data", "datasets.jld2"))
    load(filename, "gulfsharks")
end

macro gulfsharks()
    return esc(quote
       PopGen._gulfsharks()
    end)
end


text:
@testset "$(rpad("Tabulated Quadrature Rules",80))" begin

    @test typeof(RiemannQuadratureLeft()) <: QuadratureRule
    @test typeof(RiemannQuadratureRight()) <: QuadratureRule
    @test typeof(MidpointQuadrature()) <: QuadratureRule
    @test typeof(TrapezoidalQuadrature()) <: QuadratureRule

end


text:
# module for representing motion edge

using LinearAlgebra

include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))

mutable struct MotionEdge
  t1
  t2
  x1
  x2
  omega_upper_right
  omega_upper_left
  omega_lower_right
  omega_lower_left
  xi_upper
  xi_lower

  # matrix of input noise
  function mat_M(speed, yaw_rate, time, stds)
    return diagm(0 => [stds["nn"]^2*abs(speed)/time + stds["no"]^2*abs(yaw_rate)/time,
                       stds["on"]^2*abs(speed)/time + stds["oo"]^2*abs(yaw_rate)/time])
  end

  # jacobian of input
  function mat_A(speed, yaw_rate, time, theta)
    st, ct = sin(theta), cos(theta)
    stw, ctw = sin(theta + yaw_rate * time), cos(theta + yaw_rate * time)
    return [(stw - st)/yaw_rate  -speed/(yaw_rate^2)*(stw - st) + speed/yaw_rate*time*ctw;
            (-ctw + ct)/yaw_rate -speed/(yaw_rate^2)*(-ctw + ct) + speed/yaw_rate*time*stw;
            0                    time]
  end

  # jacobian of pose
  function mat_F(speed, yaw_rate, time, theta)
    F = diagm(0 => [1.0, 1.0, 1.0])
    F[1, 3] = speed / yaw_rate * (cos(theta + yaw_rate * time) - cos(theta))
    F[2, 3] = speed / yaw_rate * (sin(theta + yaw_rate * time) - sin(theta))
    return F
  end

  # init
  function MotionEdge(t1, t2, pose_list, input_list, delta;
                      motion_noise_stds=Dict("nn"=>0.20, "no"=>0.001, "on"=>0.11, "oo"=>0.20))
    self = new()

    self.t1 = t1 # time
    self.t2 = t2
    self.x1 = pose_list[string(t1)] # pose [x, y, theta]
    self.x2 = pose_list[string(t2)]
    speed = input_list[string(t2)][1] # input [speed, yaw_rate]
    yaw_rate = input_list[string(t2)][2]
    if abs(yaw_rate) < 1e-5
      yaw_rate = 1e-5
    end

    # noise matrix and jacobian
    M = mat_M(speed, yaw_rate, delta, motion_noise_stds)
    A = mat_A(speed, yaw_rate, delta, self.x1[3])
    F = mat_F(speed, yaw_rate, delta, self.x1[3])

    # precision matrix of edge
    omega_edge = inv(A*M*A' + Matrix{Float64}(I,3,3).*0.0001)

    # precision matrix of graph
    self.omega_upper_left = F' * omega_edge * F
    self.omega_upper_right = -F' * omega_edge
    self.omega_lower_left = -omega_edge * F
    self.omega_lower_right = omega_edge

    # coefficient vector
    x2 = state_transition(speed, yaw_rate, delta, self.x1)
    self.xi_upper = F' * omega_edge * (self.x2 - x2)
    self.xi_lower = -omega_edge * (self.x2 - x2)

    return self
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using URIParser
using DataFrames
using NumericalDataManipulation.Data

abstract type Resource end
abstract type UriResource <: Resource end

struct FileResource <: UriResource
    uri::URI
    file_path::AbstractString

    function FileResource(a_string::AbstractString)
        given_string = strip(a_string)
        uri_string = if(ispath(given_string))
            "file://$(abspath(given_string))"
        elseif issomepath(given_string)
            if startswith(given_string, "file://")
                given_string
            else
                "file://$given_string"
            end
        end
        @assert(!isempty(uri_string), "The given file URI is an empty string.")
        uri = URI(uri_string)
        @assert(isvalid(uri), "The URI '$uri' is not valid.")
        @assert("file" == uri.scheme, "Only the 'file' scheme is supported; Given - '$(uri.scheme)'.")
        @assert(isempty(uri.host) ||
            "localhost" == uri.host ||
            "127.0.0.1" == uri.host, "Only localhost URIs are supported; Given - '$(uri)'.")
        uri_path = uri.path
        # TODO: convert uri_path to a platform specific file_path
        file_path = uri_path
        #@assert(isfile(file_path), "The path '$file_path' does not correspond to a file."
        this = new(uri, file_path)
        return this
    end
end
function Base.show(io::IO, resource::FileResource)
    compact = get(io, :compact, false)

    file_path = resource.file_path
    uri = resource.uri

    if compact
        print(io, "FileResource(path = $file_path, URI = $uri")
    else
        print(io, "FileResource(path = $file_path, URI = $uri")
    end

end
function Base.show(io::IO, ::MIME"text/plain", resource::FileResource)
    show(io, resource)
end

function exists(resource::FileResource)
    return isfile(resource.file_path)
end


text:
import DecFP
import ElectromagneticFields

using GeometricIntegrators.Equations

export guiding_center_4d_ode, guiding_center_4d_iode, guiding_center_4d_iode_λ,
       guiding_center_4d_iode_dec128, guiding_center_4d_vode,
       guiding_center_4d_dg, guiding_center_4d_formal_lagrangian

function guiding_center_4d_periodicity(q, periodic=true)
    period = zeros(eltype(q), size(q,1))

    if periodic
        period[1:3] .= periodicity(q)
    end

    return period
end

function guiding_center_4d_ode(qᵢ=qᵢ, params=parameters; periodic=true)
    ODE(guiding_center_4d_v, qᵢ;
            parameters=params, h=hamiltonian,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end

function guiding_center_4d_iode(qᵢ=qᵢ, params=parameters; periodic=true)
    pᵢ = guiding_center_4d_pᵢ(qᵢ)

    IODE(guiding_center_4d_ϑ, guiding_center_4d_f, guiding_center_4d_g, qᵢ, pᵢ;
            parameters=params, h=hamiltonian, v=guiding_center_4d_v,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end

function guiding_center_4d_iode_dec128(qᵢ=qᵢ, params=parameters; periodic=true)
    guiding_center_4d_iode(Dec128.(qᵢ), Dec128.(μ); periodic=periodic)
end

function guiding_center_4d_iode_λ(qᵢ=qᵢ, params=parameters; periodic=true)
    pᵢ = guiding_center_4d_pᵢ(qᵢ)
    λ₀ = guiding_center_4d_λ₀(qᵢ)

    IODE(guiding_center_4d_ϑ, guiding_center_4d_f, guiding_center_4d_g, qᵢ, pᵢ, λ₀;
            parameters=params, h=hamiltonian, v=guiding_center_4d_v,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end

function guiding_center_4d_vode(qᵢ=qᵢ, params=parameters; periodic=true)
    pᵢ = guiding_center_4d_pᵢ(qᵢ)

    VODE(guiding_center_4d_ϑ, guiding_center_4d_f, guiding_center_4d_g, qᵢ, pᵢ;
            parameters=params, h=hamiltonian, v=guiding_center_4d_v,
            Ω=guiding_center_4d_ω, ∇H=guiding_center_4d_dH,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end

function guiding_center_4d_dg(qᵢ=qᵢ, params=parameters; κ=0.0, periodic=true)
    guiding_center_4d_ϑ_κ(t, q, v, p, params) = guiding_center_4d_ϑ(t, q, v, p, params, κ)
    guiding_center_4d_f_κ(t, q, v, f, params) = guiding_center_4d_f(t, q, v, f, params, κ)
    guiding_center_4d_g_κ(t, q, λ, g, params) = guiding_center_4d_g(t, q, λ, g, params, κ)

    IODE(guiding_center_4d_ϑ_κ, guiding_center_4d_f_κ, guiding_center_4d_g_κ, qᵢ, qᵢ;
            parameters=params, h=hamiltonian, v=guiding_center_4d_v,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end

function guiding_center_4d_formal_lagrangian(qᵢ=qᵢ, params=parameters; periodic=true)
    pᵢ = guiding_center_4d_pᵢ(qᵢ)

    VODE(guiding_center_4d_ϑ, guiding_center_4d_f, guiding_center_4d_g, qᵢ, pᵢ;
            parameters=params, h=hamiltonian, v=guiding_center_4d_v,
            Ω=guiding_center_4d_ω, ∇H=guiding_center_4d_dH,
            periodicity=guiding_center_4d_periodicity(qᵢ, periodic))
end


text:
module TestPop

include("preamble.jl")

@testset "pop(sequence)" begin
    @test pop!!([0, 1]) == ([0], 1)
    @test pop!!((0, 1)) === ((0,), 1)
    @test pop!!((a=0, b=1)) === ((a=0,), 1)
    @test_broken pop!!(SVector(0, 1)) === (SVector(0), 1)
end

@testset "pop(assoc, key)" begin
    @test pop!!(Dict(:a => 1), :a) == (Dict(), 1)
    @test pop!!((a=1,), :a) == (NamedTuple(), 1)
end

@testset "pop(assoc, key, default)" begin
    @test pop!!(Dict(:a => 1), :a, nothing) == (Dict(), 1)
    @test pop!!(Dict(:a => 1), :b, nothing) == (Dict(:a => 1), nothing)
    @test pop!!((a=1,), :a, nothing) == (NamedTuple(), 1)
    @test pop!!((a=1,), :b, nothing) == ((a=1,), nothing)
end

@testset "mutation" begin
    @testset for args in [
        ([0],),
        (Dict(:a => 1), :a),
        (Dict(:a => 1), :a, nothing),
        (Dict(:a => 1), :b, nothing),
    ]
        @test pop!!(args...)[1] === args[1]
    end
end

end  # module


text:
"""
    calc_preds(X, Z, B)

Predict values

# Arguments

- X = 2d array of floats consisting of the row covariates, standardized as
  necessary
- Z = 2d array of floats consisting of the column covariates, standardized
  as necessary
- B = 2d array of floats consisting of coefficient estimates

# Value

2d array of floats

"""
function calc_preds(X::AbstractArray{Float64,2},
                    Z::AbstractArray{Float64,2},
                    B::AbstractArray{Float64,2})

    # Predict new values
    X*B*transpose(Z)
end

"""
    calc_preds(preds, X, Z, B)

Predict values in place

# Arguments

- preds = 2d array of floats consisting of the predicted values, to be
  updated in place
- X = 2d array of floats consisting of the row covariates, standardized as
  necessary
- Z = 2d array of floats consisting of the column covariates, standardized
  as necessary
- B = 2d array of floats consisting of coefficient estimates

# Value

None; updates predicted values in place.

"""
function calc_preds!(preds::AbstractArray{Float64,2},
                     X::AbstractArray{Float64,2},
                     Z::AbstractArray{Float64,2},
                     B::AbstractArray{Float64,2})

    # Update predicted values in place
    mul!(preds, X*B, transpose(Z))
end


text:
struct ManifoldKernelDensity{M <: MB.AbstractManifold, B <: BallTreeDensity, L, P}
  manifold::M
  """ legacy expects matrix of coordinates (as columns) """
  belief::B
  _partial::L
  """ just an example point for local access to the point data type"""
  _u0::P
  infoPerCoord::Vector{Float64}
end
const MKD{M,B,L} = ManifoldKernelDensity{M, B, L}


text:
export CUTENSORError

struct CUTENSORError <: Exception
    code::cutensorStatus_t
end

Base.convert(::Type{cutensorStatus_t}, err::CUTENSORError) = err.code

Base.showerror(io::IO, err::CUTENSORError) =
    print(io, "CUTENSORError: ", description(err), " (code $(reinterpret(Int32, err.code)), $(name(err)))")

name(err::CUTENSORError) = unsafe_string(cutensorGetErrorString(err))

## COV_EXCL_START
function description(err::CUTENSORError)
    if err.code == CUTENSOR_STATUS_SUCCESS
        "the operation completed successfully"
    elseif err.code == CUTENSOR_STATUS_NOT_INITIALIZED
        "the library was not initialized"
    elseif err.code == CUTENSOR_STATUS_ALLOC_FAILED
        "the resource allocation failed"
    elseif err.code == CUTENSOR_STATUS_INVALID_VALUE
        "an invalid value was used as an argument"
    elseif err.code == CUTENSOR_STATUS_ARCH_MISMATCH
        "an absent device architectural feature is required"
    elseif err.code == CUTENSOR_STATUS_MAPPING_ERROR
        "an access to GPU memory space failed"
    elseif err.code == CUTENSOR_STATUS_EXECUTION_FAILED
        "the GPU program failed to execute"
    elseif err.code == CUTENSOR_STATUS_INTERNAL_ERROR
        "an internal operation failed"
    elseif err.code == CUTENSOR_STATUS_NOT_SUPPORTED
        "operation not supported (yet)"
    elseif err.code == CUTENSOR_STATUS_LICENSE_ERROR
        "error detected trying to check the license"
    elseif err.code == CUTENSOR_STATUS_CUBLAS_ERROR
        "error occurred during a CUBLAS operation"
    elseif err.code == CUTENSOR_STATUS_CUDA_ERROR
        "error occurred during a CUDA operation"
    elseif err.code == CUTENSOR_STATUS_INSUFFICIENT_WORKSPACE
        "insufficient workspace memory for this operation"
    elseif err.code == CUTENSOR_STATUS_INSUFFICIENT_DRIVER
        "insufficient driver version"
    else
        "no description for this error"
    end
end
## COV_EXCL_STOP

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    throw(CUTENSORError(res))
end

function initialize_api()
    CUDA.prepare_cuda_call()
end

macro check(ex)
    quote
        res = @retry_reclaim isequal(CUTENSOR_STATUS_ALLOC_FAILED) $(esc(ex))
        if res != CUTENSOR_STATUS_SUCCESS
            throw_api_error(res)
        end

        return
    end
end


text:
module TestCounters

using Test
using ThreadLocalCounters

""" Get a tlc defined in this module by name. """
tlc_by_name(name) = only(
    c for c in ThreadLocalCounters.list(; all = true) if
    (c._module === @__MODULE__) && c.name === name
)

function test_tlc()
    tlc = tlc_by_name(:test_tlc)
    empty!(tlc)
    @test tlc ∉ ThreadLocalCounters.list()  # empty counter is not included

    nhits = 4 * Threads.nthreads()
    Threads.@threads for _ in 1:nhits
        @tlc test_tlc
    end
    @test tlc ∈ ThreadLocalCounters.list()

    @test sum(tlc) == nhits
    ThreadLocalCounters.clear()
    @test sum(tlc) == 0
end

end  # module


text:
struct Bounds
    lower::Array{Float64, 1}
    upper::Array{Float64, 1}
end

"""
    Bounds(dims)

Default constructor for a Bounds object. Returns a bounds object with a lower bound of [0...] and upper bound[1...]
with length of `dims`.
"""
Bounds(dims) = Bounds( zeros( dims ), ones( dims ) )

"""
    Bounds(dims)

Constructor for a Bounds object. Returns a bounds object with a lower bound of [lower...] and upper bound[upper...]
with length of `dims`.
"""
Bounds(lower, upper, dims) = Bounds( lower * ones( dims ), upper * ones( dims ) )

mutable struct Particle
    BestPos::Array{Float64, 1}
    BestScore::Float64
    Pos::Array{Float64, 1}
    Vel::Array{Float64, 1}
end

"""
    Particle(ProblemBounds, VelocityBounds)

Default constructor for a Particle object. It creates a random unformly distributed particle within the specified `ProblemBounds`,
and limits it's velocity to the specified `VelocityBounds`.
"""
function Particle(ProblemBounds, VelocityBounds)
    Dimen = length(ProblemBounds.upper)
    @assert Dimen == length(VelocityBounds.upper)
    BndDists = ProblemBounds.upper .- ProblemBounds.lower
    VelDists = VelocityBounds.upper .- VelocityBounds.lower
    PosI = (BndDists .* rand(Dimen)) .+ ProblemBounds.lower
    VelI = (VelDists .* rand(Dimen)) .+ VelocityBounds.lower
    return Particle( PosI, Inf, PosI, VelI )
end

#Particle( bounds(2,-1,3), bounds(0,1,3) )
"""
    PSO(fn, Bounds, VelRange, Particles; tolerance = 1e-6, maxiters = 1000, InertialDecay = 0.5, PersonalWeight = 0.5, GlobalWeight = 0.5, InternalParams = nothing)

Minimizes function `fn` with-in the user specified `Bounds` via a Particle Swarm Optimizer.
The particle velocities are limitted to the `VelRange`.
The number of particles are defined by the `Particles` parameter.

Returns a Tuple of the following form: ( GlobalBestPos, GlobalBestScore, P )
Where P is an array of the particles used in the optimization.

*Note: if the optimization function requires an additional constant parameter, please pass that parameter to InternalParams.
This will only work if the optimized parameter(o) and constant parameter(c) for the function of interest has the following format: F(o,c) *

Kennedy, J.; Eberhart, R. (1995). Particle Swarm Optimization. Proceedings of IEEE International Conference on Neural Networks. IV. pp. 1942–1948. doi:10.1109/ICNN.1995.488968
"""
function PSO(fn, Bounds::Bounds, VelRange::Bounds, Particles::Int;
                tolerance = 1e-6, maxiters = 1000,
                InertialDecay = 0.5, PersonalWeight = 0.5, GlobalWeight = 0.5, InternalParams = nothing)
    P = [ Particle(Bounds, VelRange) for x in 1 : Particles ]
    Dimen = length(Bounds.upper)
    GlobalBestPos = P[1].Pos
    GlobalBestScore = Inf
    iter = 0
    while (GlobalBestScore >= tolerance) && (iter < maxiters)
        for p in 1 : Particles
            if iter > 1
                #Update Positions
                (Wpersonal, Wglobal) = ( rand( Dimen ), rand( Dimen ) )
                P[p].Vel .=  ( InertialDecay .* P[p].Vel ) .+
                          ( Wpersonal .* PersonalWeight .* ( P[p].BestPos .-  P[p].Pos) ) .+
                          ( Wglobal   .* GlobalWeight   .* ( GlobalBestPos .- P[p].Pos ) )
                P[p].Pos .+= P[p].Vel
            end
            outbndslow = Bounds.lower .>= P[p].Pos
            if any(outbndslow)
                P[ p ].Pos[ outbndslow ] .= Bounds.lower[ outbndslow ]
            end
            outbndshigh = P[p].Pos .>= Bounds.upper
            if any(outbndshigh)
                P[ p ].Pos[ outbndshigh ] .= Bounds.upper[ outbndshigh ]
            end
            #Evaluate Fn to Obtain Score
            Score = isa(InternalParams, Nothing) ? fn(P[p].Pos) : fn(P[p].Pos, InternalParams)
            #Evaluate Scores
            if Score < P[p].BestScore
                P[p].BestPos .= P[p].Pos; P[p].BestScore = Score
                if Score < GlobalBestScore
                    GlobalBestScore = Score; GlobalBestPos = P[p].Pos
                end
            end
        end
        iter += 1
    end
    return ( GlobalBestPos, GlobalBestScore, P )
end

#parab(x) = sum(x.^2)

#PSO( parab, bounds(-3, 3, 1), bounds(0, 1, 1), 10;
#                tolerance = 1e-9, maxiters = 1000,
#                InertialDecay = 0.5, PersonalWeight = 0.5, GlobalWeight = 0.5)


text:
# joinpath(Pkg.dir("Brainstorm"), "src", "Math", "Arbs") |> cd

import Base: convert, one, show, zero

using Brainstorm._Math._MPFR: MPFR_RNDD, MPFR_RNDU
using Nemo
using Nemo: _arb_set, arf_struct

const Fslong = Int
const Fulong = UInt

const WORD_MAX = typemax(Fslong)
const WORD_MIN = typemin(Fslong)
const UWORD_MAX = typemax(Fulong)
const UWORD_MIN = typemin(Fulong)

const fmpr_rnd_t = Cint

const FMPR_RND_DOWN = 0
const FMPR_RND_UP = 1
const FMPR_RND_FLOOR = 2
const FMPR_RND_CEIL = 3
const FMPR_RND_NEAR = 4

const arb_t = arb

include("arf.jl")
include("arf_interval.jl")
include("arb.jl")
include("arb_calc.jl")

#ArbRoots
#- https://github.com/fredrik-johansson/arb/blob/master/examples/real_roots.c
#- https://wiki.sagemath.org/days4schedule?action=AttachFile&do=get&target=witty.pdf

func(x, o) = begin
    α = x^995(x^2 - 9999)^2 - 1
    β = o > 1 ? 4(x^996)(-9999 + x^2) + 995(x^994)(-9999 + x^2)^2 : zero(x)
    α, β
end

arbParent = ArbField(8192)

interv = (big(98.0), big(101.0))
maxdepth = 100
maxeval = 100000
maxfound = 100
iter = 10000

α = arb_calc_isolate_roots(func, interv, maxdepth, maxeval, maxfound, arbParent, true)
β = [arb_calc_refine_root_bisect(func, x, iter, arbParent) for x in α[1]]
γ = [(func(x[1], 1)[1], func(x[2], 1)[1]) for (_, x) in β]


text:
"""
    Read an XDA Mesh generated by libMesh
"""
function readXDAMesh(filename::String)

    startLog(main_perf_log, "readXDAMesh()")

    mesh = Mesh()

    # Set the sideset / nodeset numberings:
    boundary_info = mesh.boundary_info

    open(filename) do f
        lines = readlines(filename)

        lines = [split(line) for line in lines]

        current_line = 1

        # First line is a version number
        current_line += 1

        num_elem = parse(Int64, lines[current_line][1])
        current_line += 1

#        println("num_elem: ", num_elem)

        num_nodes = parse(Int64, lines[current_line][1])
        current_line += 1

#        println("num_nodes: ", num_elem)

        # BC file: 4
        current_line += 1
        # subdomain file: 5
        current_line += 1
        # proc file: 6
        current_line += 1
        # p-level file: 7
        current_line += 1
        # type size: 8
        current_line += 1
        # uid size: 9
        current_line += 1
        # pid size: 10
        current_line += 1
        # sid size: 11
        current_line += 1
        # p-level size: 12
        current_line += 1
        # eid size: 13
        current_line += 1
        # side size: 14
        current_line += 1
        # bid size: 15
        current_line += 1
        # subdomain to name map: 16
        current_line += 1

        num_elem_at_level_0 = parse(Int64, lines[current_line][1])
        current_line += 1

#        println("num_elem_at_level_0: ", num_elem_at_level_0)

        # The next lines are elements, but we don't want to read those yet
        elem_lines_start = current_line
        current_line += num_elem_at_level_0

        nodes = Array{Node}(num_nodes)

        # Let's read the nodes
        for node_id in 1:num_nodes
            x = parse(Float64, lines[current_line][1])
            y = parse(Float64, lines[current_line][2])

            node = Node{2}(node_id, Vec{2}((x,y)), [], invalid_processor_id)

            nodes[node_id] = node

            current_line += 1
        end

#        println(nodes)

        # NOW let's read the elements
        elements = Array{Element}(num_elem)

        for elem_id in 1:num_elem_at_level_0
            # Parse the line it is stored like this:
            # [ type sid (n0 ... nN-1) ]
            # All we need are the node IDs (but we need to transfer them to 1-based indexing
            node_ids = [parse(Int64, node_id) + 1 for node_id in lines[elem_lines_start + elem_id - 1][3:end]]

            element = Element(elem_id,
                              [nodes[node_id] for node_id in node_ids],
                              [],
                              invalid_processor_id)

            elements[elem_id] = element
        end

#        println(elements)

        # Unique ID presence
        current_line += 1
        # Sideset id to name map
        has_sideset_names = parse(Int64, lines[current_line][1])
        current_line += 1

        if has_sideset_names > 0
            # Vector length of sideset IDs with names
            current_line += 1

            # Sideset IDs: translate to 1-based
            sideset_ids = [parse(Int64, id) + 1 for id in lines[current_line]]
            current_line += 1

#            println("sideset_ids: ", sideset_ids)

            for bid in sideset_ids
                boundary_info.side_list[bid] = Array{ElemSidePair}(0)
            end

            # Vector length of sideset names
            current_line += 1

            # Sideset names (not currently used)
            sideset_names = lines[current_line]
            current_line += 1
        end

        num_side_bcs = parse(Int64, lines[current_line][1])
        current_line += 1

        # Read in the element/side pairs
        for i in 1:num_side_bcs
            data = [parse(Int64, val) + 1 for val in lines[current_line]]

            if !(data[3] in keys(boundary_info.side_list))
                boundary_info.side_list[data[3]] = Array{ElemSidePair}(0)
            end

            push!(boundary_info.side_list[data[3]], ElemSidePair(elements[data[1]], data[2]))

            current_line += 1
        end

        # Nodest id to name map
        has_nodeset_names = parse(Int64, lines[current_line][1])
        current_line += 1

        if has_nodeset_names > 0
            # Num nodesets with names
            current_line += 1

            nodeset_ids = [parse(Int64, id) + 1 for id in lines[current_line]]
            current_line += 1

#            println("nodeset_ids: ", nodeset_ids)

            for bid in nodeset_ids
                boundary_info.node_list[bid] = Array{Node}(0)
            end

            # Length of names
            current_line += 1

            # nodeset names (not currently used)
            nodeset_names = lines[current_line]
            current_line += 1
        end

        num_nodesets = parse(Int64, lines[current_line][1])
        current_line += 1

        for i in 1:num_nodesets
            data = [parse(Int64, val) + 1 for val in lines[current_line]]

            if !(data[2] in keys(boundary_info.node_list))
                boundary_info.node_list[data[2]] = Array{Node}(0)
            end

            push!(boundary_info.node_list[data[2]], nodes[data[1]])

            current_line += 1
        end

#        println("node_list: ", boundary_info.node_list)

        mesh.elements = elements
        mesh.nodes = nodes

        initialize!(mesh)
    end

    stopLog(main_perf_log, "readXDAMesh()")

    return mesh
end


text:
#manytoone_Array{Int}
function my_deferred_acceptance(prop_prefs::Array{Int},
                                resp_prefs::Array{Int},
                                caps::Vector{Int})
    prop_size = size(prop_prefs, 2)
    resp_size = size(resp_prefs, 2)
    prop_matched = zeros(Int64, prop_size)
    resp_matched = zeros(Int64, sum(caps))
    next_prop = zeros(Int64, prop_size)
    max_prop = Int64[]

    for p in 1:prop_size
        push!(max_prop, find(prop_prefs[:, p] .== 0)[1]-1)
    end

    indptr = Array{Int}(resp_size+1)
    indptr[1] = 1
    for i in 1:resp_size
        indptr[i+1] = indptr[i] + caps[i]
    end

    while any(prop_matched .== 0) == true
        prop_single = find(prop_matched .== 0)
        if all(next_prop[prop_single] .>= max_prop[prop_single]) == true
            break
        else
            for each_prop_single in prop_single
                proposing = prop_prefs[next_prop[each_prop_single]+1, each_prop_single]
                if proposing != 0
                    next_prop[each_prop_single] = next_prop[each_prop_single]+1
                    if sum(resp_matched[indptr[proposing]:indptr[proposing+1]-1] .== 0) != 0 && (find(resp_prefs[:, proposing] .==each_prop_single) .<  find(resp_prefs[:, proposing] .==0)) == [true]
                        prop_matched[each_prop_single] = proposing
                        matched_index = findfirst(resp_matched[indptr[proposing]:indptr[proposing+1]-1] .== 0)
                        resp_matched[matched_index + indptr[proposing] - 1] = each_prop_single
                    elseif sum(resp_matched[indptr[proposing]:indptr[proposing+1]-1] .== 0) .== 0
                        current_order = Int64[]
                        for i in 1:caps[proposing] push!(current_order, find(resp_prefs[:, proposing] .== resp_matched[indptr[proposing]+i-1])[1]) end
                        if all(findmax(current_order)[1] .> find(resp_prefs[:, proposing] .== each_prop_single)) == true
                        prop_matched[each_prop_single] = proposing
                        prop_matched[resp_matched[findmax(current_order)[2] + indptr[proposing] - 1]] = 0
                        resp_matched[findmax(current_order)[2] + indptr[proposing] - 1] = each_prop_single
                        end
                    end
                end
            end
        end
    end

    return prop_matched, resp_matched, indptr
end
#manytoone_Vector{Vector{Int}}
function my_deferred_acceptance(prop_prefs::Vector{Vector{Int}},
                                resp_prefs::Vector{Vector{Int}},
                                caps::Vector{Int})
    prop_prefs_2d = Array{Int64}(length(resp_prefs)+1, length(prop_prefs))
    for i in 1:length(prop_prefs)
        if length(prop_prefs[i]) != length(resp_prefs)
            x = vcat(prop_prefs[i], 0)
            x = vcat(x, setdiff([i for i in 1:length(resp_prefs)], prop_prefs[i]))
            prop_prefs_2d[:,i] = x
        else
            prop_prefs_2d[:,i] = vcat(prop_prefs[i], 0)
        end
    end
    resp_prefs_2d = Array{Int64}(length(prop_prefs)+1, length(resp_prefs))
    for i in 1:length(resp_prefs)
        if length(resp_prefs[i]) != length(prop_prefs)
            x = vcat(resp_prefs[i], 0)
            x = vcat(x, setdiff([i for i in 1:length(prop_prefs)], resp_prefs[i]))
            resp_prefs_2d[:,i] = x
        else
            resp_prefs_2d[:,i] = vcat(resp_prefs[i], 0)
        end
    end
    return my_deferred_acceptance(prop_prefs_2d, resp_prefs_2d, caps)
end
#onetoone_Array{Int}
function my_deferred_acceptance(prop_prefs::Array{Int},
                                resp_prefs::Array{Int})
    caps = ones(Int, size(resp_prefs, 2))
    prop_matches, resp_matches, indptr =
        my_deferred_acceptance(prop_prefs, resp_prefs, caps)
    return prop_matches, resp_matches
end
#onetoone_Vector{Vector{Int}}
function my_deferred_acceptance(prop_prefs::Vector{Vector{Int}},
                                resp_prefs::Vector{Vector{Int}})
    prop_prefs_2d = Array{Int64}(length(resp_prefs)+1, length(prop_prefs))
    for i in 1:length(prop_prefs)
        if length(prop_prefs[i]) != length(resp_prefs)
            x = vcat(prop_prefs[i], 0)
            x = vcat(x, setdiff([i for i in 1:length(resp_prefs)], prop_prefs[i]))
            prop_prefs_2d[:,i] = x
        else
            prop_prefs_2d[:,i] = vcat(prop_prefs[i], 0)
        end
    end
    resp_prefs_2d = Array{Int64}(length(prop_prefs)+1, length(resp_prefs))
    for i in 1:length(resp_prefs)
        if length(resp_prefs[i]) != length(prop_prefs)
            x = vcat(resp_prefs[i], 0)
            x = vcat(x, setdiff([i for i in 1:length(prop_prefs)], resp_prefs[i]))
            resp_prefs_2d[:,i] = x
        else
            resp_prefs_2d[:,i] = vcat(resp_prefs[i], 0)
        end
    end
    caps = ones(Int, size(resp_prefs_2d, 2))
    prop_matches, resp_matches, indptr =
        my_deferred_acceptance(prop_prefs_2d, resp_prefs_2d, caps)
    return prop_matches, resp_matches
end


text:
"""Parallel reduction: put max at u[1,1]"""
function reduceMax!(u)
    tx = blockDim().x * (blockIdx().x - 1) + threadIdx().x;
    ty = blockDim().y * (blockIdx().y - 1) + threadIdx().y;

    if tx < size(u, 1) + 1 && ty < size(u, 2) + 1
        # reduce over x
        stride = blockDim().x >> 1
        while stride > 0
            sync_threads()
            if threadIdx().x <= stride && tx + stride <= size(u, 1)
                if u[tx,ty] < u[tx + stride,ty]
                    u[tx,ty] = u[tx + stride,ty]
                end
            end
            stride >>= 1
        end

        # reduce over y
        stride = blockDim().y >> 1
        while stride > 0
            sync_threads()
            if threadIdx().y <= stride && ty + stride <= size(u, 2)
                if u[tx,ty] < u[tx,ty + stride]
                    u[tx,ty] = u[tx,ty + stride]
                end
            end
            stride >>= 1
        end
    end

    return nothing
end

# for develop, todo: how to unit test this stuff?
function testReduceMax()
    N = 160
    u = rand(N, N)
    u_d = CuArray(u)

    ths = (16, 16)
    bls = (Int(ceil(N / ths[1])), Int(ceil(N / ths[2])))
    @cuda threads = ths blocks = bls reduceMax!(u_d)
    @test u_d[1,1] == maximum(u)
end


text:
using Test, Random, Statistics, Zygote, SpecialFunctions
using TruncatedNormal: randnt, randnt_half, ∇randnt_half

Random.seed!(20)

@inferred randnt_half(1.0, 2.0)
@inferred randnt_half(Float32(1.0), Float32(2.0))
@test randnt_half(Float32(1.0), Float32(2.0)) isa Float32

# compare exact 1st and 2nd moments to Monte Carlo estimates
m1(μ,σ) = μ + σ * √(2/π) / erfcx(-μ/σ/√2)
m2(μ,σ) = μ^2 + σ^2 + μ * σ * √(2/π) / erfcx(-μ/σ/√2)

for μ = -1:1, σ = 1:2
    samples = [randnt_half(μ,σ) for _ = 1:10^6]
    @test mean(samples.^1) ≈ m1(μ,σ) atol=1e-2
    @test mean(samples.^2) ≈ m2(μ,σ) atol=1e-2

    dμ, dσ = gradient(m1, μ, σ)
    dμmc, dσmc = gradient(μ,σ) do μ,σ
        mean([randnt_half(μ,σ) for _ = 1:10^4])
    end
    @test dμmc ≈ dμ atol=0.1
    @test dσmc ≈ dσ atol=0.1

    dμ, dσ = gradient(μ,σ) do μ,σ
        sqrt(m2(μ,σ) - m1(μ,σ)^2)
    end
    dμmc, dσmc = gradient(μ,σ) do μ,σ
        samples = [randnt_half(μ,σ) for _ = 1:10^4]
        std(samples)
    end
    @test dμmc ≈ dμ atol=0.1
    @test dσmc ≈ dσ atol=0.1
end

# broadcasted versions
μ = 3randn(2,2)
σ = 3rand(2,2)
dμ, dσ = gradient(μ,σ) do μ,σ
    mean(m1.(μ,σ))
end
dμmc, dσmc = gradient(μ,σ) do μ,σ
    m = zero(randnt_half.(μ,σ))
    for _ = 1:10^4
        m += randnt_half.(μ,σ)
    end
    mean(m / 10^4)
end
@test_broken dμmc ≈ dμ atol=0.1
@test_broken dσmc ≈ dσ atol=0.1

μ = 3randn(2,2)
σ = 3rand(2,2)
dμ, dσ = gradient(μ,σ) do μ,σ
    mean(@. m2(μ,σ))
end
dμmc, dσmc = gradient(μ,σ) do μ,σ
    m = zero(randnt_half.(μ,σ))
    for _ = 1:10^4
        m += @. randnt_half(μ, σ)^2
    end
    mean(m ./ 10^4)
end
@test_broken dμmc ≈ dμ atol=0.1
@test_broken dσmc ≈ dσ atol=0.1
# https://github.com/FluxML/Zygote.jl/issues/1088
# My rule is not being called when broadcasting because Zygote goes forward-mode.


text:
function position(trigger::Maple.Trigger)::Tuple{Int, Int}
    return floor(Int, trigger.x), floor(Int, trigger.y)
end

function renderTrigger(ctx::Cairo.CairoContext, layer::Layer, trigger::Maple.Trigger, room::Maple.Room; alpha=nothing)
    if ctx.ptr != C_NULL
        Cairo.save(ctx)

        x, y = position(trigger)
        w, h = Int(trigger.width), Int(trigger.height)

        rectangle(ctx, x, y, w, h)
        clip(ctx)

        text = humanizeVariableName(trigger.name)
        drawRectangle(ctx, x, y, w, h, colors.trigger_fc, colors.trigger_bc)
        drawCenteredText(ctx, text, x, y, w, h)

        restore(ctx)
    end
end

function renderTriggerSelection(ctx::Cairo.CairoContext, layer::Layer, trigger::Maple.Trigger, room::Maple.Room; alpha=nothing)
    x, y = Int(trigger.data["x"]), Int(trigger.data["y"])
    width, height = Int(trigger.data["width"]), Int(trigger.data["height"])
    nodes = get(trigger.data, "nodes", Tuple{Int, Int}[])
    offsetCenterX, offsetCenterY = floor(Int, width / 2), floor(Int, height / 2)

    text = humanizeVariableName(trigger.name)

    for node in nodes
        nx, ny = Int.(node)

        drawArrow(ctx, x + offsetCenterX, y + offsetCenterY, nx + 4, ny + 4, colors.selection_selected_fc, headLength=6)
        drawRectangle(ctx, nx, ny, 8, 8, colors.trigger_fc, colors.trigger_bc)
    end
end

nodeLimits(trigger::Maple.Trigger) = 0, 0

editingOptions(trigger::Maple.Trigger) = Dict{String, Any}()

minimumSize(trigger::Maple.Trigger) = 8, 8
resizable(trigger::Maple.Trigger) = true, true

function triggerSelection(trigger::Maple.Trigger, room::Maple.Room, node::Int=0)
    x, y = Int(trigger.x), Int(trigger.y)
    width, height = Int(trigger.width), Int(trigger.height)
    nodes = get(trigger.data, "nodes", Tuple{Int, Int}[])

    if isempty(nodes)
        return Rectangle(x, y, width, height)

    else
        res = Rectangle[Rectangle(x, y, width, height)]

        for node in nodes
            nx, ny = Int.(node)

            push!(res, Rectangle(nx, ny, 8, 8))
        end

        return res
    end
end

const loadedTriggers = joinpath.(abs"triggers", readdir(abs"triggers"))
loadModule.(loadedTriggers)
const triggerPlacements = PlacementDict()
const triggerPlacementsCache = Dict{String, Trigger}()

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
filename() = "day09.txt"

read_data() = parse.(Int, readlines(filename()))

# read_data() = parse.(Int, split("""35
# 20
# 15
# 25
# 47
# 40
# 62
# 55
# 65
# 95
# 102
# 117
# 150
# 182
# 127
# 219
# 299
# 277
# 309
# 576""", "\n"))

function part1(pre)
    v = read_data()
    for i in pre+1:length(v)
        twosum(v[i-pre:i-1], v[i]) || return v[i]
    end
end

function twosum(nums, target)
    # @info "Checking $target: $nums"
    set = Set{Int}()
    for v in nums
        complement = target - v
        complement in set && return true
        push!(set, v)
    end
    return false
end

# part 1
# part1(25) # 41682220

# part 2

function part2_brute_force(v, target)
    for i in 1:length(v)
        for j in i:length(v)
            sum(v[i:j]) == target && return i, j, sum(extrema(v[i:j]))
        end
    end
    return -1
end

function part2(v, target)
    i = 1
    j = 2
    s = v[i] + v[j]
    while j < length(v)
        if s < target
            j += 1
            s += v[j]
        elseif s > target
            s -= v[i]
            i += 1
        else
            return i, j , sum(extrema(v[i:j]))
        end
    end
    return -1
end

function part2()
    v = read_data()
    part2(v, 41682220)
end

#=
Brute force is very slow.

julia> @time part2(read_data(), 41682220)
  0.000209 seconds (2.02 k allocations: 72.328 KiB)
(422, 438, 5388976)

julia> @time part2_brute_force(read_data(), 41682220)
  0.181522 seconds (334.63 k allocations: 1.049 GiB, 33.86% gc time)
(422, 438, 5388976)
=#


text:
module TestIndex

using Test, DataFrames
using DataFrames: Index, SubIndex, fuzzymatch

@testset "Index indexing" begin
    i = Index()
    push!(i, :A)
    push!(i, :B)

    inds = Any[1,
               big(1),
               :A,
               [true, false],
               [1],
               [big(1)],
               big(1):big(1),
               [:A],
               Union{Int, Missing}[1],
               Union{BigInt, Missing}[big(1)],
               Union{Symbol, Missing}[:A],
               Any[1],
               Any[:A]]

    for ind in inds
        if ind == :A || ndims(ind) == 0
            @test i[ind] == 1
        else
            @test (i[ind] == [1])
        end
    end

    @test_throws MethodError i[1.0]
    @test_throws ArgumentError i[true]
    @test_throws ArgumentError i[false]
    @test_throws ArgumentError i[Union{Bool, Missing}[true, false]]
    @test_throws ArgumentError i[Any[1, missing]]
    @test_throws ArgumentError i[[1, missing]]
    @test_throws ArgumentError i[[true, missing]]
    @test_throws ArgumentError i[Any[true, missing]]
    @test_throws MethodError i[[:A, missing]]
    @test_throws MethodError i[Any[:A, missing]]
    @test_throws ArgumentError i[1.0:1.0]
    @test_throws ArgumentError i[[1.0]]
    @test_throws ArgumentError i[Any[1.0]]
    @test_throws BoundsError i[0]
    @test_throws BoundsError i[10]
    @test_throws ArgumentError i[:x]
    @test_throws BoundsError i[1:3]
    @test_throws ArgumentError i[[1,1]]
    @test_throws ArgumentError i[[:A,:A]]
    @test_throws BoundsError i[Not(0)]
    @test_throws BoundsError i[Not(10)]
    @test_throws ArgumentError i[Not(:x)]
    @test_throws BoundsError i[Not(1:3)]
    @test_throws ArgumentError i[Not([1,1])]
    @test_throws ArgumentError i[Not([:A,:A])]

    @test i[1:1] == 1:1

    @test_throws BoundsError i[[true]]
    @test_throws BoundsError i[true:true]
    @test_throws BoundsError i[[true, false, true]]

    @test_throws ArgumentError i[["a"]]
    @test_throws ArgumentError i[Any["a"]]

    @test i[[]] == Int[]
    @test i[Int[]] == Int[]
    @test i[Symbol[]] == Int[]
    @test i[:] == 1:length(i)
    @test i[Not(Not(:))] == 1:length(i)
    @test i[Not(1:0)] == 1:length(i)
end

@testset "rename!" begin
    i = Index([:A,:B])
    @test names(i) == [:A,:B]
    @test rename!(i, [:a,:a], makeunique=true) == Index([:a,:a_1])
    @test_throws ArgumentError rename!(i, [:a,:a])
    @test_throws DimensionMismatch rename!(i, [:a,:b,:c])
    @test rename!(copy(i), [:a,:b]) == Index([:a,:b])
    @test names(i) == [:a,:a_1]
    @test rename!(i, [:a,:b]) == Index([:a,:b])
    @test rename!(copy(i), [:a => :A]) == Index([:A,:b])
    @test rename!(copy(i), [:a => :a]) == Index([:a,:b])
    @test rename!(copy(i), [:a => :b, :b => :a]) == Index([:b,:a])
    @test rename!(x -> Symbol(uppercase(string(x))), copy(i)) == Index([:A,:B])
    @test rename!(x -> Symbol(lowercase(string(x))), copy(i)) == Index([:a,:b])

    @test delete!(i, :a) == Index([:b])
    push!(i, :C)
    @test delete!(i, 1) == Index([:C])

    i = Index([:A, :B, :C, :D, :E])
    i2 = copy(i)
    rename!(i2, reverse(names(i2)))
    rename!(i2, reverse(names(i2)))
    @test names(i2) == names(i)
    for name in names(i)
        i2[name] # Issue #715
    end
end

@testset "SubIndex" begin
    i = Index([:A, :B, :C, :D, :E])
    si1 = SubIndex(i, :)
    si2 = SubIndex(i, 3:5)
    si3 = SubIndex(i, [3,4,5])
    si4 = SubIndex(i, [false, false, true, true, true])
    si5 = SubIndex(i, [:C, :D, :E])
    si6 = SubIndex(i, Not(Not([:C, :D, :E])))
    si7 = SubIndex(i, Not(1:2))

    @test copy(si1) == i
    @test copy(si2) == Index([:C, :D, :E])
    @test copy(si3) == Index([:C, :D, :E])
    @test copy(si4) == Index([:C, :D, :E])
    @test copy(si5) == Index([:C, :D, :E])
    @test copy(si6) == Index([:C, :D, :E])
    @test copy(si7) == Index([:C, :D, :E])

    @test_throws ArgumentError SubIndex(i, 1)
    @test_throws ArgumentError SubIndex(i, :A)
    @test_throws ArgumentError SubIndex(i, true)
    @test si1 isa Index
    @test si2.cols == 3:5
    @test si2.remap == -1:3
    @test si3.cols == 3:5
    @test si3.remap == [0, 0, 1, 2, 3]
    @test !haskey(si3, :A)
    @test si3.remap == [0, 0, 1, 2, 3]
    @test si4.cols == 3:5
    @test si4.remap == [0, 0, 1, 2, 3]
    @test !haskey(si4, :A)
    @test si4.remap == [0, 0, 1, 2, 3]
    @test si5.cols == 3:5
    @test si5.remap == [0, 0, 1, 2, 3]
    @test !haskey(si5, :A)
    @test si5.remap == [0, 0, 1, 2, 3]
    @test si6.cols == 3:5
    @test si6.remap == [0, 0, 1, 2, 3]
    @test !haskey(si6, :A)
    @test si6.remap == [0, 0, 1, 2, 3]
    @test si7.cols == 3:5
    @test si7.remap == [0, 0, 1, 2, 3]
    @test !haskey(si7, :A)
    @test si7.remap == [0, 0, 1, 2, 3]

    @test length(si1) == 5
    @test length(si2) == 3
    @test length(si3) == 3
    @test length(si4) == 3
    @test length(si5) == 3
    @test length(si6) == 3
    @test length(si7) == 3

    @test names(si1) == keys(si1) == [:A, :B, :C, :D, :E]
    @test names(si2) == keys(si2) == [:C, :D, :E]
    @test names(si3) == keys(si3) == [:C, :D, :E]
    @test names(si4) == keys(si4) == [:C, :D, :E]
    @test names(si5) == keys(si5) == [:C, :D, :E]
    @test names(si6) == keys(si5) == [:C, :D, :E]
    @test names(si7) == keys(si5) == [:C, :D, :E]

    @test_throws ArgumentError haskey(si3, true)
    @test haskey(si3, 1)
    @test !haskey(si3, 0)
    @test !haskey(si3, 4)
    @test haskey(si3, :D)
    @test !haskey(si3, :A)
    @test si3[:C] == 1
    @test si3[names(i)] == [0, 0, 1, 2, 3]
end

@testset "selector mutation" begin
    df = DataFrame(a=1:5, b=11:15, c=21:25)
    selector1 = [3,2]
    dfv1 = view(df, :, selector1)
    dfr1 = view(df, 2, selector1)
    selector2 = [1]
    dfv2 = view(dfv1, :, selector2)
    dfr2 = view(dfr1, selector2)
    @test names(dfv1) == [:c, :b]
    @test names(dfv2) == [:c]
    @test names(dfr1) == [:c, :b]
    @test names(dfr2) == [:c]
    selector1[1] = 1
    @test names(dfv1) == [:a, :b]
    @test names(dfv2) == [:c]
    @test names(dfr1) == [:a, :b]
    @test names(dfr2) == [:c]
    selector3 = [:c, :b]
    dfv3 = view(df, :, selector3)
    dfr3 = view(df, 2, selector3)
    @test names(dfv3) == [:c, :b]
    @test names(dfr3) == [:c, :b]
    selector3[1] = :a
    @test names(dfv3) == [:c, :b]
    @test names(dfr3) == [:c, :b]
end

@testset "fuzzy matching" begin
    i = Index()
    push!(i, :x1)
    push!(i, :x12)
    push!(i, :x131)
    push!(i, :y13)
    push!(i, :yy13)
    @test_throws ArgumentError i[:x13]
    @test_throws ArgumentError i[:xx13]
    @test_throws ArgumentError i[:yy14]
    @test_throws ArgumentError i[:abcd]
    @test fuzzymatch(i.lookup, :x13) == [:x1, :x12, :x131, :y13, :yy13]
    @test fuzzymatch(i.lookup, :xx1314) == [:x131]
    @test fuzzymatch(i.lookup, :yy14) == [:yy13, :y13]
    @test isempty(fuzzymatch(i.lookup, :abcd))
end

@testset "Regex indexing" begin
    i = Index()
    push!(i, :x1)
    push!(i, :x12)
    push!(i, :x131)
    push!(i, :y13)
    push!(i, :yy13)
    @test i[r"x1."] == [2, 3]
    @test isempty(i[r"xx"])
    @test i[r""] == 1:5
    @test names(SubIndex(i, r"x1.")) == [:x12, :x131]
    @test isempty(names(SubIndex(i, r"xx")))
    @test names(SubIndex(i, r"")) == names(i)
    @test DataFrames._names(SubIndex(i, r"x1.")) == [:x12, :x131]
    @test isempty(DataFrames._names(SubIndex(i, r"xx")))
    @test DataFrames._names(SubIndex(i, r"")) == names(i)
    @test DataFrames.parentcols(SubIndex(i, r"x1.")) == [2, 3]
    @test isempty(DataFrames.parentcols(SubIndex(i, r"xx")))
    @test DataFrames.parentcols(SubIndex(i, r"")) == 1:5
    @test DataFrames.parentcols(SubIndex(i, All())) == 1:5
    @test DataFrames.parentcols(SubIndex(i, Between(:x1, :x12))) == 1:2

    i2 = SubIndex(i, r"")
    @test i2[r"x1."] == [2, 3]
    @test isempty(i2[r"xx"])
    @test i2[r""] == 1:5
    @test names(SubIndex(i2, r"x1.")) == [:x12, :x131]
    @test isempty(names(SubIndex(i2, r"xx")))
    @test names(SubIndex(i2, r"")) == names(i)
    @test DataFrames._names(SubIndex(i2, r"x1.")) == [:x12, :x131]
    @test isempty(DataFrames._names(SubIndex(i2, r"xx")))
    @test DataFrames._names(SubIndex(i2, r"")) == names(i2)
    @test DataFrames.parentcols(SubIndex(i2, r"x1.")) == [2, 3]
    @test isempty(DataFrames.parentcols(SubIndex(i2, r"xx")))
    @test DataFrames.parentcols(SubIndex(i2, r"")) == 1:5
    @test DataFrames.parentcols(SubIndex(i2, All())) == 1:5
    @test DataFrames.parentcols(SubIndex(i2, Between(:x1, :x12))) == 1:2

    i3 = SubIndex(i, r"x1.")
    @test i3[r"x1.$"] == [1]
    @test isempty(i3[r"xx"])
    @test i3[r""] == 1:2
    @test names(SubIndex(i3, r"x1.$")) == [:x12]
    @test isempty(names(SubIndex(i3, r"xx")))
    @test names(SubIndex(i3, r"")) == names(i3)
    @test DataFrames._names(SubIndex(i3, r"x1.$")) == [:x12]
    @test isempty(DataFrames._names(SubIndex(i3, r"xx")))
    @test DataFrames._names(SubIndex(i3, r"")) == names(i3)
    @test DataFrames.parentcols(SubIndex(i3, r"x1.$")) == [1]
    @test isempty(DataFrames.parentcols(SubIndex(i3, r"xx")))
    @test DataFrames.parentcols(SubIndex(i3, r"")) == 1:2
    @test DataFrames.parentcols(SubIndex(i3, All())) == 1:2
    @test_throws BoundsError DataFrames.parentcols(SubIndex(i3, Between(:x1, :x12))) == 1:2
    @test DataFrames.parentcols(SubIndex(i3, Between(:x12, :x12))) == 1:1
end

@testset "Not indexing" begin
    i = Index()
    push!(i, :x1)
    push!(i, :x12)
    push!(i, :x131)
    push!(i, :y13)
    push!(i, :yy13)
    @test i[Not(Not(r"x1."))] == [2, 3]
    @test isempty(i[Not(Not(r"xx"))])
    @test i[Not(Not(r""))] == 1:5
    @test names(SubIndex(i, Not(Not(r"x1.")))) == [:x12, :x131]
    @test isempty(names(SubIndex(i, Not(Not(r"xx")))))
    @test names(SubIndex(i, Not(Not(r"")))) == names(i)
    @test DataFrames._names(SubIndex(i, Not(Not(r"x1.")))) == [:x12, :x131]
    @test isempty(DataFrames._names(SubIndex(i, Not(Not(r"xx")))))
    @test DataFrames._names(SubIndex(i, Not(Not(r"")))) == names(i)
    @test DataFrames.parentcols(SubIndex(i, Not(Not(r"x1.")))) == [2, 3]
    @test isempty(DataFrames.parentcols(SubIndex(i, Not(Not(r"xx")))))
    @test DataFrames.parentcols(SubIndex(i, Not(Not(r"")))) == 1:5

    i2 = SubIndex(i, Not(Not(r"")))
    @test i2[Not(Not(r"x1."))] == [2, 3]
    @test isempty(i2[Not(Not(r"xx"))])
    @test i2[Not(Not(r""))] == 1:5
    @test names(SubIndex(i2, Not(Not(r"x1.")))) == [:x12, :x131]
    @test isempty(names(SubIndex(i2, Not(Not(r"xx")))))
    @test names(SubIndex(i2, Not(Not(r"")))) == names(i)
    @test DataFrames._names(SubIndex(i2, Not(Not(r"x1.")))) == [:x12, :x131]
    @test isempty(DataFrames._names(SubIndex(i2, Not(Not(r"xx")))))
    @test DataFrames._names(SubIndex(i2, Not(Not(r"")))) == names(i2)
    @test DataFrames.parentcols(SubIndex(i2, Not(Not(r"x1.")))) == [2, 3]
    @test isempty(DataFrames.parentcols(SubIndex(i2, Not(Not(r"xx")))))
    @test DataFrames.parentcols(SubIndex(i2, Not(Not(r"")))) == 1:5

    i3 = SubIndex(i, Not(Not(r"x1.")))
    @test i3[Not(Not(r"x1.$"))] == [1]
    @test isempty(i3[Not(Not(r"xx"))])
    @test i3[Not(Not(r""))] == 1:2
    @test names(SubIndex(i3, Not(Not(r"x1.$")))) == [:x12]
    @test isempty(names(SubIndex(i3, Not(Not(r"xx")))))
    @test names(SubIndex(i3, Not(Not(r"")))) == names(i3)
    @test DataFrames._names(SubIndex(i3, Not(Not(r"x1.$")))) == [:x12]
    @test isempty(DataFrames._names(SubIndex(i3, Not(Not(r"xx")))))
    @test DataFrames._names(SubIndex(i3, Not(Not(r"")))) == names(i3)
    @test DataFrames.parentcols(SubIndex(i3, Not(Not(r"x1.$")))) == [1]
    @test isempty(DataFrames.parentcols(SubIndex(i3, Not(Not(r"xx")))))
    @test DataFrames.parentcols(SubIndex(i3, Not(Not(r"")))) == 1:2
end

@testset "Between indexing" begin
    df = DataFrame(a=1, b=2, c=3)
    @test select(df, Between(1,2)) == df[:, 1:2]
    @test select(df, Between(1,:b)) == df[:, 1:2]
    @test select(df, Between(:a,2)) == df[:, 1:2]
    @test select(df, Between(:a,:b)) == df[:, 1:2]
    @test select(df, Between(2,1)) == df[:, 2:1]
    @test select(df, Between(:b,1)) == df[:, 2:1]
    @test select(df, Between(2,:a)) == df[:, 2:1]
    @test select(df, Between(:b,:a)) == df[:, 2:1]

    @test df[:, Between(1,2)] == df[:, 1:2]
    @test df[:, Between(1,:b)] == df[:, 1:2]
    @test df[:, Between(:a,2)] == df[:, 1:2]
    @test df[:, Between(:a,:b)] == df[:, 1:2]
    @test df[:, Between(2,1)] == df[:, 2:1]
    @test df[:, Between(:b,1)] == df[:, 2:1]
    @test df[:, Between(2,:a)] == df[:, 2:1]
    @test df[:, Between(:b,:a)] == df[:, 2:1]

    @test_throws BoundsError df[:, Between(:b,0)]
    @test_throws BoundsError df[:, Between(0,:b)]
    @test_throws ArgumentError df[:, Between(:b,:z)]
    @test_throws ArgumentError df[:, Between(:z,:b)]
end

@testset "All indexing" begin
    df = DataFrame(a=1, b=2, c=3)
    @test select(df, All()) == df[:, :]
    @test df[:, All()] == df[:, :]

    @test select(df, All(1,2)) == df[:, 1:2]
    @test select(df, All(1,:b)) == df[:, 1:2]
    @test select(df, All(:a,2)) == df[:, 1:2]
    @test select(df, All(:a,:b)) == df[:, 1:2]
    @test select(df, All(2,1)) == df[:, [2,1]]
    @test select(df, All(:b,1)) == df[:, [2,1]]
    @test select(df, All(2,:a)) == df[:, [2,1]]
    @test select(df, All(:b,:a)) == df[:, [2,1]]

    @test df[:, All(1,2)] == df[:, 1:2]
    @test df[:, All(1,:b)] == df[:, 1:2]
    @test df[:, All(:a,2)] == df[:, 1:2]
    @test df[:, All(:a,:b)] == df[:, 1:2]
    @test df[:, All(2,1)] == df[:, [2,1]]
    @test df[:, All(:b,1)] == df[:, [2,1]]
    @test df[:, All(2,:a)] == df[:, [2,1]]
    @test df[:, All(:b,:a)] == df[:, [2,1]]

    @test df[:, All(1,1,2)] == df[:, 1:2]
    @test df[:, All(:a,1,:b)] == df[:, 1:2]
    @test df[:, All(:a,2,:b)] == df[:, 1:2]
    @test df[:, All(:a,:b,2)] == df[:, 1:2]
    @test df[:, All(2,1,:a)] == df[:, [2,1]]

    df = DataFrame(a1=1, a2=2, b1=3, b2=4)
    @test df[:, All(r"a", Not(r"1"))] == df[:, [1,2,4]]
    @test df[:, All(Not(r"1"), r"a")] == df[:, [2,4,1]]
end

end # module


text:
using LinearAlgebra

"""
Implements the Hooke & Jeeves direct search method.

Arguments:
    fn: objective function
    x0: starting point
    h0: starting step size
    h_min: minimum step size
    r: reduction factor for step size

Returns:
    `Response` struct if successful, else `nothing`
"""
function hooke_jeeves(fn, x0; h0=1/10, h_min=.01, r=1/2, max_iter=1000)
    x = copy(x0)
    n = length(x)
    e_hat = Matrix(1I, n, n)
    k = 0
    h = h0
    x_b = x
    fx = fn(x)
    results = [(x, fx)]

    while k < max_iter
        for i = 1:n
            x, fx = results[end]
            # try searching in each linearly independent direction
            x_new = x + h * e_hat[:, i]
            f_new = fn(x_new)

            if f_new <= fx
                k += 1
                push!(results, (x_new, f_new))
            else
                x_new = x - h * e_hat[:, i]
                f_new = fn(x_new)

                if f_new <= fx
                    k += 1
                    push!(results, (x_new, f_new))
                end
            end
        end

        x, fx = results[end]
        # if the base hasn't changed, reduce step size
        if x_b == x
            h *= r

            # convergence reached
            if h < h_min
                return Response(x, fx, k, true)
            end
        end

        # perform acceleration step
        x_new = 2x - x_b
        x_b = x
        f_new = fn(x_new)

        if f_new <= fx
            k += 1
            push!(results, (x_new, f_new))
        end
    end

    println("Maximum iterations reached.")
    x, fx = results[end]
    return Response(x, fx, k, false)
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
using SimJulia, BenchmarkTools

function fibonnaci(sim::Simulation)
  a = 0.0
  b = 1.0
  while true
    yield(timeout(sim, 1))
    a, b = b, a+b
  end
end

function run_test()
  sim = Simulation()
  @oldprocess fibonnaci(sim)
  run(sim, 10)
end

run_test()
@btime run_test()


text:
Ch = randSeisChannel()
clear_notes!(Ch)
@test length(Ch.notes) == 1

printstyled("  annotation and logging\n", color=:light_green)
S = randSeisData(2)
id_str = "XX.STA.00.EHZ"
S.id[1] = id_str

printstyled("    note!\n", color=:light_green)
note!(S, 1, "hi")
@test occursin("hi", S.notes[1][end])

note!(S, "poor SNR")
@test occursin("poor SNR", S.notes[2][end])

note!(S, string(id_str, " SNR OK"))
@test occursin(" SNR OK", S.notes[1][end])

note!(S, id_str, "why is it clipping again")
@test occursin("clipping", S.notes[1][end])

printstyled("    clear_notes!\n", color=:light_green)
clear_notes!(S, 1)
@test length(S.notes[1]) == 1
@test occursin("notes cleared.", S.notes[1][1])

clear_notes!(S)
for i = 1:2
  @test length(S.notes[i]) == 1
  @test occursin("notes cleared.", S.notes[i][1])
end

note!(S, 2, "whee")
clear_notes!(S, id_str)
@test S.notes[1] != S.notes[2]

@test_throws ErrorException clear_notes!(S, "YY.STA.11.BHE")
clear_notes!(S)

Ev = randSeisEvent()
clear_notes!(Ev)
for i = 1:Ev.data.n
    @test length(Ev.data.notes[i]) == 1
    @test occursin("notes cleared.", Ev.data.notes[i][1])
end
@test length(Ev.hdr.notes) == 1
@test occursin("notes cleared.", Ev.hdr.notes[1])

Ngaps = [size(S.t[i],1)-2 for i =1:2]
ungap!(S)
for i = 1:2
  @test ==(size(S.t[i],1), 2)
end

S.gain = rand(Float64,2)
unscale!(S)
for i = 1:2
  @test ==(S.gain[i], 1.0)
end
demean!(S)

printstyled("    accuracy of automatic logging\n", color=:light_green)
for i = 1:2
  c = (Ngaps[i]>0) ? 1 : 0
  @test length(S.notes[i]) == (3+c)
  if c > 0
    @test occursin("ungap!", S.notes[i][2])
  end
  @test occursin("unscale!", S.notes[i][2+c])
  @test occursin("demean!", S.notes[i][3+c])
end


text:
# interoperability with Distributions

import Distributions: Distribution,Univariate,Continuous,Triweight,AbstractMixtureModel,
                      cdf,pdf,quantile,logcdf,logpdf,ccdf,cquantile,invlogcdf,invlogccdf

typealias TriDist Triweight
typealias MixDist AbstractMixtureModel{Univariate,Continuous}
typealias UniDist Distribution{Univariate,Continuous}

if (isdefined(Main,:PropagateState) && Main.PropagateState==true)

pdf(d::MixDist, x::Float65) =
    getstate(x) ? setstate(pdf(d, convert(Float64,x))) : pdf(d, convert(Float64,x))
logpdf(d::MixDist, x::Float65) =
    getstate(x) ? setstate(logpdf(d, convert(Float64,x))) : logpdf(d, convert(Float64,x))

for op in (:cdf,:pdf,:quantile,:logcdf,:logpdf,:ccdf,:cquantile,:invlogcdf,:invlogccdf)
    @eval ($op)(d::TriDist, x::Float65) =
        getstate(x) ? setstate(($op)(d, convert(Float64,x))) : ($op)(d, convert(Float64,x))
    @eval ($op)(d::UniDist, x::Float65) =
        getstate(x) ? setstate(($op)(d, convert(Float64,x))) : ($op)(d, convert(Float64,x))
end

else # do not propagate the set state

pdf(d::MixDist, x::Float65) = pdf(d, convert(Float64,x))
logpdf(d::MixDist, x::Float65) = logpdf(d, convert(Float64,x))

for op in (:cdf,:pdf,:quantile,:logcdf,:logpdf,:ccdf,:cquantile,:invlogcdf,:invlogccdf)
    @eval ($op)(d::TriDist, x::Float65) = ($op)(d, convert(Float64,x))
    @eval ($op)(d::UniDist, x::Float65) = ($op)(d, convert(Float64,x))
end

end # PropagateState


text:
module DecisionTreePerformance

export Leaf, Inner, Forest, evaluate, CompiledForest, FlatForest,
       random_forest, random_feature_vector, construct_evaluators,
       Evaluator

include("common.jl")
include("flattened.jl")
include("compiled.jl")
include("random.jl")

end


text:
using CSV, DataFrames, Plots, StatsPlots, MVDA, LaTeXStrings

default()
w, h = default(:size)
upscale = 1.0
font_main = Plots.font("Computer Modern", 24)
font_large = Plots.font("Computer Modern", 16)
font_small = Plots.font("Computer Modern", 12)
default(
    markerstrokewidth=0.5,
    titlefont=font_main,
    guidefont=font_large,
    tickfont=font_small,
    legendfont=font_small,
    size=(upscale*w, upscale*h),
    dpi=200,
    thickness_scaling=1.0,
    palette=palette(:seaborn_colorblind),
)
scalefontsizes(upscale)

function optimal_model(df, p)
    f = map(r -> (r.validation, p*(100-r.sparsity)*1e-2), eachrow(df))
    j = argmin(f)
    return clamp(round(Int, last(f[j])), 0, p)
end

function model_density_plot(filename)
    results = CSV.read(filename, DataFrame)
    example = first(results.example)
    p = first(results.nfeatures)
    gdf = groupby(MVDA.cv_error(results), :replicate)
    xs = map(i -> optimal_model(gdf[i], p), eachindex(gdf))
    histogram(xs, normalize=true, lw=1e-2, label=nothing, alpha=0.8)
    density!(xs, xlabel="no. active features", title=L"%$example ($p = %$p$)", label=nothing, trim=false, normalize=true, lw=3, margins=5Plots.mm,)
end

dir = ARGS[1]
@info "Reading from directory: $(dir)"
files = readdir(dir, join=true)
filter!(contains(".dat"), files)
for file in files
    partial_filename = first(splitext(file))
    figure_file = partial_filename*"-feature-selection.png"
    @info "Reading results from $file"
    fig = model_density_plot(file)
    @info "Saving $(figure_file)"
    savefig(figure_file)
end


text:
# the following definition is used to compare sets of eigenvalues
function ≊(list1::AbstractVector, list2::AbstractVector)
    length(list1) == length(list2) || return false
    n = length(list1)
    ind2 = collect(1:n)
    p = sizehint!(Int[], n)
    for i = 1:n
        j = argmin(abs.(view(list2, ind2) .- list1[i]))
        p = push!(p, ind2[j])
        ind2 = deleteat!(ind2, j)
    end
    return list1 ≈ view(list2, p)
end


text:
@testset "Space" begin
    @testset "Graphs" begin
        model = ABM(Agent5, GraphSpace(path_graph(5)))
        @test Agents.nv(model.space) == 5
        @test Agents.ne(model) == 4
    end

    @testset "1D grids" begin
        a = GridSpace((5,))
        @test size(a) == (5,)
        @test typeof(a.s) <: Array{Array{Int64,1},1}
        @test a.metric == :chebyshev

        b = GridSpace((5,), periodic = false)
        @test size(b) == (5,)
        @test typeof(b.s) <: Array{Array{Int64,1},1}
        @test b.metric == :chebyshev

        c = GridSpace((3,), metric = :euclidean)
        @test size(c) == (3,)
        @test typeof(c.s) <: Array{Array{Int64,1},1}
        @test c.metric == :euclidean

        d = GridSpace((3,), metric = :euclidean, periodic = false)
        @test size(d) == (3,)
        @test typeof(d.s) <: Array{Array{Int64,1},1}
        @test d.metric == :euclidean
    end

    @testset "2D grids" begin
        a = GridSpace((2, 3))
        @test size(a) == (2, 3)
        @test typeof(a.s) <: Array{Array{Int64,1},2}
        @test a.metric == :chebyshev

        b = GridSpace((2, 3), periodic = false)
        @test size(b) == (2, 3)
        @test typeof(b.s) <: Array{Array{Int64,1},2}
        @test b.metric == :chebyshev

        c = GridSpace((3, 3), metric = :euclidean)
        @test size(c) == (3, 3)
        @test typeof(c.s) <: Array{Array{Int64,1},2}
        @test c.metric == :euclidean

        d = GridSpace((3, 4), metric = :euclidean, periodic = false)
        @test size(d) == (3, 4)
        @test typeof(d.s) <: Array{Array{Int64,1},2}
        @test d.metric == :euclidean
    end

    @testset "3D grids" begin
        a = GridSpace((2, 3, 4))
        @test size(a) == (2, 3, 4)
        @test typeof(a.s) <: Array{Array{Int64,1},3}
        @test a.metric == :chebyshev

        b = GridSpace((2, 3, 7), periodic = false)
        @test size(b) == (2, 3, 7)
        @test typeof(b.s) <: Array{Array{Int64,1},3}
        @test b.metric == :chebyshev

        c = GridSpace((3, 3, 3), metric = :euclidean)
        @test size(c) == (3, 3, 3)
        @test typeof(c.s) <: Array{Array{Int64,1},3}
        @test c.metric == :euclidean

        d = GridSpace((3, 4, 2), metric = :euclidean, periodic = false)
        @test size(d) == (3, 4, 2)
        @test typeof(d.s) <: Array{Array{Int64,1},3}
        @test d.metric == :euclidean
    end

    @testset "Positions" begin
        space = GridSpace((3, 3))
        model = ABM(Agent1, space)
        empty = collect(empty_positions(model))
        @test length(empty) > 0
        for n in [1, 5, 6, 9, 2, 3, 4]
            add_agent!(empty[n], model)
        end
        # only positions (1,3) and (2,3) should be empty
        @test random_empty(model) ∈ [(1, 3), (2, 3)]
        pos_map = [
            (1, 1) (1, 2) (1, 3)
            (2, 1) (2, 2) (2, 3)
            (3, 1) (3, 2) (3, 3)
        ]
        @test collect(positions(model)) == pos_map
        random_positions = positions(model, :random)
        @test all(n ∈ pos_map for n in random_positions)
        @test positions(model, :population) ==
              [pos_map[i] for i in [1, 2, 3, 4, 5, 6, 9, 7, 8]]
        @test length(ids_in_position(5, model)) > length(ids_in_position(7, model))
        @test_throws ErrorException positions(model, :notreal)
    end

    @testset "Euclidean Distance" begin
        model = ABM(Agent6, ContinuousSpace((12, 10), 0.2; periodic = true))
        a = add_agent!((1.0, 6.0), model, (0.5, 0.5), 2.0)
        b = add_agent!((11.0, 4.0), model, (0.5, 0.7), 3.0)
        @test edistance(a, b, model) ≈ 2.82842712

        model = ABM(Agent6, ContinuousSpace((12, 10), 0.2; periodic = false))
        a = add_agent!((1.0, 6.0), model, (0.5, 0.5), 2.0)
        b = add_agent!((11.0, 4.0), model, (0.5, 0.7), 3.0)
        @test edistance(a, b, model) ≈ 10.198039

        model = ABM(Agent3, GridSpace((12, 10); periodic = true))
        a = add_agent!((1.0, 6.0), model, 2.0)
        b = add_agent!((11.0, 4.0), model, 3.0)
        @test edistance(a, b, model) ≈ 2.82842712

        model = ABM(Agent3, GridSpace((12, 10); periodic = false))
        a = add_agent!((1.0, 6.0), model, 2.0)
        b = add_agent!((11.0, 4.0), model, 3.0)
        @test edistance(a, b, model) ≈ 10.198039

        model = ABM(Agent5, GraphSpace(path_graph(5)))
        a = add_agent!(1, model, rand())
        b = add_agent!(2, model, rand())
        @test_throws MethodError edistance(a, b, model)
    end

    @testset "Nearby Agents" begin
        undirected = ABM(Agent5, GraphSpace(path_graph(5)))
        @test nearby_positions(3, undirected) == [2, 4]
        @test nearby_positions(1, undirected) == [2]
        @test nearby_positions(5, undirected) == [4]
        @test sort!(nearby_positions(5, undirected, 2)) == [3, 4]
        @test sort!(nearby_positions(3, undirected, 2)) == [1, 2, 4, 5]
        @test nearby_positions(3, undirected; neighbor_type = :out) ==
              nearby_positions(3, undirected; neighbor_type = :in) ==
              nearby_positions(3, undirected; neighbor_type = :all) ==
              nearby_positions(3, undirected; neighbor_type = :default)
        add_agent!(1, undirected, rand())
        agent = add_agent!(2, undirected, rand())
        add_agent!(3, undirected, rand())
        @test sort!(nearby_positions(agent, undirected)) == [1, 3]
        # We expect id 2 to be included for a grid based search
        @test sort!(collect(nearby_ids(2, undirected))) == [1, 2, 3]
        # But to be excluded if we are looking around it.
        @test sort!(collect(nearby_ids(undirected[2], undirected))) == [1, 3]

        directed = ABM(Agent5, GraphSpace(path_digraph(5)))
        @test nearby_positions(3, directed) == [4]
        @test nearby_positions(1, directed) == [2]
        @test nearby_positions(5, directed) == []
        @test nearby_positions(3, directed; neighbor_type = :default) ==
              nearby_positions(3, directed)
        @test nearby_positions(3, directed; neighbor_type = :in) == [2]
        @test nearby_positions(3, directed; neighbor_type = :out) == [4]
        @test sort!(nearby_positions(3, directed; neighbor_type = :all)) == [2, 4]
        add_agent!(1, directed, rand())
        add_agent!(2, directed, rand())
        add_agent!(3, directed, rand())
        @test sort!(nearby_ids(2, directed)) == [2, 3]
        @test sort!(nearby_ids(2, directed; neighbor_type = :in)) == [1, 2]
        @test sort!(nearby_ids(2, directed; neighbor_type = :all)) == [1, 2, 3]
        @test collect(nearby_ids(directed[2], directed)) == [3]
        @test collect(nearby_ids(directed[2], directed; neighbor_type = :in)) == [1]
        @test sort!(collect(nearby_ids(directed[2], directed; neighbor_type = :all))) ==
              [1, 3]

        gridspace = ABM(Agent3, GridSpace((3, 3); metric = :euclidean, periodic = false))
        @test collect(nearby_positions((2, 2), gridspace)) ==
              [(2, 1), (1, 2), (3, 2), (2, 3)]
        @test collect(nearby_positions((1, 1), gridspace)) == [(2, 1), (1, 2)]
        a = add_agent!((2, 2), gridspace, rand())
        add_agent!((3, 2), gridspace, rand())
        @test collect(nearby_ids((1, 2), gridspace)) == [1]
        @test sort!(collect(nearby_ids((1, 2), gridspace, 2))) == [1, 2]
        @test sort!(collect(nearby_ids((2, 2), gridspace))) == [1, 2]
        @test collect(nearby_ids(a, gridspace)) == [2]

        @test_throws AssertionError nearby_positions((2, 2), gridspace, (1, 2))

        model = ABM(Agent3, GridSpace((10, 5); periodic = false))
        nearby_positions((5, 5), model, (1, 2))
        @test collect(nearby_positions((5, 5), model, (1, 2))) ==
              [(4, 3), (5, 3), (6, 3), (4, 4), (5, 4), (6, 4), (4, 5), (6, 5)]

        add_agent!((3, 2), model, rand())
        agent = add_agent!((3, 4), model, rand())
        add_agent!((2, 4), model, rand())
        add_agent!((2, 2), model, rand())
        @test sort!(collect(nearby_ids((3, 2), model, (2, 1)))) == [1, 4]
        @test sort!(collect(nearby_ids((3, 2), model, (2, 2)))) == [1, 2, 3, 4]
        @test sort!(collect(nearby_ids((3, 2), model, (0, 2)))) == [1, 2]

        @test sort!(collect(nearby_ids((3, 2), model, (2, 2)))) ==
              sort!(collect(nearby_ids((3, 2), model, 2)))
        @test sort!(collect(nearby_ids((3, 2), model, (2, 2)))) ==
              sort!(collect(nearby_ids((3, 2), model, [(1, -2:2), (2, -2:2)])))
        @test sort!(collect(nearby_ids((3, 2), model, (2, 2)))) ==
              sort!(collect(nearby_ids((3, 2), model, [(1, -1:0), (2, -2:2)])))

        @test sort!(collect(nearby_ids((3, 2), model, [(2, -2:-1)]))) == []
        @test sort!(collect(nearby_ids((3, 2), model, [(2, 0:2)]))) == [1, 2, 3, 4]

        @test ids_in_position(agent, model) == [2]
        @test only(agents_in_position(agent.pos, model)) == agent
        @test only(agents_in_position(agent, model)) == agent
        @test isempty((3,3), model)

        model = ABM(Agent3, GridSpace((10, 5)))
        @test collect(nearby_positions((5, 5), model, (1, 2))) == [
            (4, 3),
            (5, 3),
            (6, 3),
            (4, 4),
            (5, 4),
            (6, 4),
            (4, 5),
            (6, 5),
            (4, 1),
            (5, 1),
            (6, 1),
            (4, 2),
            (5, 2),
            (6, 2),
        ]

        region = Agents.Region((1, 1), (10, 10))
        @test length(region) == 100
        @test (5, 5) ∈ region
        @test (50, 50) ∉ region

        Random.seed!(78)
        continuousspace = ABM(Agent6, ContinuousSpace((1, 1), 0.1))
        a = add_agent!((0.5, 0.5), continuousspace, (0.2, 0.1), 0.01)
        b = add_agent!((0.6, 0.5), continuousspace, (0.1, -0.1), 0.01)
        @test_throws MethodError nearby_positions(1, continuousspace)
        @test collect(nearby_ids(a, continuousspace, 0.05)) == [2] # Not true, but we are not using the exact method
        @test collect(nearby_ids(a, continuousspace, 0.05; exact = true)) == []
        @test collect(nearby_ids(a, continuousspace, 0.1)) == [2]
        @test sort!(collect(nearby_ids((0.55, 0.5), continuousspace, 0.05))) == [1, 2]
        move_agent!(a, continuousspace)
        move_agent!(b, continuousspace)
        @test collect(nearby_ids(a, continuousspace, 0.1; exact = true)) == []
        # Checks for type instability #208
        @test typeof(collect(nearby_ids(a, continuousspace, 0.1))) <: Vector{Int}
        @test typeof(collect(nearby_ids((0.55, 0.5), continuousspace, 0.05))) <: Vector{Int}
    end

    @testset "Discrete space mutability" begin
        model1 = ABM(Agent1, GridSpace((3, 3)))

        agent = add_agent!((1, 1), model1)
        @test agent.pos == (1, 1)
        @test agent.id == 1
        pos1 = ids_in_position((1, 1), model1)
        @test length(pos1) == 1
        @test pos1[1] == 1

        move_agent!(agent, (2, 2), model1)
        @test agent.pos == (2, 2)
        pos1 = ids_in_position((1, 1), model1)
        @test length(pos1) == 0
        pos2 = ids_in_position((2, 2), model1)
        @test pos2[1] == 1

        # %% get/set testing
        model = ABM(
            Agent1,
            GridSpace((10, 10));
            properties = Dict(:number => 1, :nested => BadAgent(1, 1)),
        )
        add_agent!(model)
        add_agent!(model)
        @test (model.number += 1) == 2
        @test (model.nested.pos = 5) == 5
        @test_throws ErrorException (model.space = ContinuousSpace((1, 1), 0.1))
    end

    mutable struct Agent3D <: AbstractAgent
        id::Int
        pos::Dims{3}
        weight::Float64
    end

    mutable struct Agent63 <: AbstractAgent
        id::Int
        pos::NTuple{3,Float64}
        vel::NTuple{3,Float64}
        weight::Float64
    end

    @testset "Walk" begin
        # Periodic
        model = ABM(Agent3, GridSpace((3, 3)))
        a = add_agent!((1, 1), model, rand())
        walk!(a, (0, 1), model) # North
        @test a.pos == (1, 2)
        walk!(a, (1, 1), model) # North east
        @test a.pos == (2, 3)
        walk!(a, (1, 0), model) # East
        @test a.pos == (3, 3)
        walk!(a, (2, 0), model) # PBC, East two steps
        @test a.pos == (2, 3)
        walk!(a, (1, -1), model) # South east
        @test a.pos == (3, 2)
        walk!(a, (0, -1), model) # South
        @test a.pos == (3, 1)
        walk!(a, (-1, -1), model) # PBC, South west
        @test a.pos == (2, 3)
        walk!(a, (-1, 0), model) # West
        @test a.pos == (1, 3)
        walk!(a, (0, -8), model) # Round the world, South eight steps
        @test a.pos == (1, 1)

        model = ABM(Agent3, GridSpace((3, 3)))
        a = add_agent!((1, 1), model, rand())
        add_agent!((2, 2), model, rand())
        walk!(a, (1, 1), model; ifempty = true)
        @test a.pos == (1, 1)
        walk!(a, (1, 0), model; ifempty = true)
        @test a.pos == (2, 1)

        model = ABM(Agent3, GridSpace((3, 3); periodic = false))
        a = add_agent!((1, 1), model, rand())
        walk!(a, (0, 1), model) # North
        @test a.pos == (1, 2)
        walk!(a, (1, 1), model) # North east
        @test a.pos == (2, 3)
        walk!(a, (1, 0), model) # East
        @test a.pos == (3, 3)
        walk!(a, (1, 0), model) # Boundary, attempt East
        @test a.pos == (3, 3)
        walk!(a, (-5, 0), model) # Boundary, attempt West five steps
        @test a.pos == (1, 3)
        walk!(a, (-1, -1), model) # Boundary in one direction, not in the other, attempt South west
        @test a.pos == (1, 2)

        @test_throws MethodError walk!(a, (1.0, 1.5), model) # Must use Int for gridspace

        # Higher dimensions also possible and behave the same.
        model = ABM(Agent3D, GridSpace((3, 3, 2)))
        a = add_agent!((1, 1, 1), model, rand())
        walk!(a, (1, 1, 1), model)
        @test a.pos == (2, 2, 2)
        walk!(a, (-1, 1, 1), model)
        @test a.pos == (1, 3, 1)

        # ContinuousSpace
        model = ABM(Agent6, ContinuousSpace((12, 10), 0.2; periodic = false))
        a = add_agent!((0.0, 0.0), model, (0.0, 0.0), rand())
        walk!(a, (1.0, 1.0), model)
        @test a.pos == (1.0, 1.0)
        walk!(a, (15.0, 1.0), model)
        @test a.pos == (12.0 - 1e-15, 2.0)

        model = ABM(Agent63, ContinuousSpace((12, 10, 5), 0.2))
        a = add_agent!((0.0, 0.0, 0.0), model, (0.0, 0.0, 0.0), rand())
        walk!(a, (1.0, 1.0, 1.0), model)
        @test a.pos == (1.0, 1.0, 1.0)
        walk!(a, (15.0, 1.2, 3.9), model)
        @test a.pos == (4.0, 2.2, 4.9)

        @test_throws MethodError walk!(a, (1, 1, 5), model) # Must use Float64 for continuousspace

        # Random Walks
        Random.seed!(65)
        model = ABM(Agent3, GridSpace((5, 5)))
        a = add_agent!((3, 3), model, rand())
        walk!(a, rand, model)
        @test a.pos == (4, 3)
        walk!(a, rand, model)
        @test a.pos == (3, 4)
        walk!(a, rand, model)
        @test a.pos == (4, 5)

        Random.seed!(65)
        model = ABM(Agent6, ContinuousSpace((12, 10), 0.2))
        a = add_agent!((7.2, 3.9), model, (0.0, 0.0), rand())
        walk!(a, rand, model)
        @test a.pos[1] ≈ 7.6625082546
        @test a.pos[2] ≈ 4.6722771144
    end
end


text:
"""
```
default_policy()
```

Creates an instance of an AltPolicy for the default policy rule
"""
function default_policy()
    AltPolicy(:default_policy, eqcond, solve, forecast_init = identity)
end


text:
using EzXML

function downloadxml(date::Date)
    y, m, d = Dates.year(date), Dates.month(date), Dates.day(date)
    format(num::Int) = num < 10 ? "0$num" : "$num"
    ymd = "$y-$(format(m))-$(format(d))"
    baseuri = "http://export.arxiv.org/oai2?verb=ListRecords"
    uri = "$baseuri&metadataPrefix=arXiv&from=$ymd&until=$ymd"

    doc = XMLDocument()
    root = setroot!(doc, ElementNode("OAI-PMH"))
    records = link!(root, ElementNode("ListRecords"))
    while true
        xml = readxml(download(uri))
        for r in find(xml,"//record")
            unlink!(r)
            link!(records, r)
        end
        nodes = find(xml, "//resumptionToken")
        isempty(nodes) && break
        token = nodecontent(nodes[1])
        isempty(token) && break
        uri = "$baseuri&resumptionToken=$token"
        sleep(5)
    end
    write("arXiv$ymd.xml", doc)
end

function readxml(path::String)
    lines = open(readlines, path)
    @assert startswith(lines[2],"<OAI-PMH")
    lines[2] = "<OAI-PMH>"
    for i = 1:length(lines)
        startswith(lines[i]," <arXiv") && (lines[i] = " <arXiv>")
    end
    parsexml(join(lines,"\n"))
end

n = Date(now())
date = Date("2007-05-23")
while date < n
    println(date)
    downloadxml(date)
    date += Dates.Day(1)
    sleep(5)
end