text:
############################################################################################
# Constants
"RNG for sampling based solutions"
const _rng = Random.MersenneTwister(1)   # shorthand
Random.seed!(_rng, 1)

"Tolerance for stochastic solutions"
const _TOL = 1.0e-6

"Number of samples"
N = 10^3

############################################################################################
# Kernel and AD backends
kernels = [NUTS, HMC, MALA, Metropolis, Custom]
gradientkernels = [NUTS, HMC, MALA]

backends = [:ForwardDiff, :ReverseDiff, :ReverseDiffUntaped, :Zygote]
massmatrices = [MDense(), MDiagonal(), MUnit()]
generating = [UpdateFalse(), UpdateTrue()]

############################################################################################
# Initiate Base Model to check sampler

######################################## Model 1
data_uv = randn(_rng, 1000)

struct MyBaseModel <: ModelName end
myparameter = (μ = Param(0.0, Normal()), σ = Param(1.0, Gamma()))
mymodel = ModelWrapper(MyBaseModel(), myparameter)

#Create objective for both μ and σ and define a target function for it
myobjective = Objective(mymodel, data_uv, (:μ, :σ))
function (objective::Objective{<:ModelWrapper{MyBaseModel}})(θ::NamedTuple)
	@unpack data = objective
	lprior = Distributions.logpdf(Distributions.Normal(),θ.μ) + Distributions.logpdf(Distributions.Exponential(), θ.σ)
    llik = sum(Distributions.logpdf( Distributions.Normal(θ.μ, θ.σ), data[iter] ) for iter in eachindex(data))
	return lprior + llik
end
myobjective(myobjective.model.val)

function ModelWrappers.generate(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{MyBaseModel}})
    @unpack model, data = objective
    @unpack μ, σ = model.val
    return μ[1]
end

function ModelWrappers.predict(_rng::Random.AbstractRNG, objective::Objective{<:ModelWrapper{MyBaseModel}})
    @unpack model, data = objective
    @unpack μ, σ = model.val
	return rand(_rng, Normal(μ, σ))
end


text:
import Tables

function table(t::Any; args...) :: String
  rows = Tables.rows(t)

  grh.table(class="table"; args...) do
    grh.thead() do
      grh.tr() do
        grh.collection(Tables.columnnames(Tables.columns(t))) do item
          grh.th(scope="col") do
            item |> string
          end
        end
      end
    end *
    grh.tbody() do
      grh.collection(rows |> collect) do row
        grh.tr() do
          grh.collection([Tables.getcolumn(row, col) for col in Tables.columnnames(row)]) do item
            grh.td() do
              item |> string
            end
          end
        end
      end
    end
  end
end

ext: jl


text:
get_codigo_elemento(::Moeda)::String = "83"
get_valor_elemento(m::Moeda) = String(m.codigo)

const CONST_MOEDAS = [
	:999,
	:ADP,
	:AED,
	:AFN,
	:ALL,
	:AMD,
	:ANG,
	:AOA,
	:ARS,
	:ATS,
	:AUD,
	:AWG,
	:AZN,
	:BAM,
	:BBD,
	:BDT,
	:BEF,
	:BGN,
	:BHD,
	:BIF,
	:BMD,
	:BND,
	:BOB,
	:BOV,
	:BRL,
	:BSD,
	:BTN,
	:BWP,
	:BYR,
	:BZD,
	:CAD,
	:CDF,
	:CHE,
	:CHF,
	:CHW,
	:CLF,
	:CLP,
	:CNY,
	:COP,
	:COU,
	:CRC,
	:CUC,
	:CUP,
	:CVE,
	:CZK,
	:DEM,
	:DJF,
	:DKK,
	:DOP,
	:DZD,
	:EEK,
	:EGP,
	:ERN,
	:ESP,
	:ETB,
	:EUR,
	:FIM,
	:FJD,
	:FKP,
	:FRF,
	:GBP,
	:GEL,
	:GHS,
	:GIP,
	:GMD,
	:GNF,
	:GRD,
	:GTQ,
	:GYD,
	:HKD,
	:HNL,
	:HRK,
	:HTG,
	:HUF,
	:IDR,
	:IEP,
	:ILS,
	:INR,
	:IQD,
	:IRR,
	:ISK,
	:ITL,
	:JMD,
	:JOD,
	:JPY,
	:KES,
	:KGS,
	:KHR,
	:KMF,
	:KPW,
	:KRW,
	:KWD,
	:KYD,
	:KZT,
	:LAK,
	:LBP,
	:LKR,
	:LRD,
	:LSL,
	:LTL,
	:LUF,
	:LVL,
	:LYD,
	:MAD,
	:MDL,
	:MGA,
	:MKD,
	:MMK,
	:MNT,
	:MOP,
	:MRO,
	:MUR,
	:MVR,
	:MWK,
	:MXN,
	:MXV,
	:MYR,
	:MZN,
	:NAD,
	:NGN,
	:NIO,
	:NLG,
	:NOK,
	:NPR,
	:NZD,
	:OMR,
	:PAB,
	:PEN,
	:PGK,
	:PHP,
	:PKR,
	:PLN,
	:PTE,
	:PYG,
	:QAR,
	:RON,
	:RSD,
	:RUB,
	:RWF,
	:SAR,
	:SBD,
	:SCR,
	:SDG,
	:SEK,
	:SGD,
	:SHP,
	:SLL,
	:SOS,
	:SRD,
	:STD,
	:SVC,
	:SYP,
	:SZL,
	:THB,
	:TJS,
	:TMT,
	:TND,
	:TOP,
	:TRY,
	:TTD,
	:TWD,
	:TZS,
	:UAH,
	:UGX,
	:USD,
	:UYU,
	:UZS,
	:VEF,
	:VES,
	:VND,
	:VUV,
	:WST,
	:XAF,
	:XAU,
	:XBB,
	:XCD,
	:XDR,
	:XEU,
	:XFU,
	:XOF,
	:XPF,
	:YER,
	:ZAR,
	:ZMK,
	:ZMW,
	:ZWL
]


text:
using Unitful
using NLsolve

include("fricción.jl")
include("tuberias.jl")

mutable struct flujo
    L::typeof(1.0u"m")
    sh::String
    Dᵢ::typeof(1.0u"mm")
    Dₒ::typeof(1.0u"mm")
    ΔD::typeof(1.0u"mm")
    NPS::Float64
    Q::typeof(1.0u"l/hr")
    ṁ::typeof(1.0u"kg/hr")
    v::typeof(1.0u"m/s")
    ρ::typeof(1.0u"kg/m^3")
    𝜈::typeof(1.0u"m^2/s")
    μ::typeof(1.0u"cP")
    Re::Float64
    ε::typeof(1.0u"mm")
    f::Float64
    hf::typeof(1.0u"m")
end

function flujo(L,sh;Dᵢ,Dₒ,ΔD,NPS,Q,ṁ,v,ρ,𝜈,μ,Re,ε,f,hf)
    flujo(L,sh,Dᵢ,Dₒ,ΔD,NPS,Q,ṁ,v,ρ,𝜈,μ,Re,f,hf)
end

function flujo(L,sh;ṁ,ρ,μ,Dₒ,ε)
    NPS,Dᵢ,Dₒ,ΔD = tuberia(sh,Do=Dₒ)
    Q = ṁ/ρ
    v = Q/(π/4*Dᵢ^2)
    𝜈 = μ/ρ
    Re = upreferred(v*Dᵢ/𝜈)
    f = factor_friccion(Re,ε/Dᵢ)
    hf = upreferred(f*(L/Dᵢ)*(v^2 /(2*9.80665u"m/s^2")))
    flujo(L,sh,Dᵢ,Dₒ,ΔD,NPS,Q,ṁ,v,ρ,𝜈,μ,Re,ε,f,hf)
end

function flujo(L,sh;Q,ρ,𝜈,Dᵢ,ε)
    NPS,Dᵢ,Dₒ,ΔD = tuberia(sh,Di=Dᵢ)
    ṁ = Q*ρ
    v = Q/(π/4*Dᵢ^2)
    μ = 𝜈*ρ
    Re = upreferred(v*Dᵢ/𝜈)
    f = factor_friccion(Re,ε/Dᵢ)
    hf = upreferred(f*(L/Dᵢ)*(v^2 /(2*9.80665u"m/s^2")))
    flujo(L,sh,Dᵢ,Dₒ,ΔD,NPS,Q,ṁ,v,ρ,𝜈,μ,Re,ε,f,hf)
end

#=flujo(500u"m","40",0.2u"m^3/s",900u"kg/m^3",0.00001u"m^2/s",Dᵢ=200u"mm",ε=.26u"mm")

flujo(120u"m","80",40000u"kg/hr",12.493325597325967u"kg/m^3",2.655716090653184e-05u"Pa*s",Dₒ=150u"mm",ε=.0018u"inch")

.0018u"inch"
flujo(5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m")
flujo(5u"cm",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m",5u"m")

tuberia("40",Do=200u"mm")

factor_friccion(1, 1E-4)
factor_friccion(128000, .0013)
factor_friccion(1, 1E-4,función=Wood)
factor_friccion(128000, .0013,función=Moody)

laminar(100)=#

"""
    Reynolds(D,v,ρ,μ)
This is an example of Docstring. This function receives two
numbers x and y and returns the sum of the squares.
```math
x^2 + y^2
```
"""
Reynolds(D,v,ρ,μ) = upreferred(ρ*v*D/μ)
Reynolds(D,v,𝜈) = upreferred(v*D/𝜈)

#=Reynolds(2.5u"m", 0.25u"m/s", 1.1613u"kg/m^3", 1.9E-5u"Pa*s")
Reynolds(2.5u"m", 0.25u"m/s", 1.1613u"kg/m^3", 1.9E-5u"cP")
Reynolds(2.5u"m", 0.25u"m/s", 1.1613u"kg/m^3", 1.9E-5u"N*s/m^2")
Reynolds(2.5u"m", 0.25u"m/s", 1.636e-05u"m^2/s")=#

function hf(L::Unitful.Length,D::Unitful.Length,v::Unitful.Velocity,ρ::Unitful.Density,μ::Unitful.DynamicViscosity,eD;método::Function=Colebrook)
    upreferred(factor_friccion(Reynolds(D,v,ρ,μ), eD,función=método)*(L/D)*(v^2 /(2*9.80665u"m/s^2")))
end
function hf(L::Unitful.Length,D::Unitful.Length,v::Unitful.Velocity,𝜈::Unitful.KinematicViscosity,eD;método::Function=Colebrook)
    upreferred(factor_friccion(Reynolds(D,v,𝜈), eD,función=método)*(L/D)*(v^2 /(2*9.80665u"m/s^2")))
end

#=hf(400u"m",150u"mm",0.0157u"m/s",1.1613u"kg/m^3", 1.9E-3u"Pa*s",0)
hf(500u"m",200u"mm",6.4u"m/s",.00001u"m^2/s",.0013)=#


text:
struct VoidPlotter <: Plotter
end

plot( ::VoidPlotter, world::World, time ) = nothing
plot!( io, p, ::VoidPlotter, world::World, time ) = nothing


text:
mutable struct Soln
    plan::Vector{Tuple{Symbol,Vector{Any}}}
    robot_poses::Vector{Vector{Float64}}
    graph::Graph
    size::Int32
end


text:
### ============== ### ============== ###
### SIMPLE VICSEK MODEL 2D
### PERIODIC BOUNDARY CONDITION
### ============== ### ============== ###

function set_output_data_structure_vsk(path, N, ρ)

    parent_folder_path = "$(homedir())/art_DATA/$(path)"
    folder_path        = parent_folder_path * "/DATA/data_N_$(N)"
    reps_path          = folder_path * "/data_N_$(N)_rho_$(ρ)"

    try
        mkdir("$(homedir())/art_DATA")
    catch error
        println("Main data folder already exists")
    end

    try
        mkdir(parent_folder_path)
    catch error
        println("Parent folder already exists")
    end

    try
        mkdir(parent_folder_path * "/DATA")
    catch error
        println("Parent folder already exists")
    end

    try
        mkdir(folder_path)
    catch error
        println("Folder already exists")
    end

    try
        mkdir(reps_path)
    catch error
        println("Parameter folder already exists")
    end

    return reps_path
end

### ============== ### ============== ###
"""
    Box(L, M)
Simulation Box
# Constructor Arguments
* L -> size of box
* M -> number of cells per dimension

# Fields
* L -> linear size
* M -> number of cells per dimension

* bulk_cells   -> discrete coordinates
* bottom_cells -> discrete coordinates
* top_cells    -> discrete coordinates
* left_wall    -> discrete coordinates
* right_wall   -> discrete coordinates
* front_wall   -> discrete coordinates
* back_wall   -> discrete coordinates
* corners -> discrete coordinates

* center_cell -> real coordinates of cell center
* vel_cell -> average velocity of particles within cell

* p_per_cell -> Dict: cell_id => [particles id within cell]
"""
type Box
    L :: Float64
    M :: Float64

    bulk_cells   :: Array{Array{Int64,1},2}
    bottom_cells :: Array{Array{Int64,1},1}
    top_cells    :: Array{Array{Int64,1},1}
    left_cells    :: Array{Array{Int64,1},1}
    right_cells   :: Array{Array{Int64,1},1}
    corners   :: Array{Array{Int64,1},2}

    p_per_cell :: Dict{Array{Int,1}, Array{Int, 1}}

    vel_cell :: Array{Array{Float64,1},2}
    center_cell :: Array{Array{Float64,1},2}

    function Box(L, M)

        ### ============== ### ============== ###
        ## indices of cells

        bulk_cells = [[i,j] for i in 2:M-1, j in 2:M-1]

        bottom_cells = [[i,M] for i in 2:M-1]
        top_cells = [[i,1] for i in 2:M-1]

        left_cells = [[1,j] for j in 2:M-1]
        right_cells = [[M,j] for j in 2:M-1]

        corners = [[i,j] for i in [1,M], j in [1,M]]

        ### ============== ### ============== ###

        # particles id's in each cell
        p_per_cell = Dict{Array{Int,1}, Array{Int, 1}}()

        # average velocity per cell
        vel_cell = Array{Array{Float64,1}}(M,M)

        # Initialize velocity vector per cell and compute center of each cell
        for i in 1:length(vel_cell)
            vel_cell[i] = zeros(Float64, 2)
        end

        # position of center of each box
        center_cell = Array{Array{Float64,1}}(M,M)

        r = linspace(0., L, M+1)
        mp = collect(r[1:length(r) - 1] + 0.5 * step(r))

        for i in 1:M, j in 1:M
            # println(i,"\t", j, "\t", k)
            center_cell[i,j] = [mp[i], mp[j]]
        end

        new(L, float(M), bulk_cells, bottom_cells, top_cells, left_cells, right_cells, corners, p_per_cell, vel_cell, center_cell)
    end
end

### ============== ### ============== ###
"""
    Flock(N, L, v0, p, d)
Inertial Flock type
# Constructor Arguments
* N -> number of particles
* L -> size of box
* dt -> integration time step
* v0 -> particles speed
* r0 -> local interaction range
* η -> noise intensity

# Fields
* N -> number of particles
* v0 -> particles speed
* r0 -> local interaction range
* dt -> integration time step
* η -> noise intensity

* pos -> particles' positions
* vel -> particles' velocities
* v_r -> interaction vector
* k_r -> local connectivity

* p_cell_id -> particle's cell id
"""
type Flock
    N  :: Int64
    v0 :: Float64
    r0 :: Float64
    dt :: Float64
    η :: Float64

    pos :: Array{Array{Float64, 1}}
    vel :: Array{Array{Float64, 1}}
    v_r :: Array{Array{Float64, 1}}
    k_r :: Array{Float64, 1}

    p_cell_id :: Array{Array{Int64,1}}

    function Flock(N, L, dt, v0, r0, η)

        # particles' positions and velocities
        pos = [ [L*rand(), L*rand()] for i in 1:N ]
        vel = v0 * [ normalize([2*rand() - 1, 2*rand() - 1]) for i in 1:N ]

        # interaction vector
        v_r = [zeros(2) for i in 1:N]

        # connectivity per particle
        k_r = zeros(Float64, N)

        # box index for each particle in each dimension
        p_cell_id = [zeros(Int,2) for i in 1:N]

        new(N, v0, r0, dt, η, pos, vel, v_r, k_r, p_cell_id)
    end
end

## ============== ### ============== ###
function av_vel_cell(box, vel)
    for k in keys(box.p_per_cell)
        box.vel_cell[k[1], k[2]] = mean([vel[i] for i in box.p_per_cell[k]])
    end

end

### ============== ### ============== ###
function assign_cell(flock, box, cell_size)

    box.p_per_cell = Dict{Array{Int,1}, Array{Int, 1}}()

    for i in 1:flock.N
        flock.p_cell_id[i] = convert(Array{Int}, div.(floor.(flock.pos[i]), cell_size)) + 1

        haskey(box.p_per_cell, flock.p_cell_id[i]) ? push!(box.p_per_cell[flock.p_cell_id[i]], i) : box.p_per_cell[flock.p_cell_id[i]] = [i]
    end

end

### ============== ### ============== ###
function check_bulk_cells(flock, p_id, cell_id, p_per_cell)

    k_t = Vector{Float64}(9)
    v_t = Vector{Vector{Float64}}(9)

    c = 1

    # println("center cell: ", cell_id)

    for i in -1:1, j in -1:1
        key = [cell_id[1]+i, cell_id[2]+j]
        k_t[c], v_t[c] = part_bulk_cell_interaction(p_id, flock, get!(p_per_cell, key, []))
        c += 1
        # println(key, "\t", get!(p_per_cell, key, []))
    end

    sum(k_t) > 0. ? v_t = sum(v_t) ./ sum(k_t) : v_t = zeros(2)

    return v_t
end

### ============== ### ============== ###
function check_corners(flock, p_id, cell_id, p_per_cell, L, M)

    k_t = Vector{Float64}(9)
    v_t = Vector{Vector{Float64}}(9)

    c = 1
    # println("corner cell: ", cell_id)

    if cell_id[1] == 1 && cell_id[2] == 1
        for i in [M, 1, 2], j in [M, 1, 2]
            key = [i, j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    elseif cell_id[1] == 1 && cell_id[2] == M
        for i in [M, 1, 2], j in [M-1, M, 1]
            key = [i, j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    elseif cell_id[1] == M && cell_id[2] == 1
        for i in [M-1, M, 1], j in [M, 1, 2]
            key = [i, j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    elseif cell_id[1] == M && cell_id[2] == M
        for i in [M-1, M, 1], j in [M-1, M, 1]
            key = [i, j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    end

    sum(k_t) > 0. ? v_t = sum(v_t) ./ sum(k_t) : v_t = zeros(2)

    return v_t
end

### ============== ### ============== ###
function check_bottom_cells(flock, p_id, cell_id, p_per_cell, L, M)

    k_t = Vector{Float64}(9)
    v_t = Vector{Vector{Float64}}(9)

    c = 1
    # println("center cell: ", cell_id)

    # if in(cell_id[1], 2:M-1) && cell_id[2] == M
        for i in -1:1, j in [M-1, M, 1]
            key = [cell_id[1] + i, j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    # end

    sum(k_t) > 0. ? v_t = sum(v_t) ./ sum(k_t) : v_t = zeros(2)

    return v_t
end

### ============== ### ============== ###
function check_top_cells(flock, p_id, cell_id, p_per_cell, L, M)
    k_t = Vector{Float64}(9)
    v_t = Vector{Vector{Float64}}(9)

    c = 1
    # println("center cell: ", cell_id)

    # if in(cell_id[1], 2:M-1) && cell_id[2] == 1
        for i in -1:1, j in [M, 1, 2]
            key = [cell_id[1] + i, j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    # end

    sum(k_t) > 0. ? v_t = sum(v_t) ./ sum(k_t) : v_t = zeros(2)

    return v_t
end

### ============== ### ============== ###
function check_left_cells(flock, p_id, cell_id, p_per_cell, L, M)

    k_t = Vector{Float64}(9)
    v_t = Vector{Vector{Float64}}(9)

    c = 1
    # println("center cell: ", cell_id)

    # if cell_id[1] == 1 && in(cell_id[2], 2:M-1)
        for i in [M, 1, 2], j in -1:1
            key = [i, cell_id[2] + j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    # end

    sum(k_t) > 0. ? v_t = sum(v_t) ./ sum(k_t) : v_t = zeros(2)

    return v_t
end

### ============== ### ============== ###
function check_right_cells(flock, p_id, cell_id, p_per_cell, L, M)

    k_t = Vector{Float64}(9)
    v_t = Vector{Vector{Float64}}(9)

    c = 1
    # println("center cell: ", cell_id)

    # if cell_id[1] == M && in(cell_id[2], 2:M-1)
        for i in [M-1, M, 1], j in -1:1
            key = [i, cell_id[2] + j]
            k_t[c], v_t[c] = part_bound_cell_interaction(p_id, flock, get!(p_per_cell, key, []), L)
            c += 1
            # println(key, "\t", get!(p_per_cell, key, []))
        end
    # end

    sum(k_t) > 0. ? v_t = sum(v_t) ./ sum(k_t) : v_t = zeros(2)

    return v_t
end

### ============== ### ============== ###

function get_neighbors(p_id, flock, box)
    if flock.p_cell_id[p_id] in box.bulk_cells
        # println("bulk_cells")
        flock.v_r[p_id] = check_bulk_cells(flock, p_id, flock.p_cell_id[p_id], box.p_per_cell)
    elseif flock.p_cell_id[p_id] in box.bottom_cells
        # println("bottom_cells")
        flock.v_r[p_id] = check_bottom_cells(flock, p_id, flock.p_cell_id[p_id], box.p_per_cell, box.L, box.M)
    elseif flock.p_cell_id[p_id] in box.top_cells
        # println("top_cells")
        flock.v_r[p_id] = check_top_cells(flock, p_id, flock.p_cell_id[p_id], box.p_per_cell, box.L, box.M)
    elseif flock.p_cell_id[p_id] in box.left_cells
        # println("left_cells")
        flock.v_r[p_id] = check_left_cells(flock, p_id, flock.p_cell_id[p_id], box.p_per_cell, box.L, box.M)
    elseif flock.p_cell_id[p_id] in box.right_cells
        # println("right_cells")
        flock.v_r[p_id] = check_right_cells(flock, p_id, flock.p_cell_id[p_id], box.p_per_cell, box.L, box.M)
    elseif flock.p_cell_id[p_id] in box.corners
        # println("corners")
        flock.v_r[p_id] = check_corners(flock, p_id, flock.p_cell_id[p_id], box.p_per_cell, box.L, box.M)
    end
end

### ============== ### ============== ###
function part_bulk_cell_interaction(p_id, flock, cell_parts)

    k = 0.
    v_r = zeros(Float64, 2)

    for j in cell_parts
        if norm(flock.pos[p_id] - flock.pos[j]) > 0. && norm(flock.pos[p_id] - flock.pos[j]) <= flock.r0
            k += 1.
            @. v_r = flock.vel[j]
        end
    end

    return k, v_r
end

### ============== ### ============== ###
function part_bound_cell_interaction(p_id, flock, cell_parts, L)

    k = 0.
    v_r = zeros(Float64, 2)
    d_v = zeros(Float64, 2)

    for j in cell_parts

        # compute relative vector  taking into account PBC
        for i in eachindex(flock.pos[p_id])
            δ = flock.pos[p_id][i] - flock.pos[j][i]
            δ > 0.5 * L ? d_v[i] = δ - L : d_v[i] = δ
        end

        # check if particles are within interaction range
        if  norm(d_v) > 0. && norm(d_v) <= flock.r0
            k += 1.
            @. v_r = flock.vel[j]
        end
    end

    return k, v_r
end

### ============== ### ============== ###

function update_part(pos, vel, v_r, η, L)

    prop_angle = atan2(vel[2], vel[1])
    signal_angle = 0.0
    u_vel = zeros(Float64, 2)

    v_r[1] != 0.0 || v_r[2] != 0.0 ? signal_angle = atan2(v_r[2], v_r[1]) - prop_angle : signal_angle = 0.0

    total_angle = signal_angle + η * (2.0 * rand() * pi - pi);

    c = cos(total_angle)
    s = sin(total_angle)

    u_vel[1] = vel[1]*c - vel[2]*s;
    u_vel[2] = vel[1]*s + vel[2]*c;

    for i in eachindex(vel)

        vel[i] = u_vel[i]
        pos[i] += u_vel[i]

        # if(pos[i] > L * 0.5) pos[i] -= L end
        # if(pos[i] <= -L * 0.5) pos[i] += L end

        if(pos[i] > L ) pos[i] -= L end
        if(pos[i] <= 0.0) pos[i] += L end
    end

end

### ============== ### ============== ###
###          SYSTEM EVOLUTION         ###
### ============== ### ============== ###

function evolve_system(flock, box, cell_size)

    # compute each particle cell_id and assgin particles id to each box
    assign_cell(flock, box, cell_size)

    # for i in eachindex(flock.v_r)
    #     println(i, "\t",flock.vel[i])
    # end

    # Search interactions in adjacent cells
    for i in eachindex(flock.pos)
        get_neighbors(i, flock, box)
        # println("pass")
    end

    # for i in eachindex(flock.v_r)
    #     println(i, "\t",flock.v_r[i])
    # end

    # update particles' position
    broadcast(update_part, flock.pos, flock.vel, flock.v_r, flock.η, box.L)

end


text:
using OrderedCollections

using REPL.TerminalMenus

"""
    line_count_dir(directory::String)

Count lines of code in directory `directory`.
"""
function line_count_dir(directory, extensions, sort_results=true)
    extensions isa AbstractString && (extensions = (extensions,))
    all_files = [joinpath(root, f) for (root, dirs, files) in Base.Filesystem.walkdir(directory) for f in files]
    all_files_no_root = [f for (root, dirs, files) in Base.Filesystem.walkdir(directory) for f in files]
    filter!(x->any([endswith(x, y) for y in extensions]), all_files)
    n_lines = Dict([f=>0 for f in all_files])

    for f in all_files
        if isfile(f)
            n_lines[f] += countlines(f)
        else
            @warn "missing file $f"
        end
    end

    if sort_results
        n_lines = OrderedDict(sort(collect(n_lines), by=x->x[2]))
    end

    println("----------------- Number of lines per file")
    for (f,v) in n_lines
        println("$v \t $f")
    end
    println("----------------- Total number of lines of code:")
    println("n_lines = $(sum(values(n_lines)))")
    return n_lines
end

"""
    remove_dropbox_conflicts(code_dir::String)

Recursively delete files that match with "... conflicted copy ..."
"""
function remove_dropbox_conflicts(code_dir::AbstractString)
    all_files = [joinpath([root,f]...) for (root, dirs, files) in Base.Filesystem.walkdir(code_dir) for f in files];

    names_to_match = (
        "(LAPTOP-HN9BG775's conflicted copy",
        "(dhcp-103-207.caltech.edu's conflicted copy",
    )

    file_with_conflicts = filter(x->any(occursin(y, x) for y in names_to_match), all_files);

    println("Matched files found:")
    for file in file_with_conflicts
        println("$file")
    end

    options = ["yes", "no"]
    menu = RadioMenu(options, pagesize=4)
    choice = request("Would you like to delete these files?", menu)
    choice == -1 && error("Menu canceled")
    @assert choice == 1 || choice == 2
    if choice == 1
        for file in file_with_conflicts
          println("Deleting file $file")
          rm(file, force=true)
        end
    end
end

using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--directory"
            help = "Directory"
            arg_type = String
        "--extensions"
            help = "Extensions to include"
            arg_type = String
        "--sort_results"
            help = "Sort results by number of lines"
            arg_type = Bool
            default = true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    line_count_dir(parsed_args["directory"],
                   parsed_args["extensions"],
                   parsed_args["sort_results"])
end

if isinteractive()
    main()
end


text:
module aIFun

export greet

# Write your package code here.
greet() = println("Hi there.")

end


text:
## Computing marginal posterior probability density
## More on defining prior_vector, iterations, likelihood_vector to be included
likelihood = []
X = []
w = []
Z = []
function marginal_pd(N, prior_vector, iterations, likelihood_vector_current)
    for i in 1:1:iterations
        likelihood[i] = minimum(likelihood_vector_current) ## re-check this step
        X[i] = exp(-i/N)
        w[i] = X[i - 1] - X[i]
        Z = Z + likelihood[i]*w[i]
        ## save point with least likelihood as sample point with weight wi.
        ## update point with least likelihood with some MCMC steps according to prior,
        ## accepting only steps that keep the likelihood about Li
        end
    end
    return Z
end


text:
import Base: map
export Transformation
"""
    Transformation(T::Type, parser)
    Transformation{T}(transform, parser) where {T}
    Base.map(f::Function, Tc::Type, p::CombinedParser, a...)
    Base.map(f::Function, p::CombinedParser, a...)

Parser transforming result of a wrapped parser.
`a...` is passed as additional arguments to `f` (at front .

If `parser isa NamedParser`, transformation is done within the wrapped parser
(i.e. name applies to result-transforming parser).
"""
@auto_hash_equals struct Transformation{F,P,S,T} <: WrappedParser{P,S,T}
    transform::F
    parser::P
    Transformation{T}(transform, p_) where {T} =
        let p = parser(p_)
            new{typeof(transform),typeof(p),state_type(p),T}(transform, p)
        end
    function Transformation{T}(transform, p::NamedParser) where {T}
        tp = new{typeof(transform),typeof(p.parser),state_type(p.parser),T}(transform, p.parser)
        with_name(p.name, tp, p.doc)
    end
end
Transformation(T::Type, p) =
    Transformation{T}(T, p)

_deepmap_parser(f::Function,mem::AbstractDict,x::Transformation,a...;kw...) =
    Transformation{result_type(x)}(
        x.transform,
        deepmap_parser(f,mem,x.parser,a...;kw...))

children(x::Transformation) = children(x.parser)

function print_constructor(io::IO,x::Transformation)
    print_constructor(io,x.parser)
    print(io," |> map(")
    printstyled(io,x.transform, color=:bold)
    print(io,")")
end

"""
    Base.get(parser::Transformation{<:Function}, a...)
    Base.get(parser::Transformation{<:Type}, a...)

Function call `parser.transform(get(parser.parser,a...))`.
"""
function Base.get(parser::Transformation{<:Function}, sequence, till, after, i, state)
    v = get(parser.parser, sequence, till, after, i, state)
    parser.transform(v)
end
function Base.get(parser::Transformation{<:Type}, sequence, till, after, i, state)
    v = get(parser.parser, sequence, till, after, i, state)
    v isa parser.transform ? v : parser.transform(v)
end

export MatchedSubSequence
import Base: (!)
struct MatchedSubSequence end

function print_constructor(io::IO,x::Transformation{MatchedSubSequence})
    print_constructor(io,x.parser)
    printstyled(io," |> !", color=:bold)
end

@deprecate MatchedSubSequence(x) map(MatchedSubSequence(), parser(x))
export JoinSubstring
"""
    JoinSubstring = MatchedSubSequence

Deprecated but kept (because legacy `join(...; wrap=JoinSubstring)` syntax does not word with `@deprecate`).
"""
JoinSubstring = MatchedSubSequence

Base.map(::Type{MatchedSubSequence}, x::CombinedParser) =
    map(MatchedSubSequence(), x)

Base.map(::MatchedSubSequence, x::CombinedParser) = Base.map(MatchedSubSequence(), SubString{String}, x)
function Base.map(::MatchedSubSequence, result::Type, x::CombinedParser)
    Transformation{result}(MatchedSubSequence(),x)
end

function Base.get(x::Union{Transformation{MatchedSubSequence},
                           ConstantParser{<:AbstractString}},
                  sequence, till, after, i, state)
    li = _prevind(sequence,after)
    li<i ? "" : @inbounds SubString(sequence,i,li)
end

(!)(x::CombinedParser) = map(MatchedSubSequence,x)

using InternedStrings
intern(v::AbstractString) =
    InternedStrings.intern(v)::String

intern(v::Nothing) = ""

(!)(x::CombinedParser{<:Any,<:AbstractString}) =
    map(intern, x)

map(::typeof(intern), x::Transformation{typeof(intern)}) = x

"""
    MatchedSubSequence(x)
    (!)(x::CombinedParser)
    (!)(x::CombinedParser{<:Any,<:AbstractString})

Parser [`Base.map`](@ref) transformation getting
- either the matched `SubString`
- or an `InternedStrings.intern`ed copy thereof iif `result_type<:AbstractString` already.
Transformation does not evaluate `get(parser.transform,...)`.

```jldoctest
julia> Repeat(AnyChar())
.* AnyValue |> Repeat
::Vector{Char}

julia> !Repeat(AnyChar())
.* AnyValue |> Repeat |> !
::SubString{String}

julia> !!Repeat(AnyChar())
.* AnyValue |> Repeat |> ! |> map(intern)
::String

```

!!! warn
    MatchedSubSequence currently only support `SubString`, extension to Vector views requires minor effort.
    - `map(MatchedSubSequence(), SubArray{Int,1}, parser)` constructor is required (currently no sequence type is known at construction).
    - `get` function currently missing

"""
(!), MatchedSubSequence

"""
    map_constant(constant, p::CombinedParser)
    parser((p,constant)::Pair)

Construct a [`Base.map`](@ref) `Transformation{<:Constant}` resulting in `p` when calling [`get`](@ref) fast,
instead of computing result from state,
if `parser(p)` matches.

```jldocs
julia> parser("constant" | "fixed" => :constant)
|🗄 Either => :constant
├─ constant
└─ fixed
::Symbol
```

!!! note
    If the `Pair` key is a symbol, a [`NamedParser`](@ref) is created.
    ```jldocs
    julia> parser(:constant => "constant" | "fixed")
    |🗄 Either |> with_name(:constant)
    ├─ constant
    └─ fixed
    ::SubString{String}
    ```
"""
struct Constant{T}
    value::T
end
Base.show(io::IO, x::Constant) = show(io,x.value)

function print_constructor(io::IO,x::Transformation{<:Constant})
    print_constructor(io,x.parser)
    printstyled(IOContext(io, :compact => true)," => ",x.transform, color=:bold)
end

function map_constant(transform, p::CombinedParser)
    T=typeof(transform)
    Transformation{T}(Constant(transform), p)
end

parser(constant::Pair) =
    map_constant(constant.second, parser(constant.first))

function Base.get(parser::Transformation{<:Constant}, sequence, till, after, i, state)
    parser.transform.value
end

function _string(io::IO,x::Constant)
    print(io,"Constant(")
    show(io,x.value)
    print(io,")")
end
_string(io::IO,x::Function) = print(io,x)
_string(io::IO,x::Type) = print(io,x)

export MatchRange
"""
    MatchRange(p::CombinedParser)

Construct a [`Base.map`](@ref)`Transformation{UnitRange{Int}}` resulting in `p` when calling [`get`](@ref) fast,
Succeed iif `p` succeeds, if so results in sequence match index `UnitRange`.
Transformation does not evaluate `get(parser.transform,...)`.
"""
struct MatchRange
end
MatchRange(p::CombinedParser) =
    Transformation{UnitRange{Int}}(MatchRange(), p)

Base.show(io::IO, x::MatchRange) = print(io,"@")

function Base.get(parser::Transformation{MatchRange}, sequence, till, after, i, state)
    i:_prevind(sequence,after)
end

export IndexAt
"""
Struct for fast access to an index of a `Transformation`.

```jldoctest
julia> using CombinedParsers.Regexp

julia> p = re"(?:a|b*)."[1]
🗄 Sequence[1]
├─ |🗄 Either
│  ├─ a
│  └─ b*  |> Repeat
└─ [^\\n] ValueNotIn
::Union{Char, Vector{Char}}
```

See also [`getindex`](@ref), [`Sequence`](@ref).
"""
struct IndexAt{I}
    i::I
end

"""
    Base.get(parser::Transformation{<:IndexAt}, a...)

`getindex(get(parser.parser,a...).parser.transform)`
"""
function Base.get(parser::Transformation{IndexAt{I}}, sequence, till, after, i, state) where {I <: Integer}
    v = get(parser.parser,sequence, till, after, i, state)
    v[parser.transform.i]
end
function Base.get(parser::Transformation{IndexAt{Is}}, sequence, till, after, i, state) where {Is <: Union{Tuple, Vector, UnitRange}}
    tuple(get(parser.parser,sequence, till, after, i, state)[parser.transform.i]...)
end

function print_constructor(io::IO,x::Transformation{<:IndexAt})
    print_constructor(io,x.parser)
    printstyled(io,"[",x.transform.i,"]", color=:bold)
end

"""
    map(index::IndexAt, p::CombinedParser, a...)
    map(constant, p::CombinedParser, a...)

Parser matching `p`, transforming `p`s parsing results to `getindex(x,index)` or `constant`.

See also: [`get`](@ref), [`deepmap`](@ref)

"""
function Base.map(index::IndexAt{<:Integer}, p::CombinedParser)
    T=result_type(p)
    Transformation{fieldtype(T,index.i)}(index, p)
end
function Base.map(index::IndexAt{<:UnitRange}, p::CombinedParser)
    T=Tuple{fieldtypes(result_type(p))[index.i]...}
    Transformation{T}(index, p)
end

"""
    map(f::Function, p::CombinedParser, a...)

Parser matching `p`, transforming parsing results (`x`) with function `f(x,a...)`.

See also: [`get`](@ref), [`deepmap`](@ref)
"""
function Base.map(f::Function, p::CombinedParser, a...;
                  throw_empty_union=true)
    T = infer_result_type(f,Any,p,"call seq(function,type,parts...)",typeof.(a)...;
                          throw_empty_union=throw_empty_union)
    Transformation{T}(isempty(a) ? f : v -> f(v, a...), p)
end

function Base.map(f::Function, Tc::Type, p::CombinedParser, a...)
    T = infer_result_type(f,Tc,p,"call seq(function,type,parts...)",typeof.(a)...)
    Transformation{Tc}(isempty(a) ? f : v -> f(v, a...), p)
end

"""
    map(T::Type, p::CombinedParser, a...)

Parser matching `p`, transforming `p`s parsing result with constructor `T(x,a...)`.

See also: [`get`](@ref), [`deepmap`](@ref)
"""
function Base.map(Tc::Type, p::CombinedParser, a...)
    Transformation{Tc}(isempty(a) ? Tc : v -> Tc(a..., v), p)
end

function instance(Tc::Type, p::CombinedParser, a...)
    Transformation{Tc}((v) -> Tc(a..., v), p)
end

function instance(Tc::Type, p::CombinedParser)
    Transformation(Tc, p)
end

function Base.map(inner::CombinedParser, p::CombinedParser)
    Transformation{result_type(inner)}(s -> parse(inner,s), p)
end

Base.map(f::typeof(identity), p::CombinedParser) = p

@deprecate map(T::Type, f::Function, p::CombinedParser, a...) map(f,T,p,a...)
@deprecate instance(f::Function,p,a...) map(f,parser(p),a...)

"""
    infer_result_type(f::Function,Tc::Type,p::CombinedParser,onerror::AbstractString,ts::Type...; throw_empty_union=true)

Used by Parser Transformations to infer result type of a parser.
Throws error if type inference fails, if throw_empty_union=true.
"""
function infer_result_type(f::Function,Tc::Type,p::CombinedParser,onerror::AbstractString,ts::Type...; throw_empty_union=true)
    Ts = Base.return_types(f, tuple(result_type(p),ts...))
    isempty(Ts) && error("transformation type signature mismatch $f$(tuple(result_type(p),ts...))::$Ts<:$Tc")
    ( length(Ts) > 1 || Any <: first(Ts) ) && return Tc ##error(onerror*"  $f$(tuple(result_type(p),ts...))::$Ts<:$Tc")
    T = first(Ts)
    if throw_empty_union && T <: Union{}
        error("transformation type signature mismatch $f$(tuple(result_type(p),ts...))::$Ts<:$Tc")
    elseif T <: Tc
        T
    else
        @warn "type mismatch $f$(tuple(result_type(p),ts...))::$T<:$Tc"
        Tc
    end
end

export deepmap
"""
    deepmap(f, parser, predicate, a...; kw...)

Substitute all `sub_parser`s with [`Base.map`](@ref)`(f,sub_parser, a...; kw...)` iif
`dodeepmap(parser, predicate)`;
otherwise keep `sub_parser`

```jldocs
julia> p = re"(a+)b+"
🗄 Sequence |> regular expression combinator with 1 capturing groups
├─ (a+)  |> Repeat |> Capture 1
└─ b+  |> Repeat
::Tuple{Vector{Char}, Vector{Char}}

julia> p("aaabb")
(['a', 'a', 'a'], ['b', 'b'])

julia> deepmap(MatchedSubSequence, p, Capture)("aaabb")
("aaa", ['b', 'b'])

julia> deepmap(length, p, re"b+")("abb")
(['a'], 2)
```

Implementation is an example when the a custom leaf [`_deepmap`](@ref) method is useful and sufficient for [`deepmap_parser`](@ref).
"""
deepmap(f, parser, predicate, a...; kw...) =
    deepmap_parser(_deepmap, parser, predicate, f, a...; kw...)

"""
    _deepmap(parser, predicate, f, a...; kw...)

Implementation example when the a custom leaf [`_deepmap`](@ref) method is useful and sufficient for [`deepmap_parser`](@ref).
```julia
if dodeepmap(parser, predicate)
    map(f,parser, a...; kw...)
else
    parser
end
```
"""
function _deepmap(parser, predicate, f, a...; kw...)
    if dodeepmap(parser, predicate)
        map(f,parser, a...; kw...)
    else
        parser
    end
end

"""
    dodeepmap(parser, predicate)

`parser == predicate`.
Specialize for custom predicate type.

    dodeepmap(parser, predicate::Type)

`sub_parser isa predicate`

    dodeepmap(parser, predicate::Function)

`predicate(sub_parser)`

    dodeepmap(parser::NamedParser, predicate::Symbol)

`sub_parser.name == predicate`
"""
dodeepmap(parser, predicate) = parser == predicate
dodeepmap(parser, predicate::Function) = predicate(parser)
dodeepmap(parser, predicate::Type) = parser isa predicate
dodeepmap(parser::NamedParser, predicate::Symbol) = parser.name == predicate

export reinfer
"""
    reinfer(parser)

Run julia type inference again on a parser for optimization.
[`Either`](@ref)`{<:Vector}` parsers are converted to `Either{<:Tuple}`.

Implementation is an example when the a custom [`deepmap_parser`](@ref) method is useful.
"""
reinfer(parser) = deepmap_parser(_reinfer, parser)
_reinfer(parser) = parser
function deepmap_parser(::typeof(_reinfer),mem::AbstractDict,x::Either{<:Vector},a...;kw...)
    if haskey(mem,x)
        mem[x]
    else
        mem[x] = r = Either{result_type(x)}(Any[])
        for p in x.options
            push!(r,deepmap_parser(_reinfer,mem,p,a...;kw...))
        end
        mem[x] = Either(r.options...)
    end
end


text:
@testset "traces" begin
    @testset "Trajectory" begin
        t = Trajectory(; state = Vector{Int}(), reward = Vector{Bool}())
        @test (:state, :reward) == keys(t)
        @test haskey(t, :state)
        @test haskey(t, :reward)
        push!(t; state = 3, reward = true)
        push!(t; state = 4, reward = false)
        @test t[:state] == [3, 4]
        @test t[:reward] == [true, false]
        pop!(t)
        @test t[:state] == [3]
        @test t[:reward] == [true]
        empty!(t)
        @test t[:state] == Int[]
        @test t[:reward] == Bool[]
    end

    @testset "SharedTrajectory" begin
        t = SharedTrajectory(Int[], :state)
        @test (:state, :next_state, :full_state) == keys(t)
        @test haskey(t, :state)
        @test haskey(t, :next_state)
        @test haskey(t, :full_state)
        @test t[:state] == Int[]
        @test t[:next_state] == Int[]
        @test t[:full_state] == Int[]
        push!(t; state = 1, next_state = 2)
        @test t[:state] == [1]
        @test t[:next_state] == [2]
        @test t[:full_state] == [1, 2]
        empty!(t)
        @test t[:state] == Int[]
        @test t[:next_state] == Int[]
        @test t[:full_state] == Int[]
    end

    @testset "EpisodicTrajectory" begin
        t = EpisodicTrajectory(
            Trajectory(; state = Vector{Int}(), reward = Vector{Bool}()),
            :reward,
        )

        @test isfull(t) == false

        @test (:state, :reward) == keys(t)
        @test haskey(t, :state)
        @test haskey(t, :reward)
        push!(t; state = 3, reward = true)

        @test isfull(t) == true

        push!(t; state = 4, reward = false)
        @test t[:state] == [3, 4]
        @test t[:reward] == [true, false]
        pop!(t)
        @test t[:state] == [3]
        @test t[:reward] == [true]
        empty!(t)
        @test t[:state] == Int[]
        @test t[:reward] == Bool[]
    end

    @testset "CombinedTrajectory" begin
        t = CircularCompactPSALRTSALTrajectory(;
            capacity = 3,
            legal_actions_mask_size = (2,),
        )
        push!(t; state = 1, action = 1, legal_actions_mask = [false, false])
        push!(
            t;
            reward = 0.0f0,
            terminal = false,
            priority = 100,
            state = 2,
            action = 2,
            legal_actions_mask = [false, true],
        )

        @test t[:state] == [1]
        @test t[:action] == [1]
        @test t[:legal_actions_mask] == [false false]'
        @test t[:reward] == [0.0f0]
        @test t[:terminal] == [false]
        @test t[:priority] == [100]
        @test t[:next_state] == [2]
        @test t[:next_action] == [2]
        @test t[:next_legal_actions_mask] == [false true]'
        @test t[:full_state] == [1, 2]
        @test t[:full_action] == [1, 2]
        @test t[:full_legal_actions_mask] == [
            false false
            false true
        ]

        push!(
            t;
            reward = 1.0f0,
            terminal = true,
            priority = 200,
            state = 3,
            action = 3,
            legal_actions_mask = [true, true],
        )

        @test t[:state] == [1, 2]
        @test t[:action] == [1, 2]
        @test t[:legal_actions_mask] == [
            false false
            false true
        ]
        @test t[:reward] == [0.0f0, 1.0f0]
        @test t[:terminal] == [false, true]
        @test t[:priority] == [100, 200]
        @test t[:next_state] == [2, 3]
        @test t[:next_action] == [2, 3]
        @test t[:next_legal_actions_mask] == [
            false true
            true true
        ]
        @test t[:full_state] == [1, 2, 3]
        @test t[:full_action] == [1, 2, 3]
        @test t[:full_legal_actions_mask] == [
            false false true
            false true true
        ]

        pop!(t)

        @test t[:state] == [1]
        @test t[:action] == [1]
        @test t[:legal_actions_mask] == [false false]'
        @test t[:reward] == [0.0f0]
        @test t[:terminal] == [false]
        @test t[:priority] == [100]
        @test t[:next_state] == [2]
        @test t[:next_action] == [2]
        @test t[:next_legal_actions_mask] == [false true]'
        @test t[:full_state] == [1, 2]
        @test t[:full_action] == [1, 2]
        @test t[:full_legal_actions_mask] == [
            false false
            false true
        ]

        empty!(t)

        @test t[:state] == []
        @test t[:action] == []
        @test t[:reward] == []
        @test t[:terminal] == []
        @test t[:next_state] == []
        @test t[:next_action] == []
        @test t[:full_state] == []
        @test t[:full_action] == []
    end

    @testset "VectCompactSARTSATrajectory" begin
        t = VectCompactSARTSATrajectory(;
            state_type = Vector{Float32},
            action_type = Int,
            reward_type = Float32,
            terminal_type = Bool,
        )
        push!(t; state = Float32[1, 1], action = 1)
        push!(t; reward = 1.0f0, terminal = false, state = Float32[2, 2], action = 2)
        push!(t; reward = 2.0f0, terminal = true, state = Float32[3, 3], action = 3)

        @test t[:state] == [Float32[1, 1], Float32[2, 2]]
        @test t[:action] == [1, 2]
        @test t[:reward] == [1.0f0, 2.0f0]
        @test t[:terminal] == [false, true]
        @test t[:next_state] == [Float32[2, 2], Float32[3, 3]]
        @test t[:next_action] == [2, 3]
    end

    @testset "ElasticCompactSARTSATrajectory" begin
        t = ElasticCompactSARTSATrajectory(;
            state_type = Float32,
            state_size = (2,),
            action_type = Int,
            reward_type = Float32,
            terminal_type = Bool,
        )
        push!(t; state = Float32[1, 1], action = 1)
        push!(t; reward = 1.0f0, terminal = false, state = Float32[2, 2], action = 2)
        push!(t; reward = 2.0f0, terminal = true, state = Float32[3, 3], action = 3)

        @test t[:state] == Float32[1 2; 1 2]
        @test t[:action] == [1, 2]
        @test t[:reward] == [1.0f0, 2.0f0]
        @test t[:terminal] == [false, true]
        @test t[:next_state] == Float32[2 3; 2 3]
        @test t[:next_action] == [2, 3]

        @test pop!(t) ==
              (reward = 2.0f0, terminal = true, state = Float32[3.0, 3.0], action = 3)
        push!(t; reward = 1.0f0, terminal = false, state = Float32[2, 2], action = 2)
        @test t[:state] == Float32[1 2; 1 2]
        @test t[:action] == [1, 2]
        @test t[:reward] == [1.0f0, 1.0f0]
        @test t[:terminal] == [false, false]
        @test t[:next_state] == Float32[2 2; 2 2]
        @test t[:next_action] == [2, 2]

        empty!(t)
        @test length(t[:state]) == 0
    end
end


text:
# Autogenerated wrapper script for NLopt_jll for aarch64-apple-darwin
export libnlopt

JLLWrappers.@generate_wrapper_header("NLopt")
JLLWrappers.@declare_library_product(libnlopt, "@rpath/libnlopt.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnlopt,
        "lib/libnlopt.0.11.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Ledgers

using Base: UUID
using DemeNet: datadir

import Synchronizers
using Synchronizers: Record

import PeaceVote: load

const SyncLedger = Synchronizers.Ledger

abstract type AbstractLedger end

record!(ledger::AbstractLedger,fname::String,bytes::Vector{UInt8}) = error("Not impl")
records(ledger::AbstractLedger) = error("Not impl")

struct Ledger <: AbstractLedger
    dir::AbstractString
    ledger::SyncLedger
end

function Ledger(dir::AbstractString)
    ledger = SyncLedger(dir)
    return Ledger(dir,ledger)
end

Ledger(uuid::UUID) = Ledger(datadir(uuid) * "/ledger/")

record!(ledger::Ledger,fname::String,data::Vector{UInt8}) = push!(ledger.ledger, Record(fname,data))

### A dirty hack
function record!(ledger::Ledger,fname::String,data::Vector{UInt8},pid::Int)
    ledger.ledger.records[pid] = Record(fname,data)
end

function binary(x)
    io = IOBuffer()
    serialize(io,x)
    return take!(io)
end

records(ledger::Ledger) = ledger.ledger.records # One can pass that

import Base: dirname, basename
dirname(record::Record) = dirname(record.fname)
basename(record::Record) = basename(record.fname)

serve(port,ledger::Ledger) = Synchronizers.serve(port,ledger.ledger)

sync!(ledger::Ledger,syncport) = Synchronizers.sync(Synchronizers.Synchronizer(syncport,ledger.ledger))

import DemeNet: serialize, deserialize

using Pkg.TOML
using Base: UUID
using ..Types: Vote, Proposal, Braid
using DemeNet: Notary, DemeSpec, Deme, datadir, Signer, Certificate, Contract, Intent, Consensus, Envelope, ID, DemeID, AbstractID

deserialize(record::Record,type::Type) = deserialize(IOBuffer(record.data),type)

function load(ledger::AbstractLedger)
    chain = Union{Certificate,Contract}[]

    for record in records(ledger)
        if dirname(record)=="members"

            id = deserialize(IOBuffer(record.data),Certificate{ID})
            push!(chain,id)

        elseif dirname(record)=="braids"

            braid = deserialize(IOBuffer(record.data),Contract{Braid})
            push!(chain,braid)

        elseif dirname(record)=="votes"

            vote = deserialize(IOBuffer(record.data),Certificate{Vote})
            push!(chain,vote)

        elseif dirname(record)=="proposals"

            proposal = deserialize(IOBuffer(record.data),Certificate{Proposal})
            push!(chain,proposal)

        end
    end
    return chain
end

function getrecord(ledger::AbstractLedger,fname::AbstractString)
    for record in records(ledger)
        if record.fname == fname
            return record
        end
    end
end

serialize(ledger::AbstractLedger,data,fname::AbstractString) = record!(ledger,fname,binary(data))

serialize(ledger::AbstractLedger,id::Certificate{ID}) = record!(ledger,"members/$(id.document.id)",binary(id))

function serialize(ledger::AbstractLedger,vote::Certificate{Vote})
    pid = vote.document.pid
    uuid = hash(vote.signature)
    record!(ledger,"votes/$pid-$uuid",binary(vote)) ### We may also use length of the ledger
end

function serialize(ledger::AbstractLedger,proposal::Certificate{Proposal})
    msg = proposal.document.msg
    propid = hash(msg)
    uuid = hash(proposal.signature)
    record!(ledger,"proposals/$propid-$uuid",binary(proposal))
end

function serialize(ledger::AbstractLedger,braid::Contract{Braid})
    uuid = hash(braid.signatures)
    record!(ledger,"braids/$uuid",binary(braid))
end

function writebootrecord(deme::Deme,config)
    fname = "config/" * string(0,base=16,pad=8)
    bytes = binary(config)
    path = datadir(deme.spec.uuid) * "/ledger/" * fname
    mkpath(dirname(path))
    write(path, bytes)
end

function readbootrecord(ledger::Ledger,type::Type)
    rec = ledger.ledger.records[1]
    deserialize(rec,type)
end

export Ledger, record!, dirname, basename, serve

end


text:
@testset "cl6" begin
    data = open(readdlm, joinpath(@__DIR__, "data", "Cl6.txt"))

    for r in 1:size(data, 1)
        row      = data[r, :]
        x        = row[1]
        expected = row[2]

        @test ClausenFunctions.cl6(x) ≈ expected atol=1e-14
        @test ClausenFunctions.cl6(-x) ≈ -expected atol=1e-14
        @test ClausenFunctions.cl6(x - 2.0*pi) ≈ expected atol=1e-13
        @test ClausenFunctions.cl6(x + 2.0*pi) ≈ expected atol=1e-13
    end
end


text:
module TestNamespace

using Omega

function testtrackerror()
  x = normal(0.0, 1.0)
  x_ = cond(x, x <ₛ 1.0)
  lt = x <ₛ 1.0
  Omega.applytrackerr(x_, Omega.rand(defΩ()))
  gt = x >ₛ -1.0
  x__ = cond(x_, x >ₛ -1.0)
  Omega.applytrackerr(x__, Omega.rand(defΩ()))
end
testtrackerror()

end


text:
module Units

E  = 1e18
P  = 1e15
T  = 1e12
G  = 1e9
M  = 1e6
k  = 1e3
h  = 1e2
da = 1e1
d  = 1e-1
c  = 1e-2
m  = 1e-3
u  = 1e-6
n  = 1e-9
p  = 1e-12
f  = 1e-15
a  = 1e-18

# package code goes here
export E,
    P,
    T,
    G,
    M,
    k,
    h,
    da,
    d,
    c,
    m,
    u,
    n,
    p,
    f,
    a

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   International Geomagnetic Field Model v12.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
#   [2] https://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export igrf12syn

"""
    igrf12syn(isv::Int, date::Number, itype::Int, alt::Number, colat::Number, elong::Number; show_warns = true)

This is a Julia implementation of the official IGRF source code, which was
written in Fortran [2]. The input and output variables are exactly the same as
the ones described in the function `igrf12syn` in [2].

# Args

* `isv`: `0` if main-field values are required, `1` if secular variation values
         are required.
* `date`: Year A.D.
* `itype`: `1` if geodetic (spheroid), `2` if geocentric (sphere).
* `alt`: Height above sea level [km] if `itype = 1`, or distance from the center of
         Earth [km] if `itype = 2` (must be > 3485 km).
* `colat`: Colatitude (0 - 180) [˚].
* `elong`: East-Longitude (0 - 360) [˚].

# Keywords

* `show_warns`: Show warnings about the data (**Default** = `true`).

# Returns

* The north component [nT] if `isv = 0`, or [nT/year] if `isv = 1`.
* The east component [nT] if `isv = 0`, or [nT/year] if `isv = 1`.
* The vertical component [nT] if `isv = 0`, or [nT/year] if `isv = 1`.
* The total intensity if `isv = 0`, or rubbish if `isv = 1`.

# Remarks

* The `date` must be greater or equal to 1900 and less than or equal 2025.
  Notice that a warning message is printed for dates grated than 2020.

"""
function igrf12syn(isv::Int, date::Number, itype::Int, alt::Number,
                   colat::Number, elong::Number; show_warns = true)

    # Check the data, since this model is valid for years between 1900 and 2025.
    ( (date < 1900) || (date > 2025) ) &&
    error("This IGRF version will not work for years outside the interval [1900, 2025).")

    # Warn the user that for dates after the year 2020 the accuracy maybe reduced.
    show_warns && (date > 2020) &&
    @warn("The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2020.")

    # Declaration of variables
    gh          = gh_igrf12
    fn::Int     = 0
    gn::Int     = 0
    kmx::Int    = 0
    ll::Int     = 0
    nc::Int     = 0
    nmx::Int    = 0
    x::Float64  = 0.0
    y::Float64  = 0.0
    z::Float64  = 0.0
    t::Float64  = 0.0
    tc::Float64 = 0.0
    cl          = zeros(13)
    sl          = zeros(13)
    p           = zeros(105)
    q           = zeros(105)

    if date < 2015
        t  = 0.2*(date - 1900)
        ll = floor(Int,t)
        t  = t - ll

        # SH models before 1995.0 are only to degree 10.
        if date < 1995
            nmx = 10

            # nc = nmx*(nmx+2)
            nc  = 120

            ll  = nc*ll

            # kmx = (nmx+1)*(nmx+2)/2
            kmx = 66
        else
            nmx = 13

            # nc  = nmx*(nmx+2)
            nc  = 195

            ll  = floor(Int,0.2*(date - 1995))

            # 19 is the number of SH models that extend to degree 10.
            ll  = 120*19 + nc*ll

            # kmx = (nmx+1)*(nmx+2)/2
            kmx = 105
        end

        tc = 1 - t

        if isv == 1
            t  = +0.2
            tc = -0.2
        end
    else
        t  = date - 2015
        tc = 1.0

        if isv == 1
            t  = 1.0
            tc = 0.0
        end

        ll  = 3060
        nmx = 13

        # nc  = nmx*(nmx+2)
        nc = 195

        # kmx = (nmx+1)*(nmx+2)/2
        kmx = 105
    end

    r::Float64  = alt
    ct          = cos(colat*pi/180)
    st          = sin(colat*pi/180)
    cl[1]       = cos(elong*pi/180)
    sl[1]       = sin(elong*pi/180)
    cd::Float64 = 1.0
    sd::Float64 = 0.0
    l::Int      = 1
    m::Int      = 1
    n::Int      = 0

    if itype != 2
        # Conversion between geodetic and geocentric coordinates using WGS84
        # spheroid.
        a2    = 40680631.6
        b2    = 40408296.0
        one   = a2*st^2
        two   = b2*ct^2
        three = one + two
        rho   = sqrt(three)
        r     = sqrt(alt*(alt + 2*rho) + (a2*one + b2*two)/three)
        cd    = (alt + rho)/r
        sd    = (a2 - b2)/rho*ct*st/r
        one   = ct
        ct    = ct*cd -  st*sd
        st    = st*cd + one*sd
    end

    ratio = 6371.2/r
    rr    = ratio^2

    # Computation of Schmidt quasi-normal coefficients p and x(=q).
    # =============================================================

    p[1] = 1.0
    p[3] = st
    q[1] = 0.0
    q[3] = ct

    for k = 2:kmx
        # There is no need to check bounds here. The code guarantees that
        # everything is inside the right bounds. This increased the source code
        # performance by 13%.
        @inbounds begin
            if n < m
                m  = 0
                n  = n+1
                rr = rr*ratio
                fn = n
                gn = n-1
            end

            fm = m

            if (m == n)
                if k != 3
                    one   = sqrt(1 - 0.5/fm)
                    j     = k - n - 1
                    p[k]  = one*st*p[j]
                    q[k]  = one*(st*q[j] + ct*p[j])
                    cl[m] = cl[m-1]*cl[1] - sl[m-1]*sl[1]
                    sl[m] = sl[m-1]*cl[1] + cl[m-1]*sl[1]
                end
            else
                gmm   = m^2
                one   = sqrt(fn^2 - gmm)
                two   = sqrt(gn^2 - gmm)/one
                three = (fn + gn)/one
                i     = k - n
                j     = i - n + 1
                p[k]  = three*ct*p[i] - two*p[j]
                q[k]  = three*(ct*q[i] - st*p[i]) - two*q[j]
            end

            # Synthesis of x, y, and z in geocentric coordinates.
            lm  = ll + l
            one = (tc*gh[lm] + t*gh[lm+nc])*rr

            if m != 0
                two   = (tc*gh[lm+1] + t*gh[lm+nc+1])*rr
                three = one*cl[m] + two*sl[m]
                x     = x + three*q[k]
                z     = z - (fn + 1)*three*p[k]

                if st != 0
                    y = y + (one*sl[m] - two*cl[m])*fm*p[k]/st
                else
                    y = y + (one*sl[m] - two*cl[m])*q[k]*ct
                end

                l = l + 2
            else
                x = x + one*q[k]
                z = z - (fn + 1)*one*p[k]
                l = l + 1
            end

            m = m + 1
        end
    end

    # Conversion to coordinate system specified by itype.
    # ===================================================
    one   = x
    x     = x*cd +   z*sd
    z     = z*cd - one*sd
    f     = sqrt(x^2 + y^2 + z^2)

    return x,y,z,f
end


text:
using Revise, Pkg, DataFrames, Finometrics, BenchmarkTools

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

Pkg.activate(pwd())
using Capacity

#need this since we are uing generated functions
if @isdefined(FIRST_RUN)
  revise(Capacity)
end

const FIRST_RUN = true
#=throw("1) Determine which spec is better for the tables (see tables0 vs tables)
2) verify the SEs for the AR(1) approach- NOTE- verify the time series are no longer stationary
- maybe put in a check to prevent this from occuring in thef uture? like an INF if the coef is ≥1?
3) Convert the multiple solution simulation graphs into slides
4) Outline the anlytical solutions")=#
#throw("Need to run the funds with refresh
#also, need to do full refresh runs to account for the error in 6/9 month
#momentum")
#throw("Reset the size acceptance threshold to 20% (or smaller)")
#=throw( "- why are MF assets higher than total assets? Doesn't seem right- maybe use the corridor method,
 or try changing up the factors?
   -Try the new domestic only factor
   -Determine if we really need time FE, or an intercept, or nothing. Update the deck accordingly")=#

#throw("now run the new and imrpoved AR(1)")

@warn("
check LARS todo
check timing of shocks to momentum assets

  Why are g and ng producing the same results? repkicate and test
")

@warn("
 - check the mutual fund file and make sure the betas and assets are ok- consider trimming the dates
 via the global setting if needed, may need to aquire new mf file if necessary. Check and
   compare the winsorized results. If it doesn't make a difference, may want to disable
   or shortcircuit the winsorization
 -compute standard errors on G (weighted abs std deviation betas * A / (total assets))
  - then maybe use the delta method to compute the SE of the ratio
 - try the mutual fund assets approach, either just the growth or the bootstrap approach
 - try OOS regressions
 - consider alternate measures of participation
 - check correlation of D/P ratio w/ the lmt measure
 - long run resutls are suspect due to the persistence of the log measure
 - longer horizons- try running the monthyl data on 35/60 years and see what happens
 - short run results seem pretty good- note using G instead of aG makes sense as
 prices will rise with demand
 - the rise in prices is tough to capture because contemperaneous outflows may
 coincide with crisis periods that often have high momentum returns

")
@info("
# Important note: CUDA v2.6.3 runs fastest, v3.x.x runs slow so far.
The easiest way to transition between version is to, in the main environment,
switch to the later version.
The results seem basically equivelent if allowed to run to completion,
  although in one instance I didn't get convergence with the new a version
wrt to the reason for the gap, I tried 1) checking the broadcast
  2) the regression. Implicity, the projection isn't responsible for the
  performance regression either since the gap still exists in the
  level version. Its not limit either.
  It is stable, at leas twith the old version.
  Note the latest run is microscopicalyl different from a bit back,
  but the difference is so tiny w/ correlation 1.0, beta ≈ 1.0. More improtantly,
  the answer is stable as written with a full pass.")
@info "WISHLIST (in rough priority)
REPLICATION and VERIFICATION!!! Mutual funds or momentum?
Mutual funds- some decent correlations. Could try some randomization for placebo test
Ultimately however, if the momentum emasure works well, that is strong evidence that
outweighs even placebo tests
"
#=NOTE OLD @info "I recced this against Lewellen 2015. Returns in Lewellen are a bit lower, and the b/m metric
sd is a bit high. Size and investment (using inv4) are ok. N seems out of control, but who knows
how they are calculating that."=#

#=@ NOTE OLD @info "MCMC- Seems like A₀ is mostly fixed- probably opportunities for performance improvements
  via profiling as well as additional verificaiton that the differences between the rep
  and actual fall with sample size"=#
#@ NOTE OLD @info "a GMM style EIV model will be much more efficient (check rossi slides)"
#note- the uncontrolled version greatly increases power and correlation

function capacityscript()
  #Capacity.gz2zstd("data\\crsp\\crspd35y")
  #Capacity.gz2zstd("data\\crsp\\crspd30y")
  #Capacity.gz2zstd("data\\crsp\\crspd23y")
  #Capacity.gz2zstd("data\\crsp\\crspd20y")
  #Capacity.gz2zstd("data\\crsp\\crspd10y")
  #Capacity.gz2zstd("data\\preliminary\\prelimcwarspm")5
  #Capacity.gz2zstd("data\\crsp\\crspd2y")
  #Capacity.gz2zstd("data\\comp\\compq")
  #Capacity.gz2zstd("data\\ibes\\ibes")
  #Capacity.gz2zstd("data\\hfr\\hfralive")
  #Capacity.gz2zstd("data\\hfr\\hfrdead")

  Capacity.getparameters!() #WARNING- mutates parameter dictionary to match excel parameter file
  #Capacity.formpreliminaryfunds()z
  #Capacity.formpreliminarybeta()

  #Capacity.computemeasure()
  Capacity.analyze()

  #Capacity.capacityio()

  #Capacity.comomentumtable()
  #Capacity.comomentumgraphs()
  #Capacity.fundgraphs()

end

@time capacityscript()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ c365f2f2-f361-11ea-39ae-799739acc8eb
begin
	using Pkg
	Pkg.activate("synbench_env", shared=true)
	#Pkg.add("Fairness")
	#Pkg.add("MLJBase")
	#Pkg.add("MLJ")
	#Pkg.add("StatsPlots")
	#Pkg.add("DecisionTree")
	#Pkg.add("DataFrames")
	using Fairness, MLJ, MLJBase
	using DataFrames, StatsPlots
end

# ╔═╡ f3b75fa4-f361-11ea-36f5-fdef03dddf63
# Score on training data to assess model quality
function get_scores(m, X, y, rows, set)
  yhat = predict(m, rows=rows)
  if typeof(yhat[1]) <: MLJBase.UnivariateFinite
    yhat = MLJBase.mode.(yhat)
  end
  ft = fair_tensor(yhat, y[rows], X[rows, :z])
  df = disparity(
    [accuracy, false_positive_rate],
    ft, refGrp="A",
    func=(x, y)->(x-y)
  )
  df[:,"accuracy"] .= accuracy(ft)
  df[:,"false_positive_rate"] .= false_positive_rate(ft)
  df[:,"set"] .= set
  return df
end

# ╔═╡ f4db5732-f361-11ea-3359-b72ac75cebb4
function run_synthetic_experiment(md, ntrain = 10000, genFun = genZafarData; genargs...)
  (method, mod) = md

  ntest = 1000
  repls = 30

  # Train data
  n = ntrain + (ntest * repls)
  X,y = genFun(n; genargs...)

  # Fit classifier
  train, test = partition(eachindex(y), ntrain/n, shuffle=true)
  m = machine(mod, X, y)
  fit!(m, rows=train)

  # Compute scores on train data
  df = get_scores(m, X, y, train, "insample")

  # Score on different test subsets
  tsts = [test[i:(i+ntest-1)] for i in 1:ntest:length(test)]
  for tst in tsts
    append!(df, get_scores(m, X, y, tst, "oob"))
  end
  df = df[df.labels .== "B",:]
  select!(df, Not(1))
  df[:, "n"] .= ntrain
  df[:, "method"] .= method
  return df
end

# ╔═╡ f5a67cc8-f361-11ea-0de5-33599ea9d491
# Create a list of classifiers to compare
begin
	@load RandomForestClassifier pkg=DecisionTree
	model = @pipeline ContinuousEncoder RandomForestClassifier
	m1 = ReweighingSamplingWrapper(classifier=model, grp=:z)
	m2 = LinProgWrapper(classifier=m1, grp=:z, measure=false_positive_rate)
	m3 = EqOddsWrapper(classifier=m1, grp=:z)
	clf = Dict("RF" => model, "RF_RW" => m1, "RF_RW_LP" => m2, "RF_RW_EO" => m3)
	keys(clf)
end

# ╔═╡ f65f0eb4-f361-11ea-19d8-d976d9779175
begin
	# Test on Zafar et al., 2017 Dataset
	res = DataFrame()
	for c in clf
    	append!(res, run_synthetic_experiment(c, 20000))
	end
end

# ╔═╡ 1ee229b6-f362-11ea-2ab7-671cb85e12f6
begin
	@df res[res.set .!= "insample",:] scatter(:accuracy,  :false_positive_rate_disparity, group = :method)
	xlabel!("Accuracy")
	ylabel!("Δ False Positive Rate")
end

# ╔═╡ 1fbbd274-f362-11ea-079a-35c51d9bd87c
begin
	res2 = DataFrame()
	for c in clf
    	append!(res2, run_synthetic_experiment(c, 10000, genBiasedSampleData; sampling_bias = 0.5))
	end
end

# ╔═╡ 5e2a5706-f362-11ea-0781-1d3f16ee8376
begin
	@df res2[res2.set .!= "insample",:] scatter(:accuracy,  :false_positive_rate_disparity, group = :method)
	xlabel!("Accuracy")
	ylabel!("Δ False Positive Rate")
end

# ╔═╡ Cell order:
# ╠═c365f2f2-f361-11ea-39ae-799739acc8eb
# ╠═f3b75fa4-f361-11ea-36f5-fdef03dddf63
# ╠═f4db5732-f361-11ea-3359-b72ac75cebb4
# ╠═f5a67cc8-f361-11ea-0de5-33599ea9d491
# ╠═f65f0eb4-f361-11ea-19d8-d976d9779175
# ╠═1ee229b6-f362-11ea-2ab7-671cb85e12f6
# ╠═1fbbd274-f362-11ea-079a-35c51d9bd87c
# ╠═5e2a5706-f362-11ea-0781-1d3f16ee8376


text:
function adjacency_matrix(adj::AbstractMatrix, T::DataType=eltype(adj))
    m, n = size(adj)
    (m == n) || throw(DimensionMismatch("adjacency matrix is not a square matrix: ($m, $n)"))
    T.(adj)
end

adjacency_matrix(fg::FeaturedGraph, T::DataType=eltype(graph(fg))) = adjacency_matrix(graph(fg), T)

Zygote.@nograd adjacency_matrix

function degrees(fg::FeaturedGraph, T::DataType=eltype(graph(fg)); dir::Symbol=:out)
    degrees(graph(fg), T; dir=dir)
end

function degree_matrix(fg::FeaturedGraph, T::DataType=eltype(graph(fg)); dir::Symbol=:out)
    degree_matrix(graph(fg), T; dir=dir)
end

function inv_sqrt_degree_matrix(fg::FeaturedGraph, T::DataType=eltype(graph(fg)); dir::Symbol=:out)
    inv_sqrt_degree_matrix(graph(fg), T; dir=dir)
end

function laplacian_matrix(fg::FeaturedGraph, T::DataType=eltype(graph(fg)); dir::Symbol=:out)
    laplacian_matrix(graph(fg), T; dir=dir)
end

function normalized_laplacian(fg::FeaturedGraph, T::DataType=eltype(graph(fg)); selfloop::Bool=false)
    normalized_laplacian(graph(fg), T; selfloop=selfloop)
end

function scaled_laplacian(fg::FeaturedGraph, T::DataType=eltype(graph(fg)))
    scaled_laplacian(graph(fg), T)
end

## Inplace operations

function laplacian_matrix!(fg::FeaturedGraph, T::DataType=eltype(graph(fg)); dir::Symbol=:out)
    if fg.matrix_type == :adjm
        fg.graph .= laplacian_matrix(graph(fg), T; dir=dir)
        fg.matrix_type = :laplacian
    end
    fg
end

function normalized_laplacian!(fg::FeaturedGraph, T::DataType=eltype(graph(fg)); selfloop::Bool=false)
    if fg.matrix_type == :adjm
        fg.graph .= normalized_laplacian(graph(fg), T; selfloop=selfloop)
        fg.matrix_type = :normalized
    end
    fg
end

function scaled_laplacian!(fg::FeaturedGraph, T::DataType=eltype(graph(fg)))
    if fg.matrix_type == :adjm
        fg.graph .= scaled_laplacian(graph(fg), T)
        fg.matrix_type = :scaled
    end
    fg
end


text:
using splines

println("1D interpolation")
orders = [100]
smin = [0.0]
smax = [1.0]
N_tries = 10

f(x) = sin(x*10)
# 1d test
a = linspace(0,1,orders[1])
b = f(a)
b
di = 1
M = length(a)

coefs = filter_coeffs(smin,smax,orders,b)

#s = linspace(-0.1,1.1,100000)
s = linspace(0.01,0.99,100000)

s = copy(s'')

interp_values = eval_UC_spline(smin,smax,orders,coefs,s)
val, grad = eval_UC_spline_G(smin,smax,orders,coefs,s)

interp_values = eval_UC_spline(smin,smax,orders,coefs,s)
val, grad = eval_UC_spline_G(smin,smax,orders,coefs,s)

fun = interpolant_cspline(smin,smax,orders,b)
vals = fun(s)

println("Error: ", maximum(maximum(abs(interp_values - vals )  )))

tic()
for i = 1:N_tries
    interp_values = eval_UC_spline(smin,smax,orders,coefs,s)
end

true_values = f(s)

println("Maximum interpolation error : ", norm(true_values-interp_values,Inf))
toc()

# Test multi splines:

CC = zeros(2,size(coefs,1))
CC[1,:] = coefs   # first spline
CC[2,:] = coefs/2   # second spline
interp_multi_values = eval_UC_multi_spline(smin,smax,orders,CC,s)

#using PyPlot
#plot(s,sol)
#plot(s,true_values)
#plot(a,b,"o")
#side effect of pyplot: 0.0 not recognized as a constant

# 2d test
println("2D interpolation")
smin = [0.0,0.0]
smax = [1.0,1.0]
orders = [10,20]
a1 = linspace(0,1,orders[1])
a2 = linspace(0,1,orders[2])
b = [sin(i+j) for i=a1, j=a2]
b = convert(Array{Float64},b)
#b = reshape(b, orders[1], orders[2])
#di = [1.0,1.0]
coefs = filter_coeffs(smin,smax,orders,b)

forders = [1000,1000]
N = prod(forders)
mgrid = [ [x y] for x=linspace(0,1,forders[1]), y=linspace(0,1,forders[2])]

vgrid = vcat(mgrid...)
interp_vals = eval_UC_spline(smin,smax,orders,coefs,vgrid)

tic()
for i = 1:10
interp_vals = eval_UC_spline(smin,smax,orders,coefs,vgrid)
end
toc()
true_vals = [sin(x+y) for x=linspace(0,1,forders[1]), y=linspace(0,1,forders[2])]
true_vals = convert(Array{Float64},true_vals[:])

println("Maximum interpolation error : ", norm(true_vals-interp_vals,Inf))

# 3d test
println("3D interpolation")

a = [0.0, 0.0, 0.0]
b = [1.0, 1.0, 1.0]
orders = [20,20,20]

a1 = linspace(0,1,orders[1])
a2 = linspace(0,1,orders[2])
a3 = linspace(0,1,orders[3])

mat = [(i+j+k)^2 for i=a1, j=a2, k=a3]
mat = convert(Array{Float64},mat)

coefs = filter_coeffs(a,b,orders,mat)

forders = [100,100,100]
N = prod(forders)

mgrid = [ [x y z] for x=linspace(0,1,forders[1]), y=linspace(0,1,forders[2]), z=linspace(0,1,forders[3])]

vgrid = vcat(mgrid...)
interp_vals = eval_UC_spline(a,b,orders,coefs,vgrid)

tic()
for i = 1:10
interp_vals = eval_UC_spline(a,b,orders,coefs,vgrid)
end
toc()
true_vals = [(x+y+z)^2 for x=linspace(0,1,forders[1]), y=linspace(0,1,forders[2]),  z=linspace(0,1,forders[3])]
true_vals = convert(Array{Float64},true_vals[:])

println("Maximum interpolation error : ", norm(true_vals-interp_vals,Inf))


text:
using GeoInterfaceRFC
using Test

include("test_primitives.jl")


text:
location_str(::Type{Face})    = "Face"
location_str(::Type{Center})  = "Center"
location_str(::Type{Nothing}) = "⋅"

show_location(LX, LY, LZ) = "($(location_str(LX)), $(location_str(LY)), $(location_str(LZ)))"

show_location(field::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = show_location(LX, LY, LZ)

short_show(m::Missing) = "$m"

short_show(field::AbstractField) = string(typeof(field).name.wrapper, " located at ", show_location(field))

Base.show(io::IO, field::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} =
    print(io, "$(short_show(field))\n",
          "├── architecture: $(architecture(field))\n",
          "└── grid: $(short_show(field.grid))")

Base.show(io::IO, field::Field) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field))\n",
          "├── grid: $(short_show(field.grid))\n",
          "└── boundary conditions: $(short_show(field.boundary_conditions))")

show_status(::Nothing) = "nothing"
show_status(status) = "time=$(status.time)"

Base.show(io::IO, field::ZeroField) = print(io, "ZeroField")

short_show(array::OffsetArray{T, D, A}) where {T, D, A} = string("OffsetArray{$T, $D, $A}")

Base.show(io::IO, ::MIME"text/plain", f::AbstractField) = show(io, f)


text:
using VarianceComponentsHDFE
using Random, BenchmarkTools, Test
using JLD, SparseArrays, LinearAlgebra
using AlgebraicMultigrid, IterativeSolvers
using Laplacians

#if (BLAS.vendor() == :openblas64 && !haskey(ENV, "OPENBLAS_NUM_THREADS"))
#    blas_num_threads = min(1, Int64(round(Sys.CPU_THREADS / 2)))  # or lower?
#    println("Setting BLAS threads = $blas_num_threads")
#    BLAS.set_num_threads(blas_num_threads)
#end

pkg_dir = pkgdir(VarianceComponentsHDFE)

# Obtain environment variables
use_matlab_CMG = get(ENV, "USE_MATLAB_CMG", "false") == "true" ? true : false
run_large_benchmark = get(ENV, "VCHDFE_LARGE_BENCHMARK", "false")  == "true" ? true : false
num_rhs_large = try
        parse(Int, get(ENV, "BENCHMARK_NUM_RHS_LARGE", "2"))
    catch
        2
    end

if use_matlab_CMG
    using  Pkg
    include(string(Pkg.dir("Laplacians") , "/src/matlabSolvers.jl"))
end

# Benchmarks for the no controls case in the JL algorithm
const SUITE = BenchmarkGroup()

# Load the data
include("prepare_benchmark_data.jl")
large_data = load(pkg_dir*"/benchmark/data/large_main.jld")
X = large_data["X_GroundedLaplacian"]
X_lap = large_data["X_Laplacian"]

S_xx = large_data["S_xx"]
S_xx_sparse = sparse(S_xx) # The CMG solver can't handle a symmetric S_xx

S_xx_lap = large_data["S_xx_lap"]
S_xx_sparse_lap = sparse(S_xx_lap)
A = large_data["A"]
m,k = size(X)
k_lap = size(X_lap,2)

#SDDM solver
sol_sddm = approxchol_sddm(S_xx_sparse, verbose=false)
SUITE["Large: SDDM Solver Build for S_xx_sparse"] = @benchmarkable approxchol_sddm($S_xx_sparse, verbose=false, tol=1e-6)

sol_KMPsddm = KMPSDDMSolver(S_xx_sparse, maxits=300; verbose=false, tol=1e-6)
SUITE["Large: KMPSDDM Solver Build for S_xx_sparse"] = @benchmarkable KMPSDDMSolver($S_xx_sparse, maxits=300; verbose=false, tol=1e-6)

#Create Adjacency and LAP solver
sol_lap = approxchol_lap(A; verbose=false)
SUITE["Large: Lap Solver Build for Adjacency Matrix"] = @benchmarkable approxchol_lap($A; verbose=false)

sol_cglap = Laplacians.cgLapSolver(A,maxits=300;verbose=false)
SUITE["Large: cgLap Solver Build for Adjacency Matrix"] = @benchmarkable  Laplacians.cgLapSolver($A,maxits=300;verbose=false)

sol_KMPlap = KMPLapSolver(A,maxits=300;verbose=false,tol=1e-6)
SUITE["Large: KMPLap Solver Build for Adjacency Matrix"] = @benchmarkable  KMPLapSolver($A,maxits=300;verbose=false,tol=1e-6)

# Only a single RHS, so set p = 1
R_p = convert(Array{Float64,2}, bitrand(1,m))
rademacher!(R_p)

# Setup the JLA_RHS
JLA_RHS = (R_p*X)[1,:]
JLA_RHS_sparse =  SparseMatrixCSC{Float64,Int64}(sparse(JLA_RHS))

JLA_RHS_lap = (R_p*X_lap)[1,:]
JLA_RHS_lap_sparse = SparseMatrixCSC{Float64,Int64}(sparse(JLA_RHS_lap))

#SDDM Solvers
SUITE["Large: S_xx_sparse SDDM solver JL RHS"] = @benchmarkable sol_sddm($JLA_RHS, verbose=false)
SUITE["Large: S_xx_sparse KMPSDDM solver JL RHS"] = @benchmarkable sol_KMPsddm($JLA_RHS, verbose=false, tol=1e-6)

#LAP Solvers
SUITE["Large: S_xx_sparse LAP solver JL RHS"] = @benchmarkable sol_lap($JLA_RHS_lap, verbose=false)
SUITE["Large: S_xx_sparse cgLAP solver JL RHS"] = @benchmarkable sol_cglap($JLA_RHS_lap, verbose=false)
SUITE["Large: S_xx_sparse KMPLap solver JL RHS"] = @benchmarkable sol_KMPlap($JLA_RHS_lap, verbose=false, tol=1e-6)

# Setup and benchmark the precondiioner
SUITE["Large: S_xx precondition: AMG ruge_stuben"] = @benchmarkable aspreconditioner(ruge_stuben($S_xx))
SUITE["Large: S_xx_sparse precondition: AMG ruge_stuben"] = @benchmarkable aspreconditioner(ruge_stuben($S_xx_sparse))

P = aspreconditioner(ruge_stuben(S_xx))
P_sparse = aspreconditioner(ruge_stuben(S_xx_sparse))

# AMG/CG Benchmarks
z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx iterative solve: AMG/CG"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS, Pl = $P , log=true, maxiter=300, tol=1e-6)
z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx_sparse iterative solve: AMG/CG"] = @benchmarkable cg!($z, $S_xx_sparse, $JLA_RHS, Pl = $P , log=true, maxiter=300, tol=1e-6)

z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx iterative solve RHS_sparse: AMG/CG"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS_sparse, Pl = $P , log=true, maxiter=300, tol=1e-6)
z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx_sparse iterative solve RHS_sparse: AMG/CG"] = @benchmarkable cg!($z, $S_xx_sparse, $JLA_RHS_sparse, Pl = $P , log=true, maxiter=300, tol=1e-6)

z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx iterative solve: AMG/CG, P_sparse"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS, Pl = $P_sparse, log=true, maxiter=300, tol=1e-6)
z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx_sparse iterative solve: AMG/CG, P_sparse"] = @benchmarkable cg!($z, $S_xx_sparse, $JLA_RHS, Pl = $P_sparse, log=true, maxiter=300, tol=1e-6)

z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx iterative solve RHS_sparse: AMG/CG, P_sparse"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS_sparse, Pl = $P_sparse, log=true, maxiter=300, tol=1e-6)
z = 0.1.*ones(length(JLA_RHS))
SUITE["Large: S_xx_sparse iterative solve RHS_sparse: AMG/CG, P_sparse"] = @benchmarkable cg!($z, $S_xx_sparse, $JLA_RHS_sparse, Pl = $P_sparse, log=true, maxiter=300, tol=1e-6)

if use_matlab_CMG
    SUITE["Large: S_xx iterative solve: CMG"] = @benchmarkable matlabCmgSolver($S_xx_sparse, $JLA_RHS; tol=1e-6, maxits=300)
    SUITE["Large: S_xx iterative solve: CMG, sparse RHS"] = @benchmarkable matlabCmgSolver($S_xx_sparse, $JLA_RHS_sparse; tol=1e-6, maxits=300)
end


text:
# "smart" summary and REPL printing

# Generic print for power system components (a new type that includes generation
# technology types). The long version simply prints every field of the
# device. Tailored versions can be made for specific devices, if/when
# desired. JJS 11/15/18
function printPST(pst::Component, short = false,
                  io::IO = stdout)
    if short
        pst_summary = Base.summary(pst)
        # objects of type TechnicalParams do not have a 'name' field
        if in(:name, fieldnames(typeof(pst)))
            print(io, "$pst_summary(name=\"$(pst.name)\")")
        else
            print(io, "$pst_summary")
        end
    else
        print(io, Base.summary(pst),":")
        fnames = fieldnames(typeof(pst))
        for i = 1:nfields(fnames)
            thefield = getfield(pst, i)
            # avoid warning about printing nothing
            if thefield==nothing
                print(io, "\n   ", fnames[i], ": nothing")
            else
                print(io, "\n   ", fnames[i], ": ", thefield)
            end
        end
    end
end
# overload Base.show with printPST function defined above
# single-line format
Base.show(io::IO, pst::Component) = printPST(pst, true, io)
# Multi-line format for plaintext (e.g. from repl); can specify for HTML and
# others too
Base.show(io::IO, ::MIME"text/plain", pst::Component) =
    printPST(pst, false, io)

# provide limited additional information to summary of GenClasses
function summaryGenClasses(gens::GenClasses)
    ntherm = gens.thermal isa Nothing ? 0 : length(gens.thermal)
    nrenew = gens.renewable isa Nothing ? 0 : length(gens.renewable)
    nhydro = gens.hydro isa Nothing ? 0 : length(gens.hydro)
    return "GenClasses(T:$ntherm,R:$nrenew,H:$nhydro)"
end
Base.summary(gens::GenClasses) = summaryGenClasses(gens::GenClasses)

# improve output for System

# overload Base.summary for System to be useful, e.g., in output of
# varinfo(); also used now in printPowerSystem, JJS 11/26/18
function summaryPowerSystem(system::System)
    busstr = string("buses:", length(system.buses))
    genstr = Base.summary(system.generators)
    loadstr = string("loads:", length(system.loads))
    branchstr = system.branches isa Nothing ? "nothing" :
        string("branches:", length(system.branches))
    # is storage an array? not sure yet, JJS 11/28/18
    storstr = system.storage isa Nothing ? "nothing" :
        string("storage:", length(system.storage))
    # include some output for basepower or time_periods? currently not doing so
    output = string("System($busstr,$genstr,$loadstr,$branchstr,$storstr)")
    return output
end
Base.summary(system::System) = summaryPowerSystem(system::System)

# print function for System type
function printPowerSystem(system::System, short = false,
                  io::IO = stdout)
    if short
        print(io, Base.summary(system))
    else
        print(io, "System:")
        fnames = fieldnames(typeof(system))
        for i = 1:nfields(fnames)
            fname = fnames[i]
            thefield = getfield(system, i)
            if thefield==nothing
                print(io, "\n   ", fname, ": nothing")
            elseif fname==:generators
                # print out long version for generators
                print(io, "\n   ", fname, ": \n     ")
                printPST(system.generators, false, io)
                print(io, "\n     (end generators)")
            else
                print(io, "\n   ", fname, ": ", thefield)
            end
        end
    end
end
Base.show(io::IO, system::System) = printPowerSystem(system, true, io)
Base.show(io::IO, ::MIME"text/plain", system::System) =
    printPowerSystem(system, false, io)


text:
include("Graph.jl")
include("NodeCentrality.jl")

using SparseArrays
using LinearAlgebra

function spA(G)
    Is = zeros(Int, 2*G.m)
    Js = zeros(Int, 2*G.m)
    Vs = zeros(2*G.m)
    for (ID, u, v, w) in G.E
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = w
        Is[G.m + ID] = v
        Js[G.m + ID] = u
        Vs[G.m + ID] = w
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function EdgeKatz(G; fac1 = 0.85)
    KC = KatzCentrality(G, fac=fac1)
    EK = zeros(G.m)
    for (ID, u, v, w) in G.E
        EK[ID] = (1.0 + KC[u]) * (1.0 + KC[v])
    end
    return EK
end

function getLp(G)
    L = zeros(G.n, G.n)
    for (ID, u, v, w) in G.E
        L[u, v] -= w
        L[v, u] -= w
        L[u, u] += w
        L[v, v] += w
    end
    L .+= (1.0 / G.n)
    Lp = inv(L)
    Lp .-= (1.0 / G.n)
    return Lp
end

function BiharmonicDistanceRelatedCentrality(G)
    Lp = getLp(G)
    C = zeros(G.m)
    for (ID, u, v, w) in G.E
        C[ID] = norm(Lp[:, u] - Lp[:, v])^2
    end
    return C
end

function SpanningEdgeCentrality(G)
    Lp = getLp(G)
    C = zeros(G.m)
    for (ID, u, v, w) in G.E
        C[ID] = Lp[u, u] + Lp[v, v] - 2*Lp[u, v]
    end
    return C
end

function EdgeBetweenness(G)
    g = getAdjacentList(G)
    idx = Dict{Tuple{Int32, Int32}, Int32}()
    for (ID, u, v, w) in G.E
        idx[(u, v)] = ID
        idx[(v, u)] = ID
    end
    Cb = zeros(G.m)
    p = Array{Array{Int32, 1}, 1}(undef, G.n)
    d = zeros(Int32, G.n)
    S = zeros(Int32, G.n+10)
    sigma = zeros(G.n)
    d = zeros(Int32, G.n)
    Q = zeros(Int32, G.n+10)
    delta = zeros(G.n)

    for s = 1 : G.n
        foreach(i -> p[i] = [], 1 : G.n)
        top = 0
        sigma .= 0
        sigma[s] = 1.0
        d .= -1
        d[s] = 0
        front = 1
        rear = 1
        Q[1] = s

        while front <= rear
            v = Q[front]
            front += 1
            top += 1
            S[top] = v
            for w in g[v]
                if d[w] < 0
                    rear += 1
                    Q[rear] = w
                    d[w] = d[v] + 1
                end
                if d[w] == (d[v] + 1)
                    sigma[w] += sigma[v]
                    push!(p[w], v)
                end
            end
        end

        delta .= 0

        while top > 0
            w = S[top]
            top -= 1
            for v in p[w]
                delta[v] += ((sigma[v] / sigma[w]) * (1 + delta[w]))
                Cb[idx[(v, w)]] += ((sigma[v] / sigma[w]) * (1 + delta[w]))
            end
        end

    end

    return Cb
end

function ApproxEdgeBetweenness(G; eps = 0.5)
    g = getAdjacentList(G)
    idx = Dict{Tuple{Int32, Int32}, Int32}()
    for (ID, u, v, w) in G.E
        idx[(u, v)] = ID
        idx[(v, u)] = ID
    end
    pivots = selectPivots(G, eps)
    Cb = zeros(G.m)
    p = Array{Array{Int32, 1}, 1}(undef, G.n)
    d = zeros(Int32, G.n)
    S = zeros(Int32, G.n+10)
    sigma = zeros(G.n)
    d = zeros(Int32, G.n)
    Q = zeros(Int32, G.n+10)
    delta = zeros(G.n)

    for s in pivots
        foreach(i -> p[i] = [], 1 : G.n)
        top = 0
        sigma .= 0
        sigma[s] = 1.0
        d .= -1
        d[s] = 0
        front = 1
        rear = 1
        Q[1] = s

        while front <= rear
            v = Q[front]
            front += 1
            top += 1
            S[top] = v
            for w in g[v]
                if d[w] < 0
                    rear += 1
                    Q[rear] = w
                    d[w] = d[v] + 1
                end
                if d[w] == (d[v] + 1)
                    sigma[w] += sigma[v]
                    push!(p[w], v)
                end
            end
        end

        delta .= 0

        while top > 0
            w = S[top]
            top -= 1
            for v in p[w]
                delta[v] += ((sigma[v] / sigma[w]) * (1 + delta[w]))
                Cb[idx[(v, w)]] += ((sigma[v] / sigma[w]) * (1 + delta[w]))
            end
        end

    end

    return Cb
end

function getW(G)
    IpL = zeros(G.n, G.n)
    for (ID, u, v, w) in G.E
        IpL[u, v] -= w
        IpL[v, u] -= w
        IpL[u, u] += w
        IpL[v, v] += w
    end
    foreach(i -> IpL[i, i] += 1, 1 : G.n)
    return inv(IpL)
end

function ForestCentrality(G)
    W = getW(G)
    C = zeros(G.m)
    for (ID, u, v, w) in G.E
        C[ID] = (W[u, u] + W[v, v] - 2*W[u, v]) / W[u, v]
    end
    return C
end


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

mutable struct IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus <: SwaggerModel
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1FlowSchemaCondition} } # spec name: conditions

    function IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus(;conditions=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus

const _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus = Dict{Symbol,Symbol}(Symbol("conditions")=>Symbol("conditions"))
const _property_types_IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus = Dict{Symbol,String}(Symbol("conditions")=>"Vector{IoK8sApiFlowcontrolV1alpha1FlowSchemaCondition}")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus)
    true
end

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