text:
using ClinicalTrialSampleSize
using Base.Test

# TwoSampleMean
# Constructor: TwoSampleMean()
@test_throws ErrorException TwoSampleMean(-Inf, Inf, 1)
@test_throws ErrorException TwoSampleMean(25, 20, 0)

params = (
    (TwoSampleMean(3, 2, 1, true), 16.722, 0.8, 1, 0.05, "two"),
    (TwoSampleMean(10, 7, 2), 43.780, 0.9, 5, 0.05, "two"),
    (TwoSampleMean(7, 4, 1, true), 13.842, 0.9, 3, 0.1, "one"),
    (TwoSampleMean(7, 4, 1), 13.139, 0.9, 3, 0.1, "one"),
    (TwoSampleMean(50, 45, 1.5), 137.123, 0.85, (10.4, 15.2), 0.01, "one"),
    (TwoSampleMean(50, 45, 1.5, true), 138.142, 0.85, (10.4, 15.2), 0.01, "one"),
    (TwoSampleMean(120, 114, 1), 171.548, 0.95, (13, 17.5), 0.05, "two"),
    (TwoSampleMean(120, 114, 1, true), 172.634, 0.95, (13, 17.5), 0.05, "two"),
)

for (test, n, p, std, alpha, side) in params

        p_est = power(test; n = n, std = std, alpha = alpha, side = side)
        n_est = sample_size(test; power = p, std = std, alpha = alpha, side = side)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# Two-sample non-inferiority/superiority for mean
# Constructor: TwoSampleMeanInferior/Superior()
@test_throws ErrorException TwoSampleMeanInferior(-Inf, 5, 1, -1)
@test_throws ErrorException TwoSampleMeanInferior(30, 35, -3, -3)
@test_throws ErrorException TwoSampleMeanInferior(30, 35, 1, 3)
@test_throws ErrorException TwoSampleMeanSuperior(-Inf, Inf, 1, 1)
@test_throws ErrorException TwoSampleMeanSuperior(35, 30, -2, 5)
@test_throws ErrorException TwoSampleMeanSuperior(35, 30, 1, -3)

params = (
    (TwoSampleMeanInferior(2, 2, 1, -1, true), 13.170, 0.8, 1, 0.05),
    (TwoSampleMeanInferior(1.9, 2.5, 2, -1), 16.850, 0.85, 0.5, 0.05),
    (TwoSampleMeanInferior(1.94, 2.0, 1, -2.3), 46.286, 0.95, (1.1, 4.5), 0.05),
    (TwoSampleMeanInferior(1.94, 2.0, 1, -2.3, true), 47.227, 0.95, (1.1, 4.5), 0.05),
    (TwoSampleMeanSuperior(2.5, 2, 2, 0), 45.710, 0.8, 1.3, 0.1),
    (TwoSampleMeanSuperior(20, 15, 1, 3, true), 81.810, 0.85, 5.5, 0.1),
    (TwoSampleMeanSuperior(121, 119, 1.5, 2.3), 8278.386, 0.95, (3, 9), 0.1),
    (TwoSampleMeanSuperior(121, 119, 1.5, 2.3, true), 8278.324, 0.95, (3, 9), 0.1),
)

for (test, n, p, std, alpha) in params

        p_est = power(test; n = n, std = std, alpha = alpha)
        n_est = sample_size(test; power = p, std = std, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# Two-sample equivalence for mean
# Constructor: TwoSampleMeanEqual()
@test_throws ErrorException TwoSampleMeanEqual(-Inf, Inf, 1, 1)
@test_throws ErrorException TwoSampleMeanEqual(30, 30, -1, 0)
@test_throws ErrorException TwoSampleMeanEqual(30, 30, 1, -5)

params = (
    (TwoSampleMeanEqual(2, 2, 1, 0.3), 190.308, 0.8, 1, 0.05),
    (TwoSampleMeanEqual(1.94, 1.95, 2, 0, true), 8919.364, 0.85, 0.25, 0.05),
    (TwoSampleMeanEqual(112, 110, 1, 1.7), 15542.81, 0.95, (5.2, 10.3), 0.1),
    (TwoSampleMeanEqual(59, 54, 2, 3, true), 23.843, 0.8, (4, 2.5), 0.1),
)

for (test, n, p, std, alpha) in params

        p_est = power(test; n = n, std = std, alpha = alpha)
        n_est = sample_size(test; power = p, std = std, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# Two-sample proportion
# Constructor: TwoSampleProp()
@test_throws ErrorException TwoSampleProp(-1.2, 1.2, 1)
@test_throws ErrorException TwoSampleProp(0.3, 0.25, -1)

params = (
    (TwoSampleProp(0.5, 0.2, 1), 35.756, 0.8, 0.05, "two"),
    (TwoSampleProp(0.9, 0.4, 2), 10.235, 0.85, 0.05, "two"),
    (TwoSampleProp(0.44, 0.4, 1), 2603.41, 0.95, 0.1, "one"),
)

for (test, n, p, alpha, side) in params

        p_est = power(test; n = n, alpha = alpha, side = side)
        n_est = sample_size(test; power = p, alpha = alpha, side = side)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# Two-sample non-inferiority/superiority proportion
# Constructor: TwoSamplePropInferior/Superior
@test_throws ErrorException TwoSamplePropInferior(-1.2, 1.2, 1, 0)
@test_throws ErrorException TwoSamplePropInferior(0.3, 0.2, -2, -0.3)
@test_throws ErrorException TwoSamplePropInferior(0.3, 0.2, 1, 0.1)
@test_throws ErrorException TwoSamplePropSuperior(1.1, 0.3, 1, 0)
@test_throws ErrorException TwoSamplePropSuperior(0.4, 0.2, -1, 0.2)
@test_throws ErrorException TwoSamplePropSuperior(0.4, 0.2, 1, -0.2)

params = (
    (TwoSamplePropInferior(0.25, 0.3, 1, -0.06), 24575.66, 0.8, 0.05),
    (TwoSamplePropInferior(0.71, 0.8, 2, -0.05), 883.028, 0.85, 0.1),
    (TwoSamplePropSuperior(0.7, 0.5, 1, 0.1), 393.937, 0.9, 0.05),
    (TwoSamplePropSuperior(0.32, 0.27, 1.5, 0.01), 1149.05, 0.85, 0.1),
)

for (test, n, p, alpha) in params

        p_est = power(test; n = n, alpha = alpha)
        n_est = sample_size(test; power = p, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# Two-sample equivalence proportion
# Constructor: TwoSamplePropEqual
@test_throws ErrorException TwoSamplePropEqual(1.9, 1.2, 1, 0.1)
@test_throws ErrorException TwoSamplePropEqual(0.2, 0.2, -1, 0)
@test_throws ErrorException TwoSamplePropEqual(0.2, 0.2, 1, -2)

params = (
    (TwoSamplePropEqual(0.35, 0.3, 1, 0.06), 37466.83, 0.8, 0.05),
    (TwoSamplePropEqual(0.85, 0.8, 1, 0.1), 851.494, 0.85, 0.1),
    (TwoSamplePropEqual(0.4, 0.4, 2, 0.02), 9456.681, 0.95, 0.1),
)

for (test, n, p, alpha) in params

        p_est = power(test; n = n, alpha = alpha)
        n_est = sample_size(test; power = p, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# McNemar's paired test
# Constructor: McNemarProp
@test_throws ErrorException McNemarProp(1.9, 0.3)

params = (
    (McNemarProp(0.35, 0.5), 294.144, 0.8, 0.05, "two"),
    (McNemarProp(0.85, 0.8), 4742.158, 0.85, 0.1, "two"),
    (McNemarProp(0.4, 0.2), 123.608, 0.95, 0.1, "one"),
)

for (test, n, p, alpha, side) in params

        p_est = power(test; n = n, alpha = alpha, side = side)
        n_est = sample_size(test; power = p, alpha = alpha, side = side)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# Two sample crossover test for mean
# Constructor: CrossoverMean
@test_throws ErrorException CrossoverMean(-Inf, Inf)

params = (
    (CrossoverMean(4.35, 3.5, true), 136.774, 0.8, 5, 0.05, "two"),
    (CrossoverMean(0.85, 0.8), 23487.195, 0.85, 3, 0.01, "two"),
    (CrossoverMean(10.4, 10.2, true), 8.127, 0.9, 0.3, 0.1, "one"),
    (CrossoverMean(1.4, 0.2), 11.894, 0.95, 2, 0.1, "one"),
)

for (test, n, p, std, alpha, side) in params

        p_est = power(test; n = n, std = std, alpha = alpha, side = side)
        n_est = sample_size(test; power = p, std = std, alpha = alpha, side = side)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

return true


text:
_hashpowers(u::UInt) = u
function _hashpowers(u::UInt, power::Tuple, powers...)
    if iszero(power[2])
        _hashpowers(u, powers...)
    else
        _hashpowers(hash(power, u), powers...)
    end
end
function Base.hash(m::AbstractMonomial, u::UInt)
    nnz = count(!iszero, exponents(m))
    if iszero(nnz)
        hash(1, u)
    elseif isone(nnz) && isone(degree(m))
        hash(variable(m), u)
    else
        _hashpowers(u, powers(m)...)
    end
end

function Base.hash(t::AbstractTerm, u::UInt)
    if iszero(t)
        hash(0, u)
    elseif coefficient(t) == 1
        hash(monomial(t), u)
    else
        hash(monomial(t), hash(coefficient(t), u))
    end
end

function Base.hash(p::AbstractPolynomial, u::UInt)
    if iszero(p)
        hash(0, u)
    else
        reduce((u, t) -> hash(t, u), terms(p); init=u)
    end
end

ext: jl


text:
#broad_pdf.jl
using JSON

function create_broad_pdf(n_x::Int, n_y::Int, n_z::Int, dependence::Float64=0.) :: Dict{Tuple{Int,Int,Int},BigFloat}
    @assert n_x ≥ 2 "create_broad_pdf(): n_x >= 2 needed"
    @assert n_y ≥ 2 "create_broad_pdf(): n_y >= 2 needed"
    @assert n_z ≥ 2 "create_broad_pdf(): n_z >= 2 needed"
    @assert 0. ≤ dependence ≤ 1. "create_broad_pdf(): indep in [0,1] needed"
    pdf = Dict{Tuple{Int,Int,Int},BigFloat}()
    for y = 1:n_y
        for z = 1:n_z
            for x = 1:n_x
                pdf[ (x,y,z) ] = ( 1 + dependence*(rand()-.5) )/( n_x*n_y*n_z )
            end # for z
        end # for z
    end # for y
    return pdf
end # create_broad_pdf()

function broad_pdf_to_file(filename::String, n_x::Int, n_y::Int, n_z::Int, dependence::Float64=0.) :: Void
    open(filename,"w") do file
        pdf = create_broad_pdf(n_x,n_y,n_z,dependence)
        pdflist = []
        for k in keys(pdf)
            (x,y,z) = k
            e = [[x,y,z],pdf[(x,y,z)]]
            push!(pdflist, e )
        end
        JSON.print(file,pdflist)
    end
    return nothing
end # broad_pdf_to_file()

function loop_broad_pdfs_to_files() :: Void
    for n_x = 2:5
        for n_y = n_x : 5*n_x
            n_z = n_y
            filename = "broad_indep--"*string(n_x)*"-"*string(n_y)*"-"*string(n_z)*".dens"
            broad_pdf_to_file(filename,n_x,n_y,n_z)
            filename = "broad_dep--"*string(n_x)*"-"*string(n_y)*"-"*string(n_z)*".dens"
            broad_pdf_to_file(filename,n_x,n_y,n_z,1.)
        end
    end
    return nothing;
end

;
#EOF

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
using CellListMap
using StaticArrays
import Random

#
# In this test we compute the "gravitational force", assigning to each particle
# a different mass. In this case, the closure is used to pass the masses and
# the force vector to the function that computes the potential.
#
function gravitational_force(;N=100_000,parallel=true,x=nothing)

    # Number of particles, sides and cutoff
    sides = @SVector [250,250,250]
    cutoff = 10.
    box = Box(sides, cutoff)

    # Particle positions
    Random.seed!(321)
    if x === nothing
        x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    end

    # masses
    mass = [ 5 * x[i][1] for i in 1:N ]

    # Initialize auxiliary linked lists
    cl = CellList(x, box, parallel=parallel)

    # Function to be evalulated for each pair
    function calc_forces!(x, y, i, j, d2, mass, forces)
        G = 9.8 * mass[i] * mass[j] / d2
        d = sqrt(d2)
        df = (G / d) * (x - y)
        forces[i] = forces[i] - df
        forces[j] = forces[j] + df
        return forces
    end

    # Preallocate and initialize forces
    forces = [ zeros(SVector{3,Float64}) for i in 1:N ]

    # Run pairwise computation
    map_pairwise!(
        (x, y, i, j, d2, forces) -> calc_forces!(x, y, i, j, d2, mass, forces),
        forces,box,cl,
        parallel=parallel
    )

    return forces
end


text:
using NumericalIntegration
using LinearAlgebra
using FFTW

#=
This function computes the Fourier modes either integrating "by hand" or using FFTW. We split the comlex modes, associated to the complex exponentials into a and b, associated to the cos and sin functions by a(kx,ky) = 2*Re(f(kx,ky)) and b(kx,ky) = -2*Im(f(kx,ky)).
=#

function Fourier_modes_2D(e :: Array{T,2}, x :: Array{T,1}, y :: Array{T,1}, order_x :: Integer, order_y :: Integer, integration_method :: String) where {T<:Real}
    Lx, Ly = x[end]-x[1], y[end]-y[1]
    xmid, ymid = (x[end]+x[1])/2, (y[end]+y[1])/2
    nx, ny = length(x),length(y)

    if order_x > Integer(floor(nx)) || order_y > Integer(floor(ny))
        print("\n Error, requested order is higher than the number of grid points\n")
        nothing
    end

    if integration_method == "trapezoidal"

        a, b = zeros(order_x+1,order_y+1), zeros(order_x+1,order_y+1)
        xx, yy = zeros(nx,ny), zeros(nx,ny)
        #We create the 2D grid
        for i in 1:ny
            xx[:,i] = x
        end
        for i in 1:nx
            yy[i,:] = y
        end

        for my in 0:order_y
            for mx in 0:order_x
                integral_cos, integral_sin = zeros(nx), zeros(nx)
                integrand_cos = 2/(Lx*Ly).*e.*cos.(2*π*(mx/Lx.*(xx.-xmid).+my/Ly.*(yy.-ymid)))
                integrand_sin = 2/(Lx*Ly).*e.*sin.(2*π*(mx/Lx.*(xx.-xmid).+my/Ly.*(yy.-ymid)))
                for i in 1:nx
                    integral_cos[i] = integrate(y,integrand_cos[i,:], Trapezoidal())
                    integral_sin[i] = integrate(y,integrand_sin[i,:], Trapezoidal())
                end
                a[mx+1,my+1] = integrate(x,integral_cos, Trapezoidal())
                b[mx+1,my+1] = integrate(x,integral_sin, Trapezoidal())
            end
        end

    elseif integration_method == "FFTW"

        plan = plan_rfft(e)
        fft_z = (1/(nx*ny).*(plan*e)[1:order_x+1,1:order_y+1])
        a = 2 .*real(fft_z)
        b = -2 .*imag(fft_z)

    else
        print("\n No such integration method\n")
        nothing
    end
    return a, b
end


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

@doc raw"""Eviction evicts a pod from its node subject to certain policies and safety constraints. This is a subresource of Pod.  A request to cause such an eviction is created by POSTing to .../pods/<pod name>/evictions.

    IoK8sApiPolicyV1beta1Eviction(;
        apiVersion=nothing,
        deleteOptions=nothing,
        kind=nothing,
        metadata=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - deleteOptions::IoK8sApimachineryPkgApisMetaV1DeleteOptions : DeleteOptions may be provided
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : ObjectMeta describes the pod that is being evicted.
"""
mutable struct IoK8sApiPolicyV1beta1Eviction <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    deleteOptions::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1DeleteOptions } # spec name: deleteOptions
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata

    function IoK8sApiPolicyV1beta1Eviction(;apiVersion=nothing, deleteOptions=nothing, kind=nothing, metadata=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("deleteOptions"), deleteOptions)
        setfield!(o, Symbol("deleteOptions"), deleteOptions)
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        o
    end
end # type IoK8sApiPolicyV1beta1Eviction

const _property_map_IoK8sApiPolicyV1beta1Eviction = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("deleteOptions")=>Symbol("deleteOptions"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"))
const _property_types_IoK8sApiPolicyV1beta1Eviction = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("deleteOptions")=>"IoK8sApimachineryPkgApisMetaV1DeleteOptions", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1Eviction }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1Eviction))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1Eviction }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiPolicyV1beta1Eviction[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1Eviction }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1Eviction[property_name]

function check_required(o::IoK8sApiPolicyV1beta1Eviction)
    true
end

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


text:
# # Battle Royale
#
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../battle.mp4" type="video/mp4">
# </video>
# ```
#
# This example illustrates how to leverage higher dimensions of a `GridSpace` to identify
# the distance from neighbors not just spatially, but also categorically. We'll also use
# the [`walk!`](@ref) function extensively.

# The `Models` module includes this example as [`Models.battle`](@ref).

# ## Rules of Engagement
# Agents wander around the map looking for opponents. When a grid space is occupied by
# two or more agents there will be a battle. With experience gained from the fight, the
# victor searches for more opponents to crush and losers scurry away defeated or possibly
# even die. This process repeats until there is a single, definitive winner.
#
# For this battle ground to exist, the following rules must be followed:
# - Agents have an experience level, starting at level 1 up to a maximum of 10.
# - Agents will search for the nearest worthy opponent (one with equal or ±1 experience level) and move towards them to attack, so long as something more important doesn't happen, which could be
#   - A tougher opponent (with experience level +2 or higher) is nearby: run!
#   - There are no worthy opponents available, but there are weak ones (with experience level -2 or lower): chase them down.
#   - Capture and taunt a weaker opponent, then kill them.
#   - Notice a tough opponent is occupied, sneak up and kill them.
#   - There is no-one worthy to fight, but also no-one left to taunt. All bets are off: THERE CAN BE ONLY ONE.

# Battles are won by weighted chance - a higher level gives an agent a larger chance of
# winning, but does not guarantee it. When a victor is chosen
#   - The difference in experience between opponents is swapped.
#   - If an agents experience reaches 0, they die.
#
# Captured opponents will be killed once taunted. The captor will gain half of their
# experience. If an opportunist manages to take the captor by surprise, they can gain
# up to half of the captor's experience. This means a level 1 agent may eliminate a
# level 10 captor and jump straight to level 6.
#
# Once all rules of engagement have been exhausted, the final showdown begins. Opponents
# fight their closest adversary regardless of experience level. Winner takes all.

# ## Model Setup

using Random # hide
using Agents
using InteractiveDynamics
using CairoMakie
CairoMakie.activate!() # hide

mutable struct Fighter <: AbstractAgent
    id::Int
    pos::Dims{3}
    has_prisoner::Bool
    capture_time::Int
    shape::Symbol # For plotting
end

# As you can see, the properties of out agent are very simple and contain only two
# parameters that are needed to store context from one time step to the next. All
# other properties needed are stored in the space. `pos` is three-dimensional, two
# for the actual space agents move within, and a third categorical dimension representing
# their level. `shape` is used solely for plotting (well, used once just for convenience).

# Now let's set up the battle field:
function battle(; fighters = 50)
    model = ABM(
        Fighter,
        GridSpace((100, 100, 10); periodic = false);
        scheduler = Schedulers.randomly,
    )

    n = 0
    while n != fighters
        pos = (rand(model.rng, 1:100, 2)..., 1) # Start at level 1
        if isempty(pos, model)
            add_agent!(pos, model, false, 0, :diamond)
            n += 1
        end
    end

    return model
end

Random.seed!(6547) # hide
model = battle()

# 50 opponents positioned randomly on a 100x100 grid, with no escape
# (`periodic = false`). To leverage categorical dimensions fully, non-periodic chebyshev
# space is necessary.

# ## Game Dynamics
#
# To implement the rules of engagement, only an `agent_step!` function is required,
# along with a few helper functions.

space(agent) = agent.pos[1:2]
level(agent) = agent.pos[3]

# `space` allows us to invoke a number of helpful utilities provided by Agents.jl
# but only operate on our spatial dimensions, `level` is a wrapper to access the agent's
# experience easily.

# Nearest agents that satisfy our search criteria can be identified via Euclidean
# distance solely on the spatial dimensions of our `GridSpace`.

function closest_target(agent::Fighter, ids::Vector{Int}, model::ABM)
    if length(ids) == 1
        closest = ids[1]
    else
        close_id = argmin(map(id -> edistance(space(agent), space(model[id]), model), ids))
        closest = ids[close_id]
    end
    return model[closest]
end

# Since our battles are only between opponents with equal, or as much as one level apart,
# the odds can be set explicitly. Stronger opponents have twice the capacity of winning a
# match.

function battle!(one::Fighter, two::Fighter, model)
    if level(one) == level(two)
        ## Odds are equivalent
        one_winner = rand(model.rng) < 0.5
    elseif level(one) > level(two)
        ## Odds are in favor of one
        one_winner = 2 * rand(model.rng) > rand(model.rng)
    else
        ## Odds are in favor of two
        one_winner = rand(model.rng) > 2 * rand(model.rng)
    end

    one_winner ? (up = one; down = two) : (up = two; down = one)

    new_lvl_up = min(level(up) + 1, 10)
    new_pos_up =
        clamp.(rand(model.rng, -1:1, 2) .+ space(up), [1, 1], size(model.space)[1:2])
    move_agent!(up, (new_pos_up..., new_lvl_up), model)
    new_lvl_down = level(down) - 1
    if new_lvl_down == 0
        kill_agent!(down, model)
    else
        move_agent!(down, (space(down)..., new_lvl_down), model)
    end
end

# If an agent has a prisoner, it will taunt it for a time, then kill it, so long as
# an opportunist doesn't sneak up on them first!
# Here we use the tuple constructor with [`nearby_ids`](@ref) to look for agents at
# the same position as the captor `(0, 0)`, and *any* level `(..., 10)`. We could also
# use the range constructor in this instance
# `nearby_ids(agent, model, [(1, 0:0), (2, 0:0)])`, meaning which is more performant but not as
# readable.

function captor_behavior!(agent, model)
    close_ids = collect(nearby_ids(agent, model, (0, 0, 10)))
    if length(close_ids) == 1
        ## Taunt prisoner or kill it
        prisoner = model[close_ids[1]]
        if prisoner.capture_time > 10
            agent.shape = :rect
            gain = ceil(Int, level(prisoner) / 2)
            new_lvl = min(level(agent) + gain, 10)
            kill_agent!(prisoner, model)
            agent.has_prisoner = false
            move_agent!(agent, (space(agent)..., new_lvl), model)
        end
    else
        ## Someone is here to kill the captor. Could be more than one opponent
        prisoner = [model[id] for id in close_ids if model[id].capture_time > 0][1]
        exploiter = rand(
            model.rng,
            [
                model[id]
                for
                id in close_ids if
                model[id].capture_time == 0 && model[id].has_prisoner == false
            ],
        )
        exploiter.shape = :rect
        gain = ceil(Int, level(agent) / 2)
        new_lvl = min(level(agent) + rand(model.rng, 1:gain), 10)
        kill_agent!(agent, model)
        move_agent!(exploiter, (space(exploiter)..., new_lvl), model)
        ## Prisoner runs away in the commotion
        prisoner.shape = :utriangle
        prisoner.capture_time = 0
        walk!(prisoner, (rand(model.rng, -1:1, 2)..., 0), model)
    end
end

# When there are only few fighters standing, the stakes are higher. Prior experience is
# paramount since there is no gain, and fights are to the death.

function endgame!(agent, model)
    origin = space(agent)
    end_ids = collect(Iterators.filter(
        id -> model[id].shape == :circle && id != agent.id,
        allids(model),
    ))
    agent.shape = :circle
    if !isempty(end_ids)
        opponent = closest_target(agent, end_ids, model)
        target = space(opponent)
        if origin == target
            ## Battle
            agent.shape = :rect
            opponent.shape = :rect
            showdown!(agent, opponent, model)
        else
            walk!(agent, (sign.(target .- origin)..., 0), model)
        end
    end
end

function showdown!(one::Fighter, two::Fighter, model)
    if level(one) == level(two)
        ## Odds are equivalent
        one_winner = rand(model.rng) < 0.5
    elseif level(one) > level(two)
        ## Odds are in favor of one
        one_winner = level(one) - level(two) * rand(model.rng) > rand(model.rng)
    else
        ## Odds are in favor of two
        one_winner = rand(model.rng) > level(two) - level(one) * rand(model.rng)
    end

    one_winner ? kill_agent!(two, model) : kill_agent!(one, model)
end

# The rest of our interactions flow down a hierarchy, so we'll place them directly in the
# `agent_step!` function. We use the tuple search for `occupied_ids` here, as we did with
# `close_ids` above. The rest of the searches however use the range search to provide a
# more precise criteria.
#
# The easiest context to explore is `worthy_ids`: all we want to do is find an agent with
# a similar level. If we used the tuple search here, we would have to search
# `(100, 100, 1)` - even though we are not at all interested in the spatial location of
# the neighbors at this time. `(3, -1:1)` is therefore more accurate representation.
#
# A more complex example is that of `strong_ids`. We are looking for agents with a level
# 2-4 points higher withing a distance of `(5, 5)`. The range search becomes a little
# verbose, but precise. An equivalent tuple search is not completely possible however.
# The closest solution is `(5, 5, 4)`, which also looks for *weaker* opponents and must
# be filtered to the correct neighbor set after the fact. In this instance the range
# search has significant performance gains.

function agent_step!(agent, model)
    if agent.capture_time > 0
        ## Captured agents are powerless, but we need to keep track of how long
        ## they have been in this state
        agent.capture_time += 1
    elseif agent.has_prisoner
        captor_behavior!(agent, model)
    else
        origin = space(agent)
        ## Find agents that have captives, they are not focused
        occupied_ids = collect(Iterators.filter(
            id -> model[id].has_prisoner,
            nearby_ids(agent, model, (7, 7, 10)),
        ))
        if !isempty(occupied_ids)
            ## Sneak up behind them
            target = space(closest_target(agent, occupied_ids, model))
            agent.shape = :pentagon
            walk!(agent, (sign.(target .- origin)..., 0), model)
        else
            ## Opponents that are greatly higher in rank that the current agent
            strong_ids = collect(nearby_ids(agent, model, [(1, -5:5), (2, -5:5), (3, 2:4)]))
            if !isempty(strong_ids)
                ## Run away from nearest
                target = space(closest_target(agent, strong_ids, model))
                agent.shape = :utriangle
                walk!(agent, (sign.(origin .- target)..., 0), model)
            else
                ## There are no distractions. Search for the closest worthy opponent
                worthy_ids = collect(nearby_ids(agent, model, [(3, -1:1)]))
                if !isempty(worthy_ids)
                    opponent = closest_target(agent, worthy_ids, model)
                    target = space(opponent)
                    if origin == target
                        ## Battle
                        agent.shape = :rect
                        opponent.shape = :rect
                        battle!(agent, opponent, model)
                    else
                        ## Move towards worthy opponent
                        agent.shape = :diamond
                        walk!(agent, (sign.(target .- origin)..., 0), model)
                    end
                else
                    ## Find any weak targets in the vicinity
                    weak_ids = collect(nearby_ids(
                        agent,
                        model,
                        [(1, -10:10), (2, -10:10), (3, -4:-2)],
                    ))
                    if !isempty(weak_ids)
                        prisoner = closest_target(agent, weak_ids, model)
                        target = space(prisoner)
                        if origin == target
                            ## Capture and taunt target
                            agent.has_prisoner = true
                            agent.shape = :vline
                            prisoner.capture_time += 1
                            prisoner.shape = :hline
                        else
                            ## Chase down nearest (can move 2 steps at a time!)
                            agent.shape = :star4
                            walk!(agent, (2 .* sign.(target .- origin)..., 0), model)
                        end
                    else
                        ## Abandon honour. This is the end
                        endgame!(agent, model)
                    end
                end
            end
        end
    end

    return nothing
end

# ## Let the Battle Begin
# Plotting is relatively straightforward. [`plotabm`](@ref) cannot be used explicitly (yet)
# since it expects our categorical dimension is actually a third spatial one.
# We start with some custom legends to easier understand the dynamics.

label_action = ["Battle", "Run", "Showdown", "Sneak", "Duel", "Captor", "Prisoner", "Chase"]
actions = [:rect, :utriangle, :circle, :pentagon, :diamond, :vline, :hline, :star4]
group_action = [
    MarkerElement(
        marker = marker,
        color = :black,
        strokecolor = :transparent,
        markersize = 15,
    ) for marker in actions
]
group_level = [
    PolyElement(color = color, strokecolor = :transparent) for color in cgrad(:tab10)[1:10]
]
nothing #hide

# And some complex internals that will be hidden away in the near future

e = size(model.space.s)[1:2] .+ 2
o = zero.(e) .- 2
clr(agent) = cgrad(:tab10)[level(agent)]
mkr(a) = a.shape
colors = Observable(to_color.([clr(model[id]) for id in Schedulers.by_id(model)]))
markers = Observable([mkr(model[id]) for id in Schedulers.by_id(model)])
pos = Observable([model[id].pos for id in Schedulers.by_id(model)])
stepper = InteractiveDynamics.ABMStepper(
    clr,
    mkr,
    15,
    nothing,
    Schedulers.by_id,
    pos,
    colors,
    Observable(15),
    markers,
    nothing,
    nothing
)
nothing #hide

# Finally, the plot:

f = Figure(resolution = (600, 700))
ax = f[1, 1] = Axis(f, title = "Battle Royale")
hidedecorations!(ax)
ax.xgridvisible = true
ax.ygridvisible = true
f[2, 1] = Legend(
    f,
    [group_action, group_level],
    [label_action, string.(1:10)],
    ["Action", "Level"],
    orientation = :horizontal,
    tellheight = true,
    tellwidth = false,
    nbanks = 5,
)

scatter!(ax, pos; color = colors, markersize = 15, marker = markers, strokewidth = 0.0)
xlims!(ax, o[1], e[1])
ylims!(ax, o[2], e[2])
record(f, "battle.mp4", 0:225; framerate = 10) do i
    Agents.step!(stepper, model, agent_step!, dummystep, 1)
end
nothing # hide
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../battle.mp4" type="video/mp4">
# </video>
# ```

# Some interesting behaviour emerges: sometimes you see a group of diamonds chasing one triangle.
# What ends up happening here is usually a close pair that wishes to fight gets caught
# out by the weaker one of the two running away from an even stronger opponent. Problem
# is that this stronger opponent is chasing the stronger of the pair, but since the
# weakest of the pair is still closer to the newcomer, there is a stalemate.
# This is usually resolved by hitting a boundary or other opponents.


text:
import ACE.SphericalHarmonics: SHBasis, index_y

@doc raw"""
`struct Ylm1pBasis <: OneParticleBasis`

One-particle basis of the form
```math
\phi_{lm}({\bm r}) = Y_l^m(\hat{\br r})
```
Fundamental building block of ACE basis sets of the form
```math
   R_{nl}^{\mu_i, \mu_j}(r_{ij}) Y_l^m(\hat{\bm r})
```
This type basically just translates the `SHBasis` into a valid one-particle
basis.
"""
mutable struct Ylm1pBasis{T, VSYM, LSYM, MSYM, TDX} <: OneParticleBasis{Complex{T}}
   SH::SHBasis{T}  # SH = Ylm
   B_pool::VectorPool{Complex{T}}
   dB_pool::VectorPool{TDX}
end

# # -------- temporary hack for 1.6, should not be needed from 1.7 onwards

# function acquire_B!(basis::Ylm1pBasis, args...)
#    VT = valtype(basis, args...)
#    return acquire!(basis.B_pool, length(basis), VT)
# end

# function release_B!(basis::Ylm1pBasis, B)
#    return release!(basis.B_pool, B)
# end

# ---------------------- Implementation of Ylm1pBasis

Ylm1pBasis(maxL::Integer, T = Float64; kwargs...) =
      Ylm1pBasis((SHBasis(maxL, T)); kwargs...)

Ylm1pBasis(SH::SHBasis{T}; varsym = :rr, lsym = :l, msym = :m)  where {T} =
      Ylm1pBasis{T, varsym, lsym, msym}(SH)

function Ylm1pBasis{T, varsym, lsym, msym}(SH::SHBasis{T}) where {T, varsym, lsym, msym}
   TDX = ACE.DState{NamedTuple{(varsym,), Tuple{SVector{3, Complex{T}}}}}
   return Ylm1pBasis{T, varsym, lsym, msym, TDX}(
            SH, VectorPool{Complex{T}}(), VectorPool{TDX}() )
end

Base.length(basis::Ylm1pBasis) = length(basis.SH)

_varsym(::Ylm1pBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = VSYM
_lsym(::Ylm1pBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = LSYM
_msym(::Ylm1pBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = MSYM
_l(b, basis::Ylm1pBasis) = getproperty(b, _lsym(basis))
_m(b, basis::Ylm1pBasis) = getproperty(b, _msym(basis))
_lm(b, basis::Ylm1pBasis) = _l(b, basis), _m(b, basis)

# -> TODO : figure out how to do this well!!!
# Base.rand(basis::Ylm1pBasis) =
#       AtomState(rand(basis.zlist.list), ACE.Random.rand_vec(basis.J))

function get_spec(basis::Ylm1pBasis, i::Integer)
   l, m = ACE.SphericalHarmonics.idx2lm(i)
   L, M = _lsym(basis), _msym(basis)
   return NamedTuple{(L, M)}((l, m))
end

get_spec(basis::Ylm1pBasis) = get_spec.(Ref(basis), 1:length(basis))

# function get_spec(basis::Ylm1pBasis{T, VS, L, M}) where {T, VS, L, M}
#    @assert length(basis) == (_maxL(basis) + 1)^2
#    spec = Vector{NamedTuple{(L, M), Tuple{Int, Int}}}(undef, length(basis))
#    for l = 0:_maxL(basis), m = -l:l
#       spec[index_y(l, m)] = NamedTuple{(L, M), Tuple{Int, Int}}( (l, m) )
#    end
#    return spec
# end

==(P1::Ylm1pBasis, P2::Ylm1pBasis) =
      ( (P1.SH == P2.SH) && (typeof(P1) == typeof(P2)) )

write_dict(basis::Ylm1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_Ylm1pBasis",
          "SH" => write_dict(basis.SH),
          "varsym" => _varsym(basis),
          "lsym" => _lsym(basis),
          "msym" => _msym(basis) )

read_dict(::Val{:ACE_Ylm1pBasis}, D::Dict) =
      Ylm1pBasis(read_dict(D["SH"]),
                 varsym = Symbol(D["varsym"]),
                 lsym = Symbol(D["lsym"]),
                 msym = Symbol(D["msym"]) )

# TODO: fix the type promotion...

valtype(basis::Ylm1pBasis{T}, args...) where T = Complex{T}

gradtype(basis::Ylm1pBasis, X::AbstractState) = dstate_type(valtype(basis, X), X)

symbols(basis::Ylm1pBasis) = [_lsym(basis), _msym(basis)]

argsyms(basis::Ylm1pBasis) = (_varsym(basis), )

_maxL(Ylm::Ylm1pBasis) = ACE.SphericalHarmonics.maxL(Ylm.SH)

function indexrange(basis::Ylm1pBasis)
   maxl = _maxL(basis)
   # note we create a stupid tensor product domain and then make up for it
   # by using an index filter during the basis generation process
   return NamedTuple{(_lsym(basis), _msym(basis))}( (0:maxl, -maxl:maxl) )
end

function isadmissible(b, basis::Ylm1pBasis)
   l, m = _lm(b, basis)
   maxL = _maxL(basis)
   return (0 <= l <= maxL) && (-l <= m <= l)
end

# ---------------------------  Evaluation code
#

_rr(X, basis::Ylm1pBasis) = getproperty(X, _varsym(basis))

evaluate!(B, basis::Ylm1pBasis, X::AbstractState) =
      evaluate!(B, basis.SH, _rr(X, basis))

function evaluate_d!(dB, basis::Ylm1pBasis, X::AbstractState)
   B = acquire_B!(basis, X)
   evaluate_ed!(B, dB, basis, X)
   release_B!(basis, B)
   return dB
end

function evaluate_ed!(B, dB, basis::Ylm1pBasis, X::AbstractState)
   TDX = eltype(dB)  # need not be the same as gradtype!!!
   RSYM = _varsym(basis)
   rr = _rr(X, basis)
   dY = acquire_dB!(basis.SH, rr)
   # spherical harmonics does only ed since values are essentially free
   evaluate_ed!(B, dY, basis.SH, rr)
   for n = 1:length(basis)
      dB[n] = TDX( NamedTuple{(RSYM,)}((dY[n],)) )
   end
   release_dB!(basis.SH, dY)
   return B, dB
end

degree(b, Ylm::Ylm1pBasis) = _l(b, Ylm)

degree(b, Ylm::Ylm1pBasis, weight::Dict) = haskey(weight,_lsym(Ylm)) ? weight[_lsym(Ylm)] * degree(b, Ylm) : degree(b, Ylm)

get_index(Ylm::Ylm1pBasis, b) = index_y(_l(b, Ylm), _m(b, Ylm))

# -------------- AD

import ChainRules: rrule, NoTangent, @not_implemented

function _rrule_evaluate(basis::Ylm1pBasis, X::AbstractState,
                         w::AbstractVector{<: Number})
   dY = evaluate_d(basis, X)
   a = sum( real(w) * real.(d.rr) + imag(w) * imag.(d.rr)
            for (w, d) in zip(w, dY)  )
   TDX = dstate_type(a, X)
   return TDX( NamedTuple{(_varsym(basis),)}( (a,) ) )
end

rrule(::typeof(evaluate), basis::Ylm1pBasis, X::AbstractState) =
      evaluate(basis, X),
      w -> (NoTangent(), NoTangent(), _rrule_evaluate(basis, X, w))

rrule(::typeof(evaluate_d), basis::Ylm1pBasis, X::AbstractState) =
      evaluate_d(basis, X),
      w -> (NoTangent(), NoTangent(),
            @not_implemented("""Ylm config gradients is currently not implemented
                                composition with vectorial features is therefore
                                not yet supported."""))


text:
function metric_head(h,T, q_o,constants, soil_parameters, soil_numerical_parameters, dt)

# initializing
h_old = h;
theta_l, theta_v, soil_hydraulic = hydraulic_variables(h, T, constants, soil_parameters);
rho_w = 1000 .-7.3e-3.*(T.-(273+4)).^2 .+ 3.79e-5.*(T.-(273+4)).^3;  #water density [kg m^-3]
q_o = q_o./rho_w[end];

O=zeros(2)
O[1] = 1;

while O[1]>soil_numerical_parameters.eps_h
    temp_h_2 = h;
    temp_h = (h+h_old)./2;

    for i = 2:soil_numerical_parameters.Ns-1
        h[i] = h_old[i] + dt/(soil_numerical_parameters.dz_s^2*soil_hydraulic.C[i])*
            (soil_hydraulic.K_h[i]*(temp_h[i+1]-2*temp_h[i]+temp_h[i-1]) +
            (soil_hydraulic.K_h[i+1]-soil_hydraulic.K_h[i])*(temp_h[i+1]-temp_h[i]) +
            soil_hydraulic.K_T[i]*(T[i+1]-2*T[i]+T[i-1]) +
            (soil_hydraulic.K_T[i+1]-soil_hydraulic.K_T[i])*(T[i+1]-T[i]) +
            (soil_hydraulic.K[i+1]-soil_hydraulic.K[i])*soil_numerical_parameters.dz_s);
    end
    h[1] = h_old[1]+dt/(soil_numerical_parameters.dz_s*soil_hydraulic.C[1])*(soil_hydraulic.K[2]-soil_hydraulic.K[1]);
    h[end] = h_old[end]-dt/(soil_hydraulic.C[end]*soil_numerical_parameters.dz_s^2)*(soil_hydraulic.K_h[end]*(temp_h[end]-temp_h[end-1])+
        soil_hydraulic.K_T[end]*(T[end]-T[end-1])+(soil_hydraulic.K[end]+q_o)*soil_numerical_parameters.dz_s);

    h[h.<soil_parameters.h_max] .= soil_parameters.h_max;

    theta_l, theta_v, soil_hydraulic, rho_vs = hydraulic_variables(h, T, constants, soil_parameters);

    O[2] = mean(abs.(temp_h_2-h)./abs.(temp_h));

    if O[2]<O[1]
        O[1] = O[2];
    else
        dt = dt/2;
        h = h_old;
        theta_l, theta_v, soil_hydraulic = hydraulic_variables(h, T, constants, soil_parameters);
        O[1] = 1;
    end
    #println(dt)
end

q_l[1] = -soil_hydraulic.K[1];
q_v[1] = 0;

for i = 2:soil_numerical_parameters.Ns
    q_l[i] = -soil_hydraulic.K_lT[i]*(T[i]-T[i-1])/soil_numerical_parameters.dz_s-
        soil_hydraulic.K[i]*((h[i]-h[i-1])/soil_numerical_parameters.dz_s+1);
    q_v[i] = -soil_hydraulic.K_vT[i]*(T[i]-T[i-1])/soil_numerical_parameters.dz_s-
        soil_hydraulic.K_vh[i]*((h[i]-h[i-1])/soil_numerical_parameters.dz_s);
end

return h, theta_l, theta_v, q_l, q_v, rho_vs, dt

end


text:
module types

export Node

#=
    Node
=#
"""
"""
mutable struct Node
    id::Int
    label::String
    parent::Int
    children::Vector{Int}
    siblings::Vector{Int}
    support::Int

    function Node(; id = 0,
                    label = "",
                    parent = 0,
                    children = Int[],
                    siblings = Int[],
                    support = 0)

        new(id, label, parent, children, siblings, support)
    end

    Node(id, label, parent, children, siblings, support) = new(id, label, parent, children, siblings, support)

end

#=
    Tree
=#
"""
"""
mutable struct Tree
    nodes::Vector{Node}
    ids::Vector{Int}
    labels::Vector{String}
    size::Int
end

#=
    Pretty print functions
=#
"""
"""
function Base.show(io::IO, tree::Tree)
    @info "FP Growth Tree"
    @info "Nodes..."
    foreach(x -> println("--", x), tree.nodes)
end

end # module

max_stars_repo_name: DopedCow/FPGrowth

max_issues_repo_name: DopedCow/FPGrowth

max_forks_repo_name: DopedCow/FPGrowth


text:
# Generate a random graph with n vertices where each edge is included with probability p.
function erdos_renyi_graph(g::GT, n::Integer, p::Real; has_self_loops=false) where {GT<:AbstractGraph}
    for i=1:n
        start_ind = is_directed(g) ? 1 : i
        for j=start_ind:n
            if(rand() <= p && (i != j || has_self_loops))
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

# Convenience function with a default graph type.
function erdos_renyi_graph(n::Integer, p::Real; is_directed=true, has_self_loops=false)
    g = simple_inclist(n, is_directed=is_directed)
    erdos_renyi_graph(g, n, p, has_self_loops=has_self_loops)
end

# Generate a 'small world' random graph based on the Watts-Strogatz model.
# Written with much reference to the implementation from GraphStream <http://graphstream-project.org>.
# The resulting graph has n vertices,
#   Each vertex has a base degree of k  (n > k, k >= 2, k must be even.)
#   There is a beta chance of each edge being 'rewired'
function watts_strogatz_graph(g::GT, n::Integer, k::Integer, beta::Real) where {GT<:AbstractGraph}
    # When attributes are restored, we can enable placing the vertices around a circle.
    # This is nice for explicit visualization of the model.
    #space = linspace(0,2*pi,n+1)
    #for v in vertices(g)
    #    attrs = attributes(v)
    #    x = round((cos(space[i]) + 1)*100,2)
    #    y = round((sin(space[i]) + 1)*100,2)
    #    attrs["pos"] = "$(x),$(y)"
    #end

    # Link each node to the k/2 nodes next to it in each direction.
    # For each vertex in the graph:
    for i in 1:n
        # For the next k/2 higher indexed nodes:
        for j in 1:div(k,2)
            # Either make an edge to the neighbor or a long edge across the graph.
            if rand() > beta
                add_edge!(g, i, ((i+j - 1) % n) + 1 )
            else
                # Add a random link across the graph and delete the link to the neighbor.
                while true
                    target = rand(1:(n-1))
                    if (target >= i)
                        target += 1
                    end
                    if !(target in out_neighbors(i,g))
                        add_edge!(g,i, target)
                        break
                    end
                end
            end
        end
    end
    g
end

# Convenience function with a default graph type.
function watts_strogatz_graph(n::Integer, k::Integer, beta::Real)
    g = simple_inclist(n, is_directed=false)
    watts_strogatz_graph(g,n,k,beta)
end


text:
using JuMP, CPLEX

#Preparing optimization model
m = Model(CPLEX.Optimizer)
model = Model(CPLEX.Optimizer)
set_optimizer_attribute(m, "CPX_PARAM_EPINT", 1e-8)

c = [ 1; 2; 5]
A = [-1  1  3;
      1  3 -7]
b = [-5; 10]

println("Size of c: ", size(c))
println("Size of A: ", size(A))
println("Size of b: ", size(b))

@variable(m, x[1:3] >= 0)
@objective(m, Max, sum( c[i]*x[i] for i=1:3) )

@constraint(m, constraint[j=1:2], sum( A[j,i]*x[i] for i=1:3 ) <= b[j] )

#@constraint(m, constraint1, sum( A[1,i]*x[i] for i=1:3) <= b[1] )
#@constraint(m, constraint2, sum( A[2,i]*x[i] for i=1:3) <= b[2] )

@constraint(m, bound, x[1] <= 10)

optimize!(m)

print(m)

println("Optimal Solutions:")
for i=1:3
  println("x[$i] = ", JuMP.value(x[i]))
end

println("Dual Variables:")
for j=1:2
  println("dual[$j] = ", JuMP.dual(constraint[j]))
end

println("Size of x: ", size(x))

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using Bokeh

x = linspace(-3pi, 3pi)
y1 = cos.(x)
plot(x,y1)

# hold works very similarly to hold in other plotting libraries
hold(true)

y2 = sin.(x) ./ x
plot(x, y2, "r")
showplot("holding.html")
# under the bonnet it's just causing the plot function not to overwrite
# CURPLOT but rather to extend it

# we can do something similar bey specifying exactly which plot we wish to
# extend
hold(false)
# becaues hold is false p1 and p2 are seperate plots
p1 = plot(x, y1, plotfile="plotobject1.html")
p2 = plot(x, y2, plotfile="plotobject2.html")

y3 = y1 - y2

# if we switched hold on this would extend p2, instead we specify that it should
# extend p1
plot(x, y3, "r--", extend = p1)
showplot(p1)
showplot(p2)


text:
"""
    calibrate(c, xy)

Calibrate a coordinate, transforming it from pixel to real-world
coordinate. `xy` must be a `StaticArrays.SVector{2, Float64}`.
"""
calibrate(c, xy::V2) = c.tform(xy)

"""
    calibrate(c, img)

Calibrate an image, transforming its coordinates from pixel to real-world.
"""
function calibrate(c, img)
    indices = ImageTransformations.autorange(img, c.tform)
    warp(img, c.itform, indices)
end


text:
using Documenter, Jack

makedocs(;
    modules=[Jack],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/TsuMakoto/Jack.jl/blob/{commit}{path}#L{line}",
    sitename="Jack.jl",
    authors="TsuMakoto",
    assets=String[],
)

deploydocs(;
    repo="github.com/TsuMakoto/Jack.jl",
)


text:
using RustModels
using Base.Test

import Markov.rsm

include("toymodels.jl")
include("classicalrustmodels.jl")
include("indexmanipulation.jl")


text:
"""
v1 but main loop until w-1 and additional end test
"""
function bloom_v2 end

function bloom_v2(t::SearchSequence)
    n = sizeof(t)
    bloom_mask = UInt64(0)
    skip = n - 1
    tlast = _nthbyte(t,n)
    for j in 1:n
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tlast && j < n
            skip = n - j -1
        end
    end
    return t,bloom_mask,skip,tlast
end

function bloom_v2(s::SearchSequence, t::SearchSequence, i::Integer)
    bloom_v2(s,bloom_v2(t),i)
end

function bloom_v2(s::SearchSequence,p::Tuple, i::Integer,sv::MaybeVector=nothing)
    (t,bloom_mask,skip,tlast) = p # replaces original code, commented out
    DOSTATS = !(sv isa Nothing)
    n = sizeof(t)
    m = sizeof(s)
    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return something(findnext(isequal(_nthbyte(t,1)), s, i), 0)
    end
    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

    if DOSTATS loops = 0 end
    if DOSTATS bloomtests = 0 end
    if DOSTATS bloomskips = 0 end
    i -= 1
    while i < w
        if DOSTATS loops += 1 end
        if _nthbyte(s,i+n) == tlast
            # check candidate
            j = 1
            while j < n
                if _nthbyte(s,i+j) != _nthbyte(t,j)
                    break
                end
                j += 1
            end
            # match found?
            if j == n
                if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
                return i+1
            end
            # no match: skip and test bloom
            i += skip
            if i>=w
                break # allows for unconditional bloom test below
            end
        end
        i += 1
        if DOSTATS bloomtests += 1 end
        # i<=w is guaranteed here, comparison eliminated
        if bloom_mask & _search_bloom_mask(_nthbyte(s,i+n)) == 0
            if DOSTATS bloomskips += 1 end
            i += n
        end
    end
    if i==w
        # test end match
        j = 1
        while j <= n
            if _nthbyte(s,i+j) != _nthbyte(t,j)
                break # not found
            end
            j += 1
        end
        if j > n
            if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
            return i+1
        end
    end
    if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
    0
end


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

JLLWrappers.@generate_main_file_header("LibUnwind")
JLLWrappers.@generate_main_file("LibUnwind", UUID("745a5e78-f969-53e9-954f-d19f2f74f4e3"))
end  # module LibUnwind_jll


text:
using TimesTables
using Test

@testset "Multiplication" begin
  m = Multiplication(3,4)
  @test compute(m) == 12
  @test check(m,12)
  @test !check(m,13)
end

@testset "Division" begin
  @test_throws ErrorException d = TimesTables.Division(63,6)

  d = Division(63,7)
  @test compute(d) == 9
  @test check(d,9)
  @test !check(d,7)
end

@testset "Addition" begin
  m = Addition(3,4)
  @test compute(m) == 7
  @test check(m,7)
  @test !check(m,8)
end

@testset "Subtraction" begin
  m = Subtraction(4,3)
  @test compute(m) == 1
  @test check(m,1)
  @test !check(m,2)
end

@testset "Generation" begin
  @test typeof(generate(Multiplication, DefaultOpGenerator())) == Multiplication
  @test typeof(generate(Division, DefaultOpGenerator())) == Division
  @test typeof(generate(Addition, DefaultOpGenerator())) == Addition
  @test typeof(generate(Subtraction, DefaultOpGenerator())) == Subtraction
  @show g = generate()
  @test check(g, compute(g))
  @test !check(g, 2*compute(g))
  println(g)
end


text:
# ---
# title: Composition utilities
# cover: utilities.gif
# description: a set of helper operations that may be useful when compositing more complex augmentation workflow
# ---

# Aside from "true" operations that specify some kind of transformation, there are also a couple of
# special utility operations used for functionality such as stochastic branching.

using Augmentor
using Random
Random.seed!(1337)

img_in = testpattern(RGB, ratio=0.5)
img_out = augment(img_in, Either(0.5=>NoOp(), 0.25=>FlipX(), 0.25=>FlipY()))
#md nothing #hide

# ![](utilities.gif)

# ## References

#md # ```@docs
#md # NoOp
#md # Either
#md # CacheImage
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
pl = Either(ShearX(-5:5), ShearY(-5:5)) |> Rotate(-10:10) |> Either(NoOp(), FlipX(), FlipY()) |> CropNative(axes(img_in)) #src
cover = make_gif(img_in, pl, 10) #src
ImageMagick.save("utilities.gif", cover; fps=2) #src


text:
"""
    calculatemoments(dist)

Computes the lattice Boltzmann moments (h,v) from a distribution function `dist`.

Based on the input either a `D1Q3` or a `D2Q9` distribution this function performs a summation over the particle speeds.
In one spatial dimension the velocity set [-1, 0, 1] allows for three speeds in two for nine speeds.
See also: [`dist2array`](@ref)

# Math
The moments are calculated according to

`` \\rho = \\sum_{i=0}^N f_i ``

`` \\rho\\mathbf{u} = \\sum_{i=0}^N \\mathbf{c}_i f_i ``

`` e = \\sum_{i=0}^N \\mathbf{c}_i \\mathbf{c}_i f_i ``

Where N is either 3 or 9. In one spatial dimension the velocity is of course not a vector of ``(x,y)``.

# Example
```jldoctest
julia> using JuSwalbe

julia> a = fill(0.1, (10, 10))
10×10 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1

julia> b = fill(0.2, (10,10))
10×10 Array{Float64,2}:
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2

julia> dist = JuSwalbe.DistributionD2Q9(f0=a,f1=b,f2=a,f3=a,f4=a,f5=a,f6=a,f7=a,f8=a)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f1: Array{Float64}((10, 10)) [0.2 0.2 … 0.2 0.2; 0.2 0.2 … 0.2 0.2; … ; 0.2 0.2 … 0.2 0.2; 0.2 0.2 … 0.2 0.2]
  f2: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f3: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f4: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f5: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f6: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f7: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f8: Array{Float64}((10, 10)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]

julia> mom = calculatemoments(dist)
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((10, 10)) [0.9999999999999999 0.9999999999999999 … 0.9999999999999999 0.9999999999999999; 0.9999999999999999 0.9999999999999999 … 0.9999999999999999 0.9999999999999999; … ; 0.9999999999999999 0.9999999999999999 … 0.9999999999999999 0.9999999999999999; 0.9999999999999999 0.9999999999999999 … 0.9999999999999999 0.9999999999999999]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((10, 10)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((10, 10)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> mom.height
10×10 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0

julia> mom.velocity.x
10×10 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1

julia> mom.velocity.y
10×10 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> dist1d = JuSwalbe.DistributionD1Q3(f0=a[1:10],f1=b[1:10],f2=a[1:10]) # Vectorize a and b
JuSwalbe.DistributionD1Q3{Array{Float64,1}}
  f0: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  f1: Array{Float64}((10,)) [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2]
  f2: Array{Float64}((10,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]

julia> mom1d = calculatemoments(dist1d)
JuSwalbe.Macroquant{Array{Float64,1},Array{Float64,1}}
  height: Array{Float64}((10,)) [0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4]
  velocity: Array{Float64}((10,)) [0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25]
  pressure: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  energy: Array{Float64}((10,)) [0.30000000000000004, 0.30000000000000004, 0.30000000000000004, 0.30000000000000004, 0.30000000000000004, 0.30000000000000004, 0.30000000000000004, 0.30000000000000004, 0.30000000000000004, 0.30000000000000004]

julia> mom1d.height
10-element Array{Float64,1}:
 0.4
 0.4
 0.4
 0.4
 0.4
 0.4
 0.4
 0.4
 0.4
 0.4

julia> mom1d.velocity
10-element Array{Float64,1}:
 0.25
 0.25
 0.25
 0.25
 0.25
 0.25
 0.25
 0.25
 0.25
 0.25

```

# References
## Code
See also: [`dist2array`](@ref)

## Two spatial dimensions
Here I go with Dellar
- [Nonhydrodynamic modes and *a priori* construction of shallow water lattice Boltzmann equations](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.65.036309)

## One spatial dimension
For the one dimensional simulation the paper by Chopard is quite nice
- [Study of the 1D lattice Boltzmann shallow water equation and its coupling to build a canal network.](https://www.sciencedirect.com/science/article/pii/S0021999110003372)

"""
function calculatemoments(dist::JuSwalbe.DistributionD1Q3{Vector{T}}) where {T<:Number}
    # Get a size for the output arrays
    len = length(dist.f0)
    # Allocate the output
    height = zeros(T, len)
    velocity = zeros(T, len)
    pressure = zeros(T, len)
    energy = zeros(T, len)
    # Lattice velocities
    c = [0 1 -1]
    csquare = [0 1 1]
    # Store all distributions in an array
    distarray = hcat(dist.f0, dist.f1, dist.f2)
    height = sum(distarray, dims=2)[:, 1]
    velocity = sum(distarray .* c, dims=2)[:, 1]
    energy = sum(distarray .* csquare, dims=2)[:, 1]

    velocity ./= height

    result = JuSwalbe.Macroquant(height, velocity, pressure, energy)
    return result
end

function calculatemoments(dist::JuSwalbe.DistributionD1Q3{Vector{T}}, force::JuSwalbe.Forces{Vector{T}}) where {T<:Number}
    # Get a size for the output arrays
    len = length(dist.f0)
    # Allocate the output
    height = zeros(T, len)
    velocity = zeros(T, len)
    pressure = zeros(T, len)
    energy = zeros(T, len)
    # Lattice velocities
    c = [0 1 -1]
    csquare = [0 1 1]
    # Store all distributions in an array
    distarray = hcat(dist.f0, dist.f1, dist.f2)
    # Store forces in array
    forcearray = hcat(force.slip, force.thermal, force.h∇p)
    # Moment summation is acutally performed here
    height = sum(distarray, dims=2)[:, 1]
    velocity = sum(distarray .* c, dims=2)[:, 1] + T(0.5) * sum(forcearray, dims=2)[:, 1]
    energy = sum(distarray .* csquare, dims=2)[:, 1]
    # Last step is to divide the velocity by h
    velocity ./= height
    # Return the output
    result = JuSwalbe.Macroquant(height, velocity, pressure, energy)
    return result
end
# Below is the two dimensional thing
function calculatemoments(dist::JuSwalbe.DistributionD2Q9{Matrix{T}}) where {T<:Number}
    # Get a size for the output arrays
    width, thick = size(dist.f0)
    # Allocate the output
    height = zeros(T, (width,thick))
    velocityx = zeros(T, (width,thick))
    velocityy = zeros(T, (width,thick))
    pressure = zeros(T, (width,thick))
    energy = zeros(T, (width,thick))
    # Lattice velocities
    lattice_vel = [0.0 0.0; 1.0 0.0; 0.0 1.0; -1.0 0.0; 0.0 -1.0; 1.0 1.0; -1.0 1.0; -1.0 -1.0; 1.0 -1.0]
    indices = [Symbol("f$(i-1)") for i in 1:9]
    # Store all distributions in an array
    distarray = dist2array(dist)
    # The moments are just a summation of distribution functions
    height = sum(distarray, dims=1)[1, :, :]
    velocityx = sum(distarray .* lattice_vel[:, 1], dims=1)[1, :, :]
    velocityy = sum(distarray .* lattice_vel[:, 2], dims=1)[1, :, :]

    velocityx ./= height
    velocityy ./= height

    velocity = JuSwalbe.Twovector(velocityx, velocityy)

    result = JuSwalbe.Macroquant(height, velocity, pressure, energy)
    return result
end

function calculatemoments(dist::JuSwalbe.DistributionD2Q9{Matrix{T}}, force::JuSwalbe.Forces{Twovector{Matrix{T}}}) where {T<:Number}
    # Get a size for the output arrays
    width, thick = size(dist.f0)
    # Allocate the output
    height = zeros(T, (width,thick))
    velocityx = zeros(T, (width,thick))
    velocityy = zeros(T, (width,thick))
    pressure = zeros(T, (width,thick))
    energy = zeros(T, (width,thick))
    # Lattice velocities
    lattice_vel = [0.0 0.0; 1.0 0.0; 0.0 1.0; -1.0 0.0; 0.0 -1.0; 1.0 1.0; -1.0 1.0; -1.0 -1.0; 1.0 -1.0]
    indices = [Symbol("f$(i-1)") for i in 1:9]
    # Store all distributions in an array
    distarray = dist2array(dist)
    forces = dist2array(force)
    # The moments are just a summation of distribution functions
    height .= sum(distarray, dims=1)[1, :, :]
    velocityx .= sum(distarray .* lattice_vel[:, 1], dims=1)[1, :, :] .+ T(0.5) * sum(forces.x, dims=3)[:, :, 1]
    velocityy .= sum(distarray .* lattice_vel[:, 2], dims=1)[1, :, :] .+ T(0.5) * sum(forces.y, dims=3)[:, :, 1]

    velocityx ./= height
    velocityy ./= height

    velocity = JuSwalbe.Twovector(velocityx, velocityy)

    result = JuSwalbe.Macroquant(height, velocity, pressure, energy)
    return result
end

"""
    dist2array(dist)

Transforms the `D2Q9` struct into a three dimensional array.

Axis are chosen such that it can be easily multiplied with the nine speed lattice speeds.
The function also excepts a `Forces` struct as input and returns a three dimensional array for both `x` and `y`.

# Example
```jldoctest
julia> using JuSwalbe

julia> a = fill(0.1, (5,5))
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> b = fill(0.2, (5,5))
5×5 Array{Float64,2}:
 0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2
 0.2  0.2  0.2  0.2  0.2

julia> dist = JuSwalbe.DistributionD2Q9(f0=a, f1=b, f2=a, f3=a, f4=a, f5=a, f6=a, f7=a, f8=a)
JuSwalbe.DistributionD2Q9{Array{Float64,2}}
  f0: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f1: Array{Float64}((5, 5)) [0.2 0.2 … 0.2 0.2; 0.2 0.2 … 0.2 0.2; … ; 0.2 0.2 … 0.2 0.2; 0.2 0.2 … 0.2 0.2]
  f2: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f3: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f4: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f5: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f6: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f7: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]
  f8: Array{Float64}((5, 5)) [0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1; … ; 0.1 0.1 … 0.1 0.1; 0.1 0.1 … 0.1 0.1]

julia> arr = dist2array(dist)
9×5×5 Array{Float64,3}:
[:, :, 1] =
 0.1  0.1  0.1  0.1  0.1
 0.2  0.2  0.2  0.2  0.2
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

[:, :, 2] =
 0.1  0.1  0.1  0.1  0.1
 0.2  0.2  0.2  0.2  0.2
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

[:, :, 3] =
 0.1  0.1  0.1  0.1  0.1
 0.2  0.2  0.2  0.2  0.2
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

[:, :, 4] =
 0.1  0.1  0.1  0.1  0.1
 0.2  0.2  0.2  0.2  0.2
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

[:, :, 5] =
 0.1  0.1  0.1  0.1  0.1
 0.2  0.2  0.2  0.2  0.2
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

```

"""
function dist2array(dist::JuSwalbe.DistributionD2Q9{Matrix{T}}) where {T<:Number}
    # Get the size to allocate memory
    width, thick = size(dist.f0)
    # Symbols to access the fields: :f0, :f1, :f2 ...
    indices = [Symbol("f$(i-1)") for i in 1:9]
    # Allocate the output array
    distarray = Array{T, 3}(undef, (9,width,thick))
    for (i, j) in enumerate(indices)
        distarray[i, :, :] = getfield(dist, j)
    end
    return distarray
end

function dist2array(force::JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{T}}}) where {T<:Number}
    # Get the size to allocate memory
    width, thick = size(force.slip.x)
    # Symbols to access the fields:
    indices = [:slip :h∇p :bathymetry :thermal]
    # Allocate the output array
    forcearrayx = Array{T, 3}(undef, (width,thick,4))
    forcearrayy = Array{T, 3}(undef, (width,thick,4))
    for (i, j) in enumerate(indices)
        forcearrayx[:, :, i] = getfield(getfield(force, j),:x)
        forcearrayy[:, :, i] = getfield(getfield(force, j),:y)
    end
    forces = JuSwalbe.Twovector(x=forcearrayx, y=forcearrayy)
    return forces
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
# -*- coding: utf-8 -*-
using SemanticModels.ModelTools

include("../src/modeltools/typegraph.jl")

edgelist = @typegraph begin
    import Base: count

    """    AgentModel

    the root type for agent based models.

    See also: StateModel
    """
    abstract type AgentModel end

    """     StateModel

    holds the components of an agent based simulation using finite state machines.

    - states: a collection of distinct states an agent can occupy
    - agents: a collection of states `aᵢ = agents[i] ∈ states` indicating that agent `i` is in state `aᵢ`
    - transitions: the functions `f: states -> states`
    """
    mutable struct StateModel <: AgentModel
        states
        agents
        transitions
        loads
    end

    # +
    # Determine number of agents in a given state
    function count(sm::StateModel, state)
        return length(filter(a->(a==state), sm.agents))
    end

    # counts the number of agents
    function count(sm::StateModel)
        return length(sm.agents)
    end
    # -

    # stateload computes the fraction of agents in each state
    # it is used by tick! to update the statemodel for computing the probability of infection.
    function stateload(sm::StateModel, state::Symbol)
        return (count(sm, state)+1)/(count(sm)+1)
    end

    #     tick!(sm::StateModel)
    #
    # performs the operations that need to happen once per time step. Things like caching values that need to be computed at the beginning of the timestep like the distribution of the agents' states.

    function tick!(sm::StateModel)
        sm.loads = map(s->stateload(sm, s), sm.states)
    end

    #     step!(sm::StateModel, n=1)
    #
    # advance the simulation by `n` ticks of time.
    # This is an in-place operation that modifies the current state of the simulation.
    #

    function step!(sm::StateModel, n)
        for s in 1:n
          tick!(sm)
          for (i, a) in enumerate(sm.agents)
              sm.agents[i] = sm.transitions[a](sm, i,a)
          end
            describe(sm)
        end
        return sm
    end

    """    describe(sm::StateModel)

    summarize the state of the simulation for presentation or analysis.
    """
    function describe(sm::StateModel)
        counts = zeros(Int, size(sm.states))
        d = Dict{eltype(sm.states), Int}()
        for (i,s) in enumerate(sm.states)
            d[s] = i
        end
        for a in sm.agents
            i = d[a]
            counts[i] += 1
        end
        return collect(map(x->Pair(x...), zip(sm.states, counts)))
    end

    # ## Run the model
    #
    # This script defines a basic agent based model of disease spread called SIRS. Each agent is in one of 3 states
    #
    # 1. $S$ Susceptible
    # 2. $I$ Infected
    # 3. $R$ Recovered
    #
    # <img src="https://docs.google.com/drawings/d/e/2PACX-1vSeA7mAQ-795lLVxCWXzbkFQaFOHMpwtB121psFV_2cSUyXPyKMtvDjssia82JvQRXS08p6FAMr1hj1/pub?w=1031&h=309">
    #
    # The agents go from `S->I`, `I-R`, and `R->S` based on random numbers. The probability of S-> is dependent on the fraction of agents in state :I. The probability of recovering is a constant ρ, and the disease confers some temporary immunity with probability μ.
    #
    # This script has an entrypoint to call it so that you can include this file and run as many simulations as you want. The intended use case is to repeatedly call `main` and accumulate the return values into an array for later analysis.

    function main(nsteps)
        n = 20
        a = fill(:S, n)
        ρ = 0.5 + randn(Float64)/4 # chance of recovery
        μ = 0.5 # chance of immunity
        T = Dict(
            :S=>(x...)->rand(Float64) < stateload(x[1], :I) ? :I : :S,
            :I=>(x...)->rand(Float64) < ρ ? :I : :R,
            :R=>(x...)->rand(Float64) < μ ? :R : :S,
        )

        sam = StateModel([:S, :I, :R], a, T, zeros(Float64,3))
        newsam = step!(deepcopy(sam), nsteps)
        @show newsam.agents
        counts = describe(newsam)
        return newsam, counts
    end

    # +
    # # An example of how to run this thing.
    main(10)
    # -

    end

E = unique((f.func, f.args, f.ret) for f in edgelist)

using Pkg
try
    using MetaGraphs
catch

    Pkg.add("MetaGraphs")
finally
    using MetaGraphs
end

g = MetaDiGraph()

# +
# TODO construct graph from E

# +
# TODO draw graph with graphviz

# +
# TOOD insert graph drawing

# +
# TODO use a different script and repeat, extraction, drawing, and rendering

# +
# TODO discuss difference in type graphs


text:
function finish()
  println("finishing...")
  output(t, 0)
end


text:
using Wikidata
using Test
using HTTP

@testset "Wikidata.jl" begin
    entity = Wikidata.WikidataEntity("Q42")
    @test  Wikidata.label(entity)=="Douglas Adams"
    @test_throws HTTP.ExceptionRequest.StatusError Wikidata.WikidataEntity("Douglas Adams")
    @test Wikidata.hasproperty(entity, "P31")== true
    @test Wikidata.hasproperty(entity, "P38")== false

    @test_throws ArgumentError Wikidata.getproperty(entity, "P38")
    propList = Wikidata.getproperty(entity, "P31")
    @test size(propList, 1) == 1
    @test Wikidata.label(propList[1]) == "human"

    cambridge = Wikidata.getproperty(entity, "P19")[1]
    @test Wikidata.label(cambridge) == "Cambridge"
    latlon = Wikidata.getproperty(cambridge, "P625")[1]
    @test latlon[1] == 52.208055555556
    @test latlon[2] == 0.1225
end

max_stars_repo_name: KimBue/Wikidata


text:
module LibFoo
using Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LibFoo not installed properly, run Pkg.build(\"LibFoo\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

# Export our two super-useful functions
export call_fooifier, call_libfoo

# Function to call the `fooifier` binary with the given arguments
function call_fooifier(a, b)
    global fooifier
    return parse(Float64, readchomp(`$fooifier $a $b`))
end

# Function to call into the `libfoo` shared library with the given arguments
function call_libfoo(a, b)
    global libfoo

    hdl = Libdl.dlopen_e(libfoo)
    @assert hdl != C_NULL "Could not open $libfoo"
    foo = Libdl.dlsym_e(hdl, :foo)
    @assert foo != C_NULL "Could not find foo() within $libfoo"
    return ccall(foo, Cdouble, (Cdouble, Cdouble), a, b)
end

end #module LibFoo


text:
module Capillary
import GSL: sf_bessel_zero_Jnu
import SpecialFunctions: besselj
import StaticArrays: SVector
import Cubature: hquadrature
using Reexport
@reexport using Luna.Modes
import Luna: Maths, Grid
import Luna.PhysData: c, ε_0, μ_0, ref_index_fun, roomtemp, densityspline, sellmeier_gas
import Luna.Modes: AbstractMode, dimlimits, neff, field, Aeff, N
import Luna.LinearOps: make_linop, conj_clamp, neff_grid, neff_β_grid
import Luna.PhysData: wlfreq
import Luna.Utils: subscript
import Base: show

export MarcatiliMode, dimlimits, neff, field, N, Aeff

"""
    MarcatiliMode

Type representing a mode of a hollow capillary as presented in:

Marcatili, E. & Schmeltzer, R.
"Hollow metallic and dielectric waveguides for long distance optical transmission and lasers
(Long distance optical transmission in hollow dielectric and metal circular waveguides,
examining normal mode propagation)."
Bell System Technical Journal 43, 1783–1809 (1964).
"""
struct MarcatiliMode{Ta, Tcore, Tclad, LT} <: AbstractMode
    a::Ta # core radius callable as function of z only, or fixed core radius if a Number
    n::Int # azimuthal mode index
    m::Int # radial mode index
    kind::Symbol # kind of mode (transverse magnetic/electric or hybrid)
    unm::Float64 # mth zero of the nth Bessel function of the first kind
    ϕ::Float64 # overall rotation angle of the mode
    coren::Tcore # callable, returns (possibly complex) core ref index as function of ω
    cladn::Tclad # callable, returns (possibly complex) cladding ref index as function of ω
    model::Symbol # if :full, includes complete influence of complex cladding ref index
    loss::LT # Val{true}() or Val{false}() - whether to include the loss
    aeff_intg::Float64 # Pre-calculated integral fraction for effective area
end

function show(io::IO, m::MarcatiliMode)
    a = radius_string(m)
    loss = "loss=" * (m.loss == Val(true) ? "true" : "false")
    model = "model="*string(m.model)
    angle = "ϕ=$(m.ϕ/π)π"
    out = "MarcatiliMode{"*join([mode_string(m), a, loss, model, angle], ", ")*"}"
    print(io, out)
end

mode_string(m::MarcatiliMode) = string(m.kind)*subscript(m.n)*subscript(m.m)
radius_string(m::MarcatiliMode{<:Number, Tco, Tcl, LT}) where {Tco, Tcl, LT} = "a=$(m.a)"
radius_string(m::MarcatiliMode) = "a(z=0)=$(radius(m, 0))"

"""
    MarcatiliMode(a, n, m, kind, ϕ, coren, cladn; model=:full, loss=true)

Create a MarcatiliMode.

# Arguments
- `a` : Either a `Number` for constant core radius, or a function `a(z)` for variable radius.
- `n::Int` : Azimuthal mode index (number of nodes in the field along azimuthal angle).
- `m::Int` : Radial mode index (number of nodes in the field along radial coordinate).
- `kind::Symbol` : `:TE` for transverse electric, `:TM` for transverse magnetic,
                   `:HE` for hybrid mode.
- `ϕ::Float` : Azimuthal offset angle (for linearly polarised modes, this is the angle
                between the mode polarisation and the `:y` axis)
- `coren` : Callable `coren(ω; z)` which returns the refractive index of the core
- `cladn` : Callable `cladn(ω; z)` which returns the refractive index of the cladding
- `model::Symbol=:full` : If `:full`, use the complete Marcatili model which takes into
                          account the dispersive influence of the cladding refractive index.
                          If `:reduced`, use the simplified model common in the literature
- `loss::Bool=true` : Whether to include loss.

"""
function MarcatiliMode(a, n, m, kind, ϕ, coren, cladn; model=:full, loss=true)
    # chkzkwarg makes sure that coren and cladn take z as a keyword argument
    aeff_intg = Aeff_Jintg(n, get_unm(n, m, kind), kind)
    MarcatiliMode(a, n, m, kind, get_unm(n, m, kind), ϕ,
                   chkzkwarg(coren), chkzkwarg(cladn),
                   model, Val(loss), aeff_intg)
end

"""
    MarcatiliMode(a, gas, P; kwargs...)

Create a MarcatiliMode for a capillary with radius `a` which is filled with `gas` to
pressure `P`.
"""
function MarcatiliMode(a, gas, P;
                        n=1, m=1, kind=:HE, ϕ=0.0, T=roomtemp, model=:full,
                        clad=:SiO2, loss=true)
    rfg = ref_index_fun(gas, P, T)
    rfs = ref_index_fun(clad)
    coren = (ω; z) -> rfg(wlfreq(ω))
    cladn = (ω; z) -> rfs(wlfreq(ω))
    MarcatiliMode(a, n, m, kind, ϕ, coren, cladn, model=model, loss=loss)
end

"""
    MarcatiliMode(a, gas, P, cladn; kwargs...)

Create a MarcatiliMode for a capillary made of a cladding material defined by the refractive
index `cladn(ω; z)` with a core radius `a` which is filled with `gas` to pressure `P`.
"""
function MarcatiliMode(a, gas, P, cladn;
                        n=1, m=1, kind=:HE, ϕ=0.0, T=roomtemp, model=:full, loss=true)
    rfg = ref_index_fun(gas, P, T)
    coren = (ω; z) -> rfg(wlfreq(ω))
    MarcatiliMode(a, n, m, kind, ϕ, coren, cladn, model=model, loss=loss)
end

"""
    MarcatiliMode(a, coren; kwargs...)

Create a MarcatiliMode for a capillary with radius `a` with `z`-dependent gas fill determined
by `coren(ω; z)`.
"""
function MarcatiliMode(a, coren;
                        n=1, m=1, kind=:HE, ϕ=0.0, model=:full, clad=:SiO2, loss=true)
    rfs = ref_index_fun(clad)
    cladn = (ω; z) -> rfs(wlfreq(ω))
    MarcatiliMode(a, n, m, kind, ϕ, coren, cladn, model=model, loss=loss)
end

"""
    MarcatiliMode(a; kwargs...)

Create a `MarcatiliMode` for a capillary with radius `a` and no gas fill.
"""
MarcatiliMode(a; kwargs...) = MarcatiliMode(a, (ω; z) -> 1; kwargs...)

"""
    neff(m::MarcatiliMode, ω; z=0)

Calculate the complex effective index of Marcatili mode with dielectric core and arbitrary
(metal or dielectric) cladding.

Adapted from:

Marcatili, E. & Schmeltzer, R.
"Hollow metallic and dielectric waveguides for long distance optical transmission and lasers
(Long distance optical transmission in hollow dielectric and metal circular waveguides,
examining normal mode propagation)."
Bell System Technical Journal 43, 1783–1809 (1964).
"""
function neff(m::MarcatiliMode, ω; z=0)
    εcl = m.cladn(ω, z=z)^2
    εco = m.coren(ω, z=z)^2
    vn = get_vn(εcl, m.kind)
    neff(m, ω, εco, vn, radius(m, z))
end

# Dispatch on loss to make neff type stable
# m.loss = Val{true}() (returns ComplexF64)
function neff(m::MarcatiliMode{Ta, Tco, Tcl, Val{true}}, ω, εco, vn, a) where {Ta, Tcl, Tco}
    if m.model == :full
        k = ω/c
        n = sqrt(complex(εco - (m.unm/(k*a))^2*(1 - im*vn/(k*a))^2))
        return (real(n) < 1e-3) ? (1e-3 + im*clamp(imag(n), 0, Inf)) : n
    elseif m.model == :reduced
        return ((1 + (εco - 1)/2 - c^2*m.unm^2/(2*ω^2*a^2))
                    + im*(c^3*m.unm^2)/(a^3*ω^3)*vn)
    else
        error("model must be :full or :reduced")
    end
end

# m.loss = Val{false}() (returns Float64)
function neff(m::MarcatiliMode{Ta, Tco, Tcl, Val{false}}, ω, εco, vn, a) where {Ta, Tcl, Tco}
    if m.model == :full
        k = ω/c
        n = real(sqrt(εco - (m.unm/(k*a))^2*(1 - im*vn/(k*a))^2))
        return (n < 1e-3) ? 1e-3 : n
    elseif m.model == :reduced
        return real(1 + (εco - 1)/2 - c^2*m.unm^2/(2*ω^2*a^2))
    else
        error("model must be :full or :reduced")
    end
end

function neff_wg(m::MarcatiliMode{Ta, Tco, Tcl, Val{true}}, ω; z=0) where {Ta, Tcl, Tco}
    εcl = m.cladn(ω, z=z)^2
    vn = get_vn(εcl, m.kind)
    a = radius(m, z)
    if m.model == :full
        k = ω/c
        return (m.unm/(k*a))^2*(1 - im*vn/(k*a))^2
    elseif m.model == :reduced
        return c^2*m.unm^2/(2*ω^2*a^2) + im*(c^3*m.unm^2)/(a^3*ω^3)*vn
    else
        error("model must be :full or :reduced")
    end
end

function neff_wg(m::MarcatiliMode{Ta, Tco, Tcl, Val{false}}, ω; z=0) where {Ta, Tcl, Tco}
    εcl = m.cladn(ω, z=z)^2
    vn = get_vn(εcl, m.kind)
    a = radius(m, z)
    if m.model == :full
        k = ω/c
        return (m.unm/(k*a))^2*(1 - im*vn/(k*a))^2
    elseif m.model == :reduced
        return c^2*m.unm^2/(2*ω^2*a^2)
    else
        error("model must be :full or :reduced")
    end
end

function neff(m::MarcatiliMode{Ta, Tco, Tcl, Val{true}}, εco, nwg) where {Ta, Tcl, Tco}
    if m.model == :full
        return sqrt(complex(εco - nwg))
    elseif m.model == :reduced
        return complex((1 + (εco - 1)/2 - nwg))
    else
        error("model must be :full or :reduced")
    end
end

function neff(m::MarcatiliMode{Ta, Tco, Tcl, Val{false}}, εco, nwg) where {Ta, Tcl, Tco}
    if m.model == :full
        return real(sqrt(complex(εco - nwg)))
    elseif m.model == :reduced
        return real((1 + (εco - 1)/2 - nwg))
    else
        error("model must be :full or :reduced")
    end
end

function get_vn(εcl, kind)
    if kind == :HE
        (εcl + 1)/(2*sqrt(complex(εcl - 1)))
    elseif kind == :TE
        1/sqrt(complex(εcl - 1))
    elseif kind == :TM
        εcl/sqrt(complex(εcl - 1))
    else
        error("kind must be :TE, :TM or :HE")
    end
end

function get_unm(n, m, kind)
    if (kind == :TE) || (kind == :TM)
        if (n != 0)
            error("n=0 for TE or TM modes")
        end
        sf_bessel_zero_Jnu(1, m)
    elseif kind == :HE
        sf_bessel_zero_Jnu(n-1, m)
    else
        error("kind must be :TE, :TM or :HE")
    end
end

radius(m::MarcatiliMode{<:Number, Tco, Tcl, LT}, z) where {Tcl, Tco, LT} = m.a
radius(m::MarcatiliMode, z) = m.a(z)

dimlimits(m::MarcatiliMode; z=0) = (:polar, (0.0, 0.0), (radius(m, z), 2π))

# we use polar coords, so xs = (r, θ)
function field(m::MarcatiliMode, xs; z=0)
    r, θ = xs
    if m.kind == :HE
        return (besselj(m.n-1, r*m.unm/radius(m, z)) .* SVector(
            cos(θ)*sin(m.n*(θ + m.ϕ)) - sin(θ)*cos(m.n*(θ + m.ϕ)),
            sin(θ)*sin(m.n*(θ + m.ϕ)) + cos(θ)*cos(m.n*(θ + m.ϕ))
            ))
    elseif m.kind == :TE
        return besselj(1, r*m.unm/radius(m, z)) .* SVector(-sin(θ), cos(θ))
    elseif m.kind == :TM
        return besselj(1, r*m.unm/radius(m, z)) .* SVector(cos(θ), sin(θ))
    end
end

function N(m::MarcatiliMode; z=0)
    np1 = (m.kind == :HE) ? m.n : 2
    π/2 * radius(m, z)^2 * besselj(np1, m.unm)^2 * sqrt(ε_0/μ_0)
end

function Aeff_Jintg(n, unm, kind)
    den, err = hquadrature(r -> r*besselj(n-1, unm*r)^4, 0, 1)
    np1 = (kind == :HE) ? n : 2
    num = 1/4 * besselj(np1, unm)^4
    return 2π*num/den
end

Aeff(m::MarcatiliMode; z=0) = radius(m, z)^2 * m.aeff_intg

"""
    gradient(gas, L, p0, p1)

Convenience function to create density and core index profiles for
simple two-point gradient fills defined by the waveguide length `L` and the pressures at
`z=0` and `z=L`.
"""
function gradient(gas, L, p0, p1)
    γ = sellmeier_gas(gas)
    dspl = densityspline(gas, Pmin=p0==p1 ? 0 : min(p0, p1), Pmax=max(p0, p1))
    p(z) =  z > L ? p1 :
            z <= 0 ? p0 :
            sqrt(p0^2 + z/L*(p1^2 - p0^2))
    dens(z) = dspl(p(z))
    coren(ω; z) = sqrt(1 + γ(wlfreq(ω)*1e6)*dens(z))
    return coren, dens
end

"""
    gradient(gas, Z, P)

Convenience function to create density and core index profiles for
multi-point gradient fills defined by positions `Z` and pressures `P`.
"""
function gradient(gas, Z, P)
    γ = sellmeier_gas(gas)
    ex = extrema(P)
    dspl = densityspline(gas, Pmin=ex[1]==ex[2] ? 0 : ex[1], Pmax=ex[2])
    function p(z)
        if z <= Z[1]
            return P[1]
        elseif z >= Z[end]
            return P[end]
        else
            i = findlast(x -> x < z, Z)
            return sqrt(P[i]^2 + (z - Z[i])/(Z[i+1] - Z[i])*(P[i+1]^2 - P[i]^2))
        end
    end
    dens(z) = dspl(p(z))
    coren(ω; z) = sqrt(1 + γ(wlfreq(ω)*1e6)*dens(z))
    return coren, dens
end

#= Avoid repeated calculation of the waveguide part of the effective index for modes with
    constant core radius.
    This is used by LinearOps.make_linop =#
function neff_β_grid(grid,
                   mode::MarcatiliMode{<:Number, Tco, Tcl, LT} where {Tco, Tcl, LT},
                   λ0)
    nwg = complex(zero(grid.ω))
    sidcs = (1:length(grid.ω))[grid.sidx]
    for iω in sidcs
        nwg[iω] = neff_wg(mode, grid.ω[iω]; z=0)
    end
    _neff = let nwg=nwg, ω=grid.ω, mode=mode
        _neff(iω; z) = neff(mode, mode.coren(ω[iω], z=z)^2, nwg[iω])
    end
    _β = let nwg=nwg, ω=grid.ω, _neff=_neff
        _β(iω; z) = ω[iω]/c*real(_neff(iω; z=z))
    end
    _neff, _β
end

# Collection of modes with fixed core radius
FixedCoreCollection = Union{
    Tuple{Vararg{MarcatiliMode{<:Number, Tco, Tcl, LT}} where {Tco, Tcl, LT}},
    AbstractArray{MarcatiliMode{<:Number, Tco, Tcl, LT} where {Tco, Tcl, LT}}
    }

function neff_grid(grid, modes::FixedCoreCollection, λ0; ref_mode=1)
    nwg = Array{ComplexF64, 2}(undef, (length(grid.ω), length(modes)))
    sidcs = (1:length(grid.ω))[grid.sidx]
    for (i, mi) in enumerate(modes)
        for iω in sidcs
            nwg[iω, i] = neff_wg(mi, grid.ω[iω]; z=0)
        end
    end
    _neff = let nwg=nwg, ω=grid.ω, modes=modes
        _neff(iω, iim; z) = neff(modes[iim], modes[iim].coren(ω[iω], z=z)^2, nwg[iω, iim])
    end
    _neff
end

"""
    transmission(a, λ, L; kind=:HE, n=1, m=1)

Calculate the transmission through a capillary with core radius `a` and length `L` at the
wavelength `λ` when propagating the `MarcatiliMode` defined by `kind`, `n` and `m`.
"""
function transmission(a, λ, L; kind=:HE, n=1, m=1)
    # TODO hardcoded fill needs to be updated if using absorbing materials
    mode = MarcatiliMode(a, :He, 0; n=n, m=m, kind=kind)
    Modes.transmission(mode, wlfreq(λ), L)
end

end


text:
function Base.show(io::IO, drv::Driver)::Nothing
    if drv.ptr == C_NULL
        print(io, "NULL Driver")
    else
        print(io, "Driver: $(shortname(drv))/$(longname(drv))")
    end
    return nothing
end

function Base.show(io::IO, dataset::AbstractDataset)::Nothing
    if dataset.ptr == C_NULL
        print(io, "NULL Dataset")
        return nothing
    end
    println(io, "GDAL Dataset ($(getdriver(dataset)))")
    println(io, "File(s): ")
    for (i, filename) in enumerate(filelist(dataset))
        println(io, "  $filename")
        if i > 5
            println(io, "  ...")
            break
        end
    end
    nrasters = nraster(dataset)
    if nrasters > 0
        print(io, "\nDataset (width x height): ")
        println(io, "$(width(dataset)) x $(height(dataset)) (pixels)")
        println(io, "Number of raster bands: $nrasters")
        for i in 1:min(nrasters, 3)
            print(io, "  ")
            summarize(io, getband(dataset, i))
        end
        nrasters > 3 && println(io, "  ...")
    end

    nlayers = nlayer(dataset)
    if nlayers > 0
        println(io, "\nNumber of feature layers: $nlayers")
        ndisplay = min(nlayers, 5) # display up to 5 layers
        for i in 1:ndisplay
            layer = getlayer(dataset, i - 1)
            layergeomtype = getgeomtype(layer)
            println(io, "  Layer $(i - 1): $(getname(layer)) ($layergeomtype)")
        end
        if nlayers > 5
            print(io, "  Remaining layers:\n    ")
            for i in 6:nlayers
                print(io, "$(getname(getlayer(dataset, i - 1))), ")
                # display up to 5 layer names per line
                if i % 5 == 0 && i < nlayers
                    print(io, "\n    ")
                end
            end
        end
    end
    return nothing
end

#Add method to avoid show from DiskArrays
Base.show(io::IO, raster::RasterDataset)::Nothing = show(io, raster.ds)

Base.show(io::IO, ::MIME"text/plain", raster::RasterDataset)::Nothing =
    show(io, raster.ds)

function summarize(io::IO, rasterband::AbstractRasterBand)::Nothing
    if rasterband.ptr == C_NULL
        print(io, "NULL RasterBand")
        return nothing
    end
    access = accessflag(rasterband)
    color = getname(getcolorinterp(rasterband))
    xsize = width(rasterband)
    ysize = height(rasterband)
    i = indexof(rasterband)
    pxtype = pixeltype(rasterband)
    println(io, "[$access] Band $i ($color): $xsize x $ysize ($pxtype)")
    return nothing
end

Base.show(io::IO, rasterband::AbstractRasterBand)::Nothing =
    show(io, "text/plain", rasterband)

function Base.show(
    io::IO,
    ::MIME"text/plain",
    rasterband::AbstractRasterBand,
)::Nothing
    summarize(io, rasterband)
    if rasterband.ptr == C_NULL
        return nothing
    end
    (x, y) = blocksize(rasterband)
    sc = getscale(rasterband)
    ofs = getoffset(rasterband)
    norvw = noverview(rasterband)
    ut = getunittype(rasterband)
    nv = getnodatavalue(rasterband)
    print(io, "    blocksize: $(x)×$(y), nodata: $nv, ")
    println(io, "units: $(sc)px + $ofs$ut")
    print(io, "    overviews: ")
    for i in 1:norvw
        ovr_band = getoverview(rasterband, i - 1)
        print(io, "($(i - 1)) $(width(ovr_band))x$(height(ovr_band)) ")
        i % 3 == 0 && print(io, "\n               ")
    end
    return nothing
end

# assumes that the layer is reset, and will reset it after display
function Base.show(io::IO, layer::AbstractFeatureLayer)::Nothing
    if layer.ptr == C_NULL
        print(io, "NULL FeatureLayer")
        return nothing
    end
    layergeomtype = getgeomtype(layer)
    println(io, "Layer: $(getname(layer))")
    featuredefn = layerdefn(layer)

    # Print Geometries
    n = ngeom(featuredefn)
    ngeomdisplay = min(n, 3)
    for i in 1:ngeomdisplay
        gfd = getgeomdefn(featuredefn, i - 1)
        display = "  Geometry $(i - 1) ($(getname(gfd))): [$(gettype(gfd))]"
        if length(display) > 75
            println(io, "$display[1:70]...")
            continue
        end
        if ngeomdisplay == 1 # only support printing of a single geom column
            for f in layer
                geomwkt = toWKT(getgeom(f))
                length(geomwkt) > 25 && (geomwkt = "$(geomwkt[1:20])...)")
                newdisplay = "$display, $geomwkt"
                if length(newdisplay) > 75
                    display = "$display, ..."
                    break
                else
                    display = newdisplay
                end
            end
        end
        println(io, display)
        resetreading!(layer)
    end
    n > 3 && println(io, "  ...\n  Number of Geometries: $n")

    # Print Features
    n = nfield(featuredefn)
    nfielddisplay = min(n, 5)
    for i in 1:nfielddisplay
        fd = getfielddefn(featuredefn, i - 1)
        display = "     Field $(i - 1) ($(getname(fd))): [$(gettype(fd))]"
        if length(display) > 75
            println(io, "$display[1:70]...")
            continue
        end
        for f in layer
            field = string(getfield(f, i - 1))
            length(field) > 25 && (field = "$(field[1:20])...")
            newdisplay = "$display, $field"
            if length(newdisplay) > 75
                display = "$display, ..."
                break
            else
                display = newdisplay
            end
        end
        println(io, display)
        resetreading!(layer)
    end
    n > 5 && print(io, "...\n Number of Fields: $n")
    return nothing
end

function Base.show(io::IO, featuredefn::AbstractFeatureDefn)::Nothing
    if featuredefn.ptr == C_NULL
        print(io, "NULL FeatureDefn")
        return nothing
    end
    n = ngeom(featuredefn)
    ngeomdisplay = min(n, 3)
    for i in 1:ngeomdisplay
        gfd = getgeomdefn(featuredefn, i - 1)
        println(io, "  Geometry (index $(i - 1)): $gfd")
    end
    n > 3 && println(io, "  ...\n  Number of Geometries: $n")

    n = nfield(featuredefn)
    nfielddisplay = min(n, 5)
    for i in 1:nfielddisplay
        fd = getfielddefn(featuredefn, i - 1)
        println(io, "     Field (index $(i - 1)): $fd")
    end
    n > 5 && print(io, "...\n Number of Fields: $n")
    return nothing
end

function Base.show(io::IO, fd::AbstractFieldDefn)::Nothing
    if fd.ptr == C_NULL
        print(io, "NULL FieldDefn")
        return nothing
    end
    print(io, "$(getname(fd)) ($(gettype(fd)))")
    return nothing
end

function Base.show(io::IO, gfd::AbstractGeomFieldDefn)::Nothing
    if gfd.ptr == C_NULL
        print(io, "NULL GeomFieldDefn")
        return nothing
    end
    print(io, "$(getname(gfd)) ($(gettype(gfd)))")
    return nothing
end

function Base.show(io::IO, feature::Feature)::Nothing
    if feature.ptr == C_NULL
        print(io, "NULL Feature")
        return nothing
    end
    println(io, "Feature")
    n = ngeom(feature)
    for i in 1:min(n, 3)
        displayname = geomname(getgeom(feature, i - 1))
        println(io, "  (index $(i - 1)) geom => $displayname")
    end
    n > 3 && println(io, "...\n Number of geometries: $n")
    n = nfield(feature)
    for i in 1:min(n, 10)
        displayname = getname(getfielddefn(feature, i - 1))
        print(io, "  (index $(i - 1)) $displayname => ")
        println(io, "$(getfield(feature, i - 1))")
    end
    n > 10 && print(io, "...\n Number of Fields: $n")
    return nothing
end

function Base.show(io::IO, spref::AbstractSpatialRef)::Nothing
    if spref.ptr == C_NULL
        print(io, "NULL Spatial Reference System")
        return nothing
    end
    projstr = toPROJ4(spref)
    if length(projstr) > 45
        projstart = projstr[1:35]
        projend = projstr[(end-4):end]
        print(io, "Spatial Reference System: $projstart ... $projend")
    else
        print(io, "Spatial Reference System: $projstr")
    end
    return nothing
end

function Base.show(io::IO, geom::AbstractGeometry)::Nothing
    if geom.ptr == C_NULL
        print(io, "NULL Geometry")
        return nothing
    end
    compact = get(io, :compact, false)

    if !compact
        print(io, "Geometry: ")
        geomwkt = toWKT(geom)
        if length(geomwkt) > 60
            print(io, "$(geomwkt[1:50]) ... $(geomwkt[(end - 4):end])")
        else
            print(io, "$geomwkt")
        end
    else
        print(io, "Geometry: $(getgeomtype(geom))")
    end
    return nothing
end

function Base.show(io::IO, ct::ColorTable)::Nothing
    if ct.ptr == C_NULL
        print(io, "NULL ColorTable")
        return nothing
    end
    palette = paletteinterp(ct)
    print(io, "ColorTable[$palette]")
    return nothing
end

Base.show(io::IO, ::MIME"text/plain", ct::ColorTable)::Nothing = show(io, ct)


text:
"""
Stores the solution of an assigment problem

Can be directly used to index the appropriate values from the cost matrix
"""
struct AssignmentSolution{T, S}
	row4col::Vector{Int}
	col4row::Vector{Int}
	cost::Union{T, Nothing}
	v::Vector{S}
	u::Vector{S}
	use_row::Bool

    # AssignmentSolution(row4col, col4row, cost::Nothing, u::Vector{S}, v::Vector{S}, use_row) where S = new{S, S}(row4col, col4row, cost, u, v, use_row)
    AssignmentSolution{T, S}(row4col, col4row, cost::T, u::Vector{S}, v::Vector{S}, use_row) where {T, S} = new{T, S}(row4col, col4row, cost, u, v, use_row)
    AssignmentSolution{T, S}(row4col, col4row, cost::Nothing, u::Vector{S}, v::Vector{S}, use_row) where {T, S} = new{T, S}(row4col, col4row, cost, u, v, use_row)
    AssignmentSolution(row4col, col4row, cost::T, u::Vector{S}, v::Vector{S}, use_row) where {T, S} = new{T, S}(row4col, col4row, cost, u, v, use_row)
end

function AssignmentSolution{T, S}(row4col, col4row, cost::Union{T, Nothing}, u::Vector{S}, v::Vector{S}) where {T, S}
	n = length(row4col)
	m = length(col4row)
	@assert n == length(v)
	@assert m == length(u)
	AssignmentSolution{T, S}(row4col, col4row, cost, u, v, n <=m)
end

AssignmentSolution(row4col::Vector{Int}, col4row::Vector{Int}, cost::T, u::Vector{S}, v::Vector{S}) where {T, S} =
    AssignmentSolution{isnothing(cost) ? S : T, S}(row4col, col4row, cost, u, v)

Base.adjoint(sol::AssignmentSolution{T, S}) where {T, S} = AssignmentSolution{T, S}(
	sol.col4row, sol.row4col, sol.cost, sol.u, sol.v, length(sol.col4row) <= length(sol.row4col))

Base.getindex(sol::AssignmentSolution, i::Integer) = sol.use_row ? CartesianIndex(i, sol.row4col[i]) : CartesianIndex(sol.col4row[i], i)

Base.eltype(sol::AssignmentSolution) = CartesianIndex{2}
Base.length(sol::AssignmentSolution) = sol.use_row ? length(sol.row4col) : length(sol.col4row)

Base.iterate(sol::AssignmentSolution) = Base.iterate(sol, 0)

function Base.iterate(sol::AssignmentSolution, i)
	i += 1
	i > length(sol) ? nothing : (sol[i], i)
end

Base.getindex(A::AbstractMatrix, sol::AssignmentSolution) = A[collect(sol)]

Base.show(io::IO, sol::AssignmentSolution) = if sol.use_row
    print(io, "AssignmentSolution(CartesianIndex.(1:$(length(sol.u)), $(sol.row4col)), $(sol.cost))")
else
    print(io, "AssignmentSolution(CartesianIndex.($(sol.col4row), 1:$(length(sol.v))), $(sol.cost))")
end


text:
include("cpu.jl")

p = read_program("inputs/21.txt")

function part_1()::CPU
    c = CPU([0,0,0,0,0,0])
    darr::Array{Int64} = []
    while true
        #push!(ips, c.registers[p.ip+1])
        try
            ins = p.lines[c.registers[p.ip+1]]
            pi = c.registers[p.ip+1]
            run_instruction(c, ins)
            c.registers[p.ip+1] += 1
            # if ins.c == 5
            #     pip = c.registers[p.ip+1]
            #     if pi > pip
            #         println("Jumping from " * string(pi) * " to " * string(pip))
            #         println(p.lines[pi])
            #         println(p.lines[pip])
            #         sleep(0.01)
            #     end
            # end
            if ins.b == 0 && ins.op != "seti"
                if c.registers[3] in darr
                    println(darr[end-5:end])
                    break
                else
                    push!(darr, c.registers[3])
                end
                if length(darr) % 100 == 0
                    println(length(darr))
                end
                #sleep(0.01)
            end
        catch
            return c
        end
    end
    return c
end

@time c = part_1()