text:
"""
    get_actions!(actions::Vector{A}, scene::EntityFrame{S,D,I}, roadway::R, models::Dict{I, M},) where {S,D,I,A,R,M<:DriverModel}
Fill in `actions` with the actions of each agent present in the scene. It calls `observe!`
and `rand` for each driver models.
`actions` will contain the actions to apply to update the state of each vehicle.
"""
function get_actions!(
    actions::Vector{A},
    scene::EntityFrame{S,D,I},
    roadway::R,
    models::Dict{I, M}, # id → model
    ) where {S,D,I,A,R,M<:DriverModel}

    for (i,veh) in enumerate(scene)
        model = models[veh.id]
        observe!(model, scene, roadway, veh.id)
        actions[i] = rand(model)
    end

    actions
end

"""
    tick!(scene::EntityFrame{S,D,I}, roadway::R, actions::Vector{A}, Δt::Float64) where {S,D,I,A,R}
update `scene` in place by updating the state of each vehicle given their current action in `actions`.
It calls the `propagate` method for each vehicle in the scene.
"""
function tick!(
    scene::EntityFrame{S,D,I},
    roadway::R,
    actions::Vector{A},
    Δt::Float64,
    ) where {S,D,I,A,R}

    for i in 1 : length(scene)
        veh = scene[i]
        state′ = propagate(veh, actions[i], roadway, Δt)
        scene[i] = Entity(state′, veh.def, veh.id)
    end

    return scene
end

"""
    reset_hidden_states!(models::Dict{Int,M}) where {M<:DriverModel}
reset hidden states of all driver models in `models`
"""
function reset_hidden_states!(models::Dict{Int,M}) where {M<:DriverModel}
    for model in values(models)
        reset_hidden_state!(model)
    end
    return models
end

"""
    simulate!(::Type{A}, rec::EntityQueueRecord{S,D,I}, scene::EntityFrame{S,D,I}, roadway::R, models::Dict{I,M}, nticks::Int) where {S,D,I,A,R,M<:DriverModel}
    simulate!(rec::EntityQueueRecord{S,D,I}, scene::EntityFrame{S,D,I}, roadway::R, models::Dict{I,M}, nticks::Int) where {S,D,I,R,M<:DriverModel}
Run nticks of simulation and place all nticks+1 scenes into the QueueRecord

    simulate!(::Type{A},rec::EntityQueueRecord{S,D,I}, scene::EntityFrame{S,D,I}, roadway::R, models::Dict{I,M}, nticks::Int, callbacks::C) where {S,D,I,A,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}
    simulate!(rec::EntityQueueRecord{S,D,I}, scene::EntityFrame{S,D,I}, roadway::R, models::Dict{I,M}, nticks::Int, callbacks::C) where {S,D,I,A,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}
Callback objects can also be passed in the simulate! function.

"""
function simulate!(
    ::Type{A},
    rec::EntityQueueRecord{S,D,I},
    scene::EntityFrame{S,D,I},
    roadway::R,
    models::Dict{I,M},
    nticks::Int,
    ) where {S,D,I,A,R,M<:DriverModel}

    empty!(rec)
    update!(rec, scene)
    actions = Array{A}(undef, length(scene))

    for tick in 1 : nticks
        println("\n\ntick: ", tick)
        get_actions!(actions, scene, roadway, models)
        tick!(scene, roadway, actions, rec.timestep)
        update!(rec, scene)
    end

    return rec
end
function simulate!(
    rec::EntityQueueRecord{S,D,I},
    scene::EntityFrame{S,D,I},
    roadway::R,
    models::Dict{I,M},
    nticks::Int,
    ) where {S,D,I,R,M<:DriverModel}

    return simulate!(Any, rec, scene, roadway, models, nticks)
end

"""
    Run a simulation and store the resulting scenes in the provided QueueRecord.
Only the ego vehicle is simulated; the other vehicles are as they were in the provided trajdata
Other vehicle states will be interpolated
"""
function simulate!(
    rec::EntityQueueRecord{S,D,I},
    model::DriverModel,
    egoid::I,
    trajdata::ListRecord{S,D,I},
    roadway::R,
    frame_start::Int,
    frame_end::Int;
    prime_history::Int=0, # no prime-ing
    scene::EntityFrame{S,D,I} =  allocate_frame(trajdata),
    ) where {S,D,I,R}

    @assert(isapprox(get_timestep(rec), get_timestep(trajdata)))

    # prime with history
    prime_with_history!(model, trajdata, roadway, frame_start, frame_end, egoid, scene)

    # add current frame
    update!(rec, get!(scene, trajdata, frame_start))
    observe!(model, scene, roadway, egoid)

    # run simulation
    frame_index = frame_start
    ego_veh = get_by_id(scene, egoid)
    while frame_index < frame_end

        # pull original scene
        get!(scene, trajdata, frame_index)

        # propagate ego vehicle and set
        ego_action = rand(model)
        ego_state = propagate(ego_veh, ego_action, roadway, get_timestep(rec))
        ego_veh = Entity(ego_veh, ego_state)
        scene[findfirst(ego_veh.id, scene)] = ego_veh

        # update record
        update!(rec, scene)

        # observe
        observe!(model, scene, roadway, ego_veh.id)

        # update time
        frame_index += 1
    end

    return rec
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 10/13/2014

# Reference: Chryssanthacopoulos, Kochenderfer - Collision Avoidance System Optimization with
# Probabilistic Pilot Response Models. American Control Conference. 2011.

module StochasticLinearPRImpl

export
    initialize,
    update,

    updatePilotResponse,

    StochasticLinearPR,
    StochasticLinearPRCommand,
    StochasticLinearPRRA

using AbstractPilotResponseImpl
using AbstractPilotResponseInterfaces
using CommonInterfaces

import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractPilotResponseInterfaces.updatePilotResponse

const statemap = Dict(:coc => 1, :first => 2, :multi => 3)
const RAmap = Dict(:none => 1, :climb => 2, :level => 3, :descend => 4)
const responsemap = Dict(:none => 1, :stay => 2, :follow => 3)
const strength_change_map = Dict(true => 1, false => 2)

function generateProbabilityTable()

  #key = (state,displayRA,response,currentRA,diffStrength)
  #value = (probability vector, output values vector)
  # where output values vector is an array of tuples (nextState,response) corresponding
  # to each probability given in probability vector
  #state = {:coc = clear of conflict, :first = first RA, :multi = change in RA}
  #displayRA = {:none = no RA, :climb, :level, :descend}
  #response = {:none = not following, :stay = following last displayRA, :follow = following latest displayRA}
  #currentRA = incoming RA, values are same as displayRA
  #diffStrength = {true=currentRA has same symbol but different target_rate, false=same target_Rate}
  A = Array(Tuple{Vector{Float64},Vector{Tuple{Symbol,Symbol}}}, length(statemap), length(RAmap), length(responsemap),
            length(RAmap), length(strength_change_map))

  # all entries going to :none go there with prob 1
  A[:, :, :, RAmap[:none], :] = ([1.0],[(:coc,:none)])

  # First RA
  for RA in map(x -> RAmap[x], [:climb, :level, :descend])
    #transition from coc to first
    A[statemap[:coc], RAmap[:none], responsemap[:none],:, strength_change_map[false]] = ([1/6,5/6],[(:first,:follow),(:first,:none)]) #first time presenting RA

    #staying in first
    A[statemap[:first], RA, responsemap[:none], RA, strength_change_map[false]] = ([1/6,5/6],[(:first,:follow),(:first,:none)]) #already active, but not following yet
    A[statemap[:first], RA, responsemap[:follow], RA, strength_change_map[false]] = ([1.0],[(:first,:follow)]) #already following
  end

  # Transition into multi
  for displayRA in map(x -> RAmap[x], [:climb, :level, :descend])
    for currentRA in map(x -> RAmap[x], [:climb, :level, :descend])
      for state in map(x -> statemap[x], [:first, :multi])
        b_same = displayRA == currentRA #same sense?
        A[state, displayRA, responsemap[:none], currentRA, strength_change_map[b_same]] = ([1/4,3/4],[(:multi,:follow),(:multi,:none)]) #change already active, but not following yet
        A[state, displayRA, responsemap[:stay], currentRA, strength_change_map[b_same]] = ([1/4,3/4],[(:multi,:follow),(:multi,:stay)]) #change already active, but not following yet
        A[state, displayRA, responsemap[:follow], currentRA, strength_change_map[b_same]] = ([1/4,3/4],[(:multi,:follow),(:multi,:stay)]) #change already active, but not following yet
      end
    end
  end

  # Staying in multi
  for RA in map(x -> RAmap[x], [:climb, :level, :descend])
    A[statemap[:multi], RA, responsemap[:none], RA, strength_change_map[false]] = ([1/4,3/4],[(:multi,:follow),(:multi,:none)]) #already active, but not following yet
    A[statemap[:multi], RA, responsemap[:stay], RA, strength_change_map[false]] = ([1/4,3/4],[(:multi,:follow),(:multi,:stay)]) #already active, but not following yet
    A[statemap[:multi], RA, responsemap[:follow], RA, strength_change_map[false]] = ([1.0],[(:multi,:follow)]) #already following
  end

  return A
end

const probabilityTable = generateProbabilityTable() #generate it once and store it

type StochasticLinearPRCommand

    t::Float64
    v_d::Float64
    h_d::Float64
    psi_d::Float64

    logProb::Float64 #log probability of generating this command
end

type StochasticLinearPRRA
  ra_active::Bool #true if RA is on
  target_rate::Float64
  dh_min::Float64 #min bound on h_d
  dh_max::Float64 #max bound on h_d
end

type StochasticLinearPR <: AbstractPilotResponse

  state::Symbol #[:coc,:first,:multi]
  displayRA::Symbol #RA currently displaying = [:none, :climb, :descend]
  target_rate::Union{Void,Float64} #target_rate of currently following, used for :stay
  response::Symbol #current pilot response = [:none, :stay, :follow]
  response_time::Int64 #for first RA
  output::StochasticLinearPRCommand

  probTable::Array{Tuple{Vector{Float64},Vector{Tuple{Symbol,Symbol}}}} #transition prob table

  function StochasticLinearPR()

    obj = new()
    obj.state = :coc
    obj.displayRA = :none
    obj.target_rate = nothing
    obj.response = :none
    obj.output = StochasticLinearPRCommand(0.0, 0.0, 0.0, 0.0, 0.0)
    obj.response_time = 0
    obj.probTable = probabilityTable

    return obj
  end
end

function updatePilotResponse(pr::StochasticLinearPR, update::StochasticLinearPRCommand, RA::StochasticLinearPRRA)

  t, v_d, h_d, psi_d = update.t, update.v_d, update.h_d, update.psi_d

  if RA.ra_active
    ra = RA.target_rate == 0 ? :level : RA.target_rate > 0 ? :climb : :descend
  else
    ra = :none
  end

  strengthening = (pr.displayRA == ra != :coc) && (nothing != pr.target_rate != RA.target_rate)
  probabilities, values = pr.probTable[statemap[pr.state], RAmap[pr.displayRA], responsemap[pr.response],
                                      RAmap[ra], strength_change_map[strengthening]]
  index = select_random(probabilities)
  pr.state,pr.response = values[index]
  pr.displayRA = ra

  if pr.displayRA != :none && pr.response == :none
    pr.response_time += 1
  end

  if pr.response == :stay
    h_d = pr.target_rate
  elseif pr.response == :follow
    h_d = pr.target_rate = RA.target_rate
  elseif pr.response == :none
    pr.target_rate = nothing
  else
    error("StochasticLinearPRImpl::updatePilotResponse: No such response!")
  end

  #Assign to output
  pr.output.t = t
  pr.output.v_d = v_d
  pr.output.h_d = h_d
  pr.output.psi_d = psi_d
  pr.output.logProb = log(probabilities[index])

  return pr.output
end

update(pr::StochasticLinearPR, command, RA) = update(pr,
                                                convert(StochasticLinearPRCommand, command),
                                                convert(StochasticLinearPRRA,RA))

update(pr::StochasticLinearPR,
     command::StochasticLinearPRCommand,
     RA::StochasticLinearPRRA) = updatePilotResponse(pr, command, RA)

function initialize(pr::StochasticLinearPR)
  pr.state = :coc
  pr.displayRA = :none
  pr.response = :none
  pr.response_time = 0
end

function select_random(weights)
# SELECT_RANDOM Randomly selects an index according to specified weights.
#   Returns a randomly selected index according to the distribution
#   specified by a vector of weights.
#
#   INDEX = SELECT_RANDOM(WEIGHTS) returns a scalar index INDEX selected
#   randomly according to the specified weights WEIGHTS represented as an
#   array.

    s = cumsum(weights)
    r = s[end] * rand()
    index = findfirst(x -> (x >= r), s)

    return index
end

end

ext: jl


text:
"""
Facet (`Vector{Int64}`) -> `BellGame`

    convert(
        ::Type{BellGame},
        facet::Vector{Int64},
        scenario::Union{BlackBox, LocalSignaling};
        rep = "normalized"::String
    )
"""
function convert( ::Type{BellGame},
    facet::Vector{Int64},
    scenario::Union{BlackBox,LocalSignaling};
    rep = "normalized"::String
)
    if !(rep in ("normalized", "generalized"))
        throw(DomainError(rep, "Argument `rep` must be either 'normalized' or 'generalized'"))
    end

    game_dims = strategy_dims(scenario)

    div_facet = facet .÷ gcd(facet...)

    bound = div_facet[end]
    game_matrix = (rep == "normalized") ? cat(
        reshape(div_facet[1:end-1], (game_dims[1]-1, game_dims[2])),
        zeros(Int64, (1,game_dims[2])),
        dims=1
    ) : reshape(div_facet[1:end-1], game_dims)

    BellGame(_reverse_game_normalization(game_matrix, bound)...)
end

"""
Facet (`Vector{Int64}`) -> `BellGame`

    convert(
        ::Type{BellGame},
        facet::Vector{Int64},
        scenario::BipartiteNonSignaling;
        rep = "non-signaling"::String
    )

Transforms LocalPolytope facets into `BellGame`  types.
"""
function convert(::Type{BellGame},
    facet::Vector{Int64},
    scenario::BipartiteNonSignaling;
    rep = "non-signaling"::String
)
    if !(rep in ["non-signaling","normalized","generalized"])
        throw(DomainError(rep, "input `rep` must be in [\"non-signaling\",\"normalized\",\"generalized\"]"))
    end

    game_dims = strategy_dims(scenario)
    game = (rep == "generalized") ? reshape(facet[1:end-1], game_dims) : zeros(Int64, game_dims)
    bound = facet[end]

    if rep == "non-signaling"
        α_dim = (scenario.A-1)*scenario.X
        β_dim = (scenario.B-1)*scenario.Y

        α_game = reshape(facet[1:α_dim], (scenario.A-1, scenario.X))
        β_game = reshape(facet[α_dim+1:α_dim+β_dim], (scenario.B-1, scenario.Y))
        αβ_game = reshape(facet[α_dim+β_dim+1:end-1], ((scenario.A-1)*(scenario.B-1), scenario.X*scenario.Y))
        αβ_col_sum = sum.(eachcol(αβ_game))

        # using non-signaling constraints to remove g_a,x
        for a in 1:scenario.A-1
            game[(a-1)*scenario.B+1:(a-1)*scenario.B + scenario.B-1,:] = αβ_game[(a-1)*(scenario.B-1)+1:a*(scenario.B-1),:]

            for x in 1:scenario.X
                if α_game[a,x] != 0
                    x_vec = zeros(Int64, scenario.X)
                    x_vec[x] = 1

                    y_vec = ones(Int64,scenario.Y)
                    αβ_col_id = findfirst(i -> i != 0, kron(x_vec,y_vec).*αβ_col_sum)

                    game[(a-1)*scenario.B+1:a*scenario.B,αβ_col_id] += α_game[a,x]*ones(Int64,scenario.B)
                end
            end
        end

        # using non-signaling constraints to remove g_b,y
        for b in 1:scenario.B-1
            game_row_ids = b:scenario.B:scenario.A*scenario.B-1

            for y in 1:scenario.Y
                if β_game[b,y] != 0
                    y_vec = zeros(Int64, scenario.Y)
                    y_vec[y] = 1

                    x_vec = ones(Int64, scenario.X)
                    αβ_col_id = findfirst(i -> i != 0, kron(x_vec,y_vec).*αβ_col_sum)

                    game[game_row_ids,αβ_col_id] += β_game[b,y]*ones(Int64,scenario.A)
                end
            end
        end
    elseif rep == "normalized"
        game[1:game_dims[1]-1,:] = reshape(facet[1:end-1], (game_dims[1]-1, game_dims[2]))
    end

    (game, bound) = _reverse_game_normalization(game,bound)

    BellGame(game,bound)
end

"""
`XPORTA.IEQ` to `BellGame`'s

    convert(
        ::Type{Vector{BellGame}},
        ieq::IEQ,
        scenario::Union{BlackBox, LocalSignaling};
        rep = "normalized" :: String
    )
"""
function convert(::Type{Vector{BellGame}},
    ieq::IEQ,
    scenario::Union{BlackBox,LocalSignaling};
    rep = "normalized"::String
)
    inequalities = convert.(Int64, ieq.inequalities)
    map( row_id -> convert(BellGame, inequalities[row_id,:], scenario, rep=rep), 1:size(inequalities,1))
end

"""
`BellGame` -> Facet (`Vector{Int64}`)

    convert(::Type{Vector{Int64}}, BG::BellGame; rep = "normalized" :: String)
"""
function convert(::Type{Vector{Int64}}, BG::BellGame; rep = "normalized"::String)
    if !(rep in ("normalized", "generalized"))
        throw(DomainError(rep, "Argument `rep` must be either 'normalized' or 'generalized'"))
    end

    bound = BG.β
    game_matrix = BG.game[:,:]

    game_dims = size(game_matrix)

    if rep == "normalized"
        (game_matrix, bound) = _apply_game_normalization!(game_matrix, bound)

        game_matrix = game_matrix[1:(end-1),:]
    end

    vcat(game_matrix[:], bound)
end

"""
BellGame -> Vector{Int64}

    convert(::Type{Vector{Int64}},
        BG::BellGame,
        scenario::BipartiteNonSignaling;
        rep = "non-signaling" :: String
    )

Transforms a `BellGame` for a `BipartiteNonSignaling` scenario into a facet vector.
"""
function convert(::Type{Vector{Int64}},
    BG::BellGame,
    scenario::BipartiteNonSignaling;
    rep = "non-signaling" :: String
)
    if !(rep in ["non-signaling","normalized","generalized"])
        throw(DomainError(rep, "input `rep` must be in [\"non-signaling\",\"normalized\",\"generalized\"]"))
    end

    game_dims = size(BG)
    v_dim = LocalPolytope.vertex_dims(scenario, rep)

    facet = (rep == "generalized") ? vcat(BG[:], BG.β) : zeros(Int64, v_dim+1)

    if rep == "normalized"
        (game_matrix, bound) = _apply_game_normalization(BG[:,:], BG.β)

        facet = vcat(game_matrix[1:game_dims[1]-1,:][:], bound)
    elseif rep == "non-signaling"
        (game_matrix, bound) = _apply_game_normalization!(BG[:,:], BG.β)

        # construct G(a|x) and G(b|y)
        # in each column, subtract off from each Alice/Bob column the values excluded from the non-signaling
        α_game = zeros(Int64, (scenario.A-1, scenario.X))
        β_game = zeros(Int64, (scenario.B-1, scenario.Y))

        # removing greatest output for Alice using non-signaling constraint
        for a in 1:scenario.A-1
            target_row = a * scenario.B
            subtract_vals = game_matrix[target_row,:]

            a_dims = (a-1)*scenario.B +1: a * scenario.B

            game_matrix[a_dims,:] -= ones(Int64, scenario.B) * subtract_vals'

            α_game_rows = map(x -> begin
                x_dims = (x-1)*scenario.Y+1:x*scenario.Y

                sum(subtract_vals[x_dims])
            end, 1:scenario.X)

            α_game[a,:] = α_game_rows
        end

        # removing greatest outputs for Bob using non-signaling constraint
        for b in 1:scenario.B-1

            target_row = (scenario.A-1) * (scenario.B) + b
            subtract_vals = game_matrix[target_row,:]

            b_dims = b:scenario.B:scenario.A * scenario.B -1

            game_matrix[b_dims,:] -= ones(Int64, scenario.A) * subtract_vals'

            β_game_rows = map(y -> begin
                y_dims = y:scenario.Y:scenario.X*scenario.Y

                sum(subtract_vals[y_dims])
            end, 1:scenario.Y)

            β_game[b,:] = β_game_rows
        end

        # All remaining terms are in the no-sig subspace and are taken as is
        αβ_game = zeros(Int64, ((scenario.A-1)*(scenario.B-1), scenario.X*scenario.Y))
        for a in 1:scenario.A-1
            αβ_game[(a-1)*(scenario.B-1) + 1:a*(scenario.B-1),:] = game_matrix[(a-1)*scenario.B+1:a*scenario.B-1,:]
        end

        facet = vcat(α_game[:], β_game[:], αβ_game[:], bound)
    end

    facet
end

"""
`BellGame`'s to `XPORTA.IEQ`

    convert(::Type{IEQ}, bell_games::Vector{BellGame}; rep = "normalized" :: String)
"""
function convert(::Type{IEQ}, bell_games::Vector{BellGame}; rep = "normalized"::String)
    ieq_vectors = map( bg -> convert(Vector{Int64}, bg, rep=rep), bell_games )
    IEQ(inequalities = hcat(ieq_vectors...)'[:,:])
end

"""

Applies the normalization constraint to remove all negative values in the provided
`game_matrix`. Returns a tuple `(new_game_matrix, new_bound)`
"""
function _reverse_game_normalization(game_matrix::Matrix{Int64}, bound::Int64) :: Tuple{Matrix{Int64}, Int64}
    new_bound = bound
    new_game_matrix = game_matrix

    for col_id in 1:size(game_matrix,2)
        col = game_matrix[:,col_id]
        col_min =  min(col...)

        if col_min != 0
            new_game_matrix[:,col_id] .-= col_min
            new_bound -= col_min
        end
    end

    (new_game_matrix, new_bound)
end

function _apply_game_normalization!(game_matrix::Matrix{Int64}, bound::Int64) :: Tuple{Matrix{Int64}, Int64}

    for col_id in 1:size(game_matrix,2)
        col = game_matrix[:,col_id]
        if col[end] !== 0
            game_matrix[:,col_id] .-= col[end]
            bound -= col[end]
        end
    end

    (game_matrix, bound)
end


text:
# # Conditional Maximum Likelihood for the Rasch Model
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [``rasch.ipynb``](@__NBVIEWER_ROOT_URL__examples/generated/rasch.ipynb)
#-
using Optim, Random #hide
#
# The Rasch model is used in psychometrics as a model for
# assessment data such as student responses to a standardized
# test. Let $X_{pi}$ be the response accuracy of student $p$
# to item $i$ where $X_{pi}=1$ if the item was answered correctly
# and $X_{pi}=0$ otherwise for $p=1,\ldots,n$ and $i=1,\ldots,m$.
# The model for this accuracy is
# ```math
#   P(\mathbf{X}_{p}=\mathbf{x}_{p}|\xi_p, \mathbf\epsilon) = \prod_{i=1}^m \dfrac{(\xi_p \epsilon_j)^{x_{pi}}}{1 + \xi_p\epsilon_i}
# ```
# where $\xi_p > 0$ the latent ability of person $p$ and $\epsilon_i > 0$
# is the difficulty of item $i$.

# We simulate data from this model:

Random.seed!(123)
n = 1000
m = 5
theta = randn(n)
delta = randn(m)
r = zeros(n)
s = zeros(m)

for i in 1:n
  p = exp.(theta[i] .- delta) ./ (1.0 .+ exp.(theta[i] .- delta))
  for j in 1:m
    if rand() < p[j] ##correct
      r[i] += 1
      s[j] += 1
    end
  end
end
f = [sum(r.==j) for j in 1:m];

# Since the number of parameters increases
# with sample size standard maximum likelihood will not provide us
# consistent estimates. Instead we consider the conditional likelihood.
# It can be shown that the Rasch model is an exponential family model and
# that the sum score $r_p = \sum_{i} x_{pi}$ is the sufficient statistic for
# $\xi_p$. If we condition on the sum score we should be able to eliminate
# $\xi_p$. Indeed, with a bit of algebra we can show
# ```math
# P(\mathbf{X}_p = \mathbf{x}_p | r_p, \mathbf\epsilon) = \dfrac{\prod_{i=1}^m \epsilon_i^{x{ij}}}{\gamma_{r_i}(\mathbf\epsilon)}
# ```
# where $\gamma_r(\mathbf\epsilon)$ is the elementary symmetric function of order $r$
# ```math
# \gamma_r(\mathbf\epsilon) = \sum_{\mathbf{y} : \mathbf{1}^\intercal \mathbf{y} = r} \prod_{j=1}^m \epsilon_j^{y_j}
# ```
# where the sum is over all possible answer configurations that give a sum
# score of $r$. Algorithms to efficiently compute $\gamma$ and its
# derivatives are available in the literature (see eg Baker (1996) for a review
# and Biscarri (2018) for a more modern approach)

function esf_sum!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  @inbounds for col in 1:n
    for r in 1:col
      row = col - r + 1
      S[row+1] = S[row+1] + x[col] * S[row]
    end
  end
end

function esf_ext!(S::AbstractArray{T,1}, H::AbstractArray{T,3}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  esf_sum!(S, x)
  H[:,:,1] .= zero(T)
  H[:,:,2] .= one(T)

  @inbounds for i in 3:n+1
    for j in 1:n
      H[j,j,i] = S[i-1] - x[j] * H[j,j,i-1]
      for k in j+1:n
        H[k,j,i] = S[i-1] - ((x[j]+x[k])*H[k,j,i-1] + x[j]*x[k]*H[k,j,i-2])
        H[j,k,i] = H[k,j,i]
      end
    end
  end
end

# The objective function we want to minimize is the negative log conditional
# likelihood
# ```math
# \begin{aligned}
# \log{L_C(\mathbf\epsilon|\mathbf{r})} &= \sum_{p=1}^n \sum_{i=1}^m x_{pi} \log{\epsilon_i} - \log{\gamma_{r_p}(\mathbf\epsilon)}\\
#   &= \sum_{i=1}^m s_i \log{\epsilon_i} - \sum_{r=1}^m f_r \log{\gamma_r(\mathbf\epsilon)}
# \end{aligned}
# ```
ϵ = ones(Float64, m)
β0 = zeros(Float64, m)
last_β = fill(NaN, m)
S = zeros(Float64, m+1)
H = zeros(Float64, m, m, m+1)

function calculate_common!(x, last_x)
  if x != last_x
    copyto!(last_x, x)
    ϵ .= exp.(-x)
    esf_ext!(S, H, ϵ)
  end
end
function neglogLC(β)
  calculate_common!(β, last_β)
  return -s'log.(ϵ) + f'log.(S[2:end])
end

# Parameter estimation is usually performed with respect to the unconstrained parameter
# $\beta_i = -\log{\epsilon_i}$. Taking the derivative with respect to $\beta_i$
# (and applying the chain rule) one obtains
# ```math
#   \dfrac{\partial\log L_C(\mathbf\epsilon|\mathbf{r})}{\partial \beta_i} = -s_i + \epsilon_i\sum_{r=1}^m \dfrac{f_r \gamma_{r-1}^{(j)}}{\gamma_r}
# ```
# where $\gamma_{r-1}^{(i)} = \partial \gamma_{r}(\mathbf\epsilon)/\partial\epsilon_i$.

function g!(storage, β)
  calculate_common!(β, last_β)
  for j in 1:m
    storage[j] = s[j]
    for l in 1:m
      storage[j] -= ϵ[j] * f[l] * (H[j,j,l+1] / S[l+1])
    end
  end
end

# Similarly the Hessian matrix can be computed
# ```math
#   \dfrac{\partial^2 \log L_C(\mathbf\epsilon|\mathbf{r})}{\partial \beta_i\partial\beta_j} = \begin{dcases} -\epsilon_i \sum_{r=1}^m \dfrac{f_r\gamma_{r-1}^{(i)}}{\gamma_r}\left(1 - \dfrac{\gamma_{r-1}^{(i)}}{\gamma_r}\right) & \text{if $i=j$}\\
#     -\epsilon_i\epsilon_j\sum_{r=1}^m \dfrac{f_r \gamma_{r-2}^{(i,j)}}{\gamma_r} - \dfrac{f_r\gamma_{r-1}^{(i)}\gamma_{r-1}^{(j)}}{\gamma_r^2} &\text{if $i\neq j$}
#    \end{dcases}
# ```
# where $\gamma_{r-2}^{(i,j)} = \partial^2 \gamma_{r}(\mathbf\epsilon)/\partial\epsilon_i\partial\epsilon_j$.

function h!(storage, β)
  calculate_common!(β, last_β)
  for j in 1:m
    for k in 1:m
      storage[k,j] = 0.0
      for l in 1:m
        if j == k
          storage[j,j] += f[l] * (ϵ[j]*H[j,j,l+1] / S[l+1]) *
            (1 - ϵ[j]*H[j,j,l+1] / S[l+1])
        elseif k > j
          storage[k,j] += ϵ[j] * ϵ[k] * f[l] *
            ((H[k,j,l] / S[l+1]) - (H[j,j,l+1] * H[k,k,l+1]) / S[l+1] ^ 2)
        else #k < j
          storage[k,j] += ϵ[j] * ϵ[k] * f[l] *
            ((H[j,k,l] / S[l+1]) - (H[j,j,l+1] * H[k,k,l+1]) / S[l+1] ^ 2)
        end
      end
    end
  end
end

# The estimates of the item parameters are then obtained via standard optimization
# algorithms (either Newton-Raphson or L-BFGS). One last issue is that the model is
# not identifiable (multiplying the $\xi_p$ by a constant and dividing the $\epsilon_i$
# by the same constant results in the same likelihood). Therefore some kind of constraint
# must be imposed when estimating the parameters. Typically either $\epsilon_1 = 0$ or
# $\prod_{i=1}^m \epsilon_i = 1$ (which is equivalent to $\sum_{i=1}^m \beta_i = 0$).

con_c!(c, x) = (c[1] = sum(x); c)
function con_jacobian!(J, x)
  J[1,:] .= ones(length(x))
end
function con_h!(h, x, λ)
    for i in 1:size(h)[1]
        for j in 1:size(h)[2]
            h[i,j] += (i == j) ? λ[1] : 0.0
        end
    end
end
lx = Float64[]; ux = Float64[]
lc = [0.0]; uc = [0.0]
df = TwiceDifferentiable(neglogLC, g!, h!, β0)
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!, lx, ux, lc, uc)
res = optimize(df, dfc, β0, IPNewton())

# Compare the estimate to the truth
delta_hat = res.minimizer
[delta delta_hat]


text:
const assetpath = "/"*relpath((@__FILE__)*"/../..","/") * "/assets/meshes/"
const mni_cortex = assetpath * "cortex/connectome-cortex.obj"
const fs_cortex = assetpath * "cortex/fs-cortex.obj"
const rh_cortex = assetpath * "cortex/rh-cortex.obj"
const lh_cortex = assetpath * "cortex/lh-cortex.obj"
const fs_cortex_mesh = load(fs_cortex)
const lh_cortex_mesh = load(lh_cortex)
const rh_cortex_mesh = load(rh_cortex)
const mni_cortex_mesh = load(mni_cortex)

function set_fig(;dimensions::Tuple{Int64, Int64}=(1600,900), view=:front)
    f = Figure(resolution = dimensions)
    ax = Axis3(f[1,1], aspect = :data, azimuth = View[view]pi, elevation=0.0pi)
    hidedecorations!(ax)
    hidespines!(ax)
    f
end

function get_hemisphere(parc, hemisphere::Symbol)
    ids = findall(x -> x == string(hemisphere), parc.Hemisphere)
    parc[ids,:ID]
end

function get_roi(parc::DataFrame, roi::String)
    findall(x -> occursin(roi, x), parc.Label)
end

Region = Dict(zip([:left, :right, :all, :connectome], [lh_cortex_mesh, rh_cortex_mesh, fs_cortex_mesh, mni_cortex_mesh]))
View = Dict(zip([:right, :front, :left, :back], [0.0, 0.5, 1.0, 1.5]))

function plot_cortex!(region::Symbol=:all; color=(:grey,0.05), transparency::Bool=true, kwargs...)
    mesh!(Region[region], color=color, transparency=transparency, kwargs...)
end

function plot_cortex(region::Symbol=:all; view=:left, color=(:grey,1.0), transparency::Bool=false, kwargs...)
    f = set_fig(view=view)
    plot_cortex!(region; color, transparency, kwargs...)
    f
end

function plot_parc!(connectome::Connectome, hemisphere::Symbol; alpha=1.0)
    h_ids = get_hemisphere(connectome.parc, hemisphere)

    colors = distinguishable_colors(length(h_ids))
    for (i, j) in enumerate(h_ids)
        roi = load(assetpath * "DKT/roi_$(j).obj")
        mesh!(roi, color=(colors[i], alpha), transparency=false, show_axis=false)
    end
end

function plot_parc(connectome::Connectome, hemisphere::Symbol; alpha=1.0, view=:left)
    f = set_fig(view=view)
    plot_parc!(connectome, hemisphere; alpha)
    f
end

function plot_parc(connectome::Connectome; alpha=1.0)
    f = set_fig()
    plot_parc!(connectome, :left; alpha)
    plot_parc!(connectome, :right; alpha)
    f
end

function plot_roi!(roi::Int, color=(:grey,1.0))
    meshpath = assetpath * "DKT/roi_$(roi).obj"
    mesh!(load(meshpath), color=color, transparency=false)
end

function plot_roi!(roi::Vector{Int64}, color)
    for i in roi
        plot_roi!(i, color)
    end
end

function plot_roi(connectome::Connectome, roi::String; cortexcolor=(:grey,0.05), color=(:blue,0.1), transparency=true)

    f = set_fig()
    plot_cortex!(:all; color=cortexcolor, transparency=transparency)

    ID = get_roi(connectome.parc, roi)

    for i in ID
        plot_roi!(i, color)
    end
    f
end

function plot_roi(connectome::Connectome, roi::String, hemisphere::Symbol; cortexcolor=(:grey,0.05), color=(:blue, 0.1), roi_alpha=1.0, transparency=true)

    f = set_fig()
    plot_cortex!(hemisphere; color=cortexcolor, transparency=transparency)
    plot_cortex!
    ID = get_roi(connectome.parc, roi)
    h_ID = get_hemisphere(connectome.parc[ID,:], hemisphere)
    for j in h_ID
        plot_roi!(j, color)
    end
    f
end

function plot_roi(connectome::Connectome, roi::Vector{String}; cortexcolor=(:grey,0.05), roi_alpha=1.0, transparency=true)

    f = set_fig()
    plot_cortex!(:all; color=cortexcolor, transparency=transparency)
    color = distinguishable_colors(length(roi))
    for (i, j) in enumerate(roi)
        ID = get_roi(connectome.parc, j)
        for k in ID
            plot_roi!(k, (color[i], roi_alpha))
        end
    end
    f
end

function plot_roi(connectome::Connectome, roi::Vector{String}, hemisphere::Symbol; cortexcolor=(:grey,0.05), roi_alpha=1.0, transparency=true)

    f = set_fig()
    plot_cortex!(:all; color=cortexcolor, transparency=transparency)
    color = distinguishable_colors(length(roi))
    for (i, j) in enumerate(roi)
        ID = get_roi(connectome.parc, j)
        h_ID = get_hemisphere(connectome.parc[ID,:], hemisphere)
        for k in h_ID
            plot_roi!(k, (color[i], roi_alpha))
        end
    end
    f
end

function plot_vertex!(connectome::Connectome, node_size, color, transparency::Bool=true)
    x, y, z = connectome.parc.x[:], connectome.parc.y[:], connectome.parc.z[:]
    meshscatter!(x, y, z, markersize=node_size, color=color, transparency=transparency)
end

function plot_vertex(connectome::Connectome; node_size=1.0, color=(:blue,0.5))
    f = set_fig()
    plot_cortex!(:connectome)
    plot_vertex!(connectome, node_size, color)
    f
end

function plot_edges!(connectome::Connectome, color)
    x, y, z = connectome.parc.x[:], connectome.parc.y[:], connectome.parc.z[:]
    coordindex = findall(x->x>0, LowerTriangular(connectome.A))

    for i ∈ 1:length(coordindex)
        j, k = coordindex[i][1], coordindex[i][2]
        weight = connectome.A[j, k]
        lines!(x[[j,k]], y[[j,k]], z[[j,k]],
               color=get(color, weight), #matter
               linewidth=clamp(20*weight,2,50),
               transparency=true)
    end
end

function plot_connectome(connectome::Connectome; node_size=1.0, node_color=(:blue, 0.5), edge_color)
    f = set_fig()
    plot_cortex!(:connectome)
    plot_vertex!(connectome, node_size, node_color)
    plot_edges!(connectome, edge_color)
    f
end


text:
# This file is supposed to do handle all specifications for separate theories like for instance `su2higgs()` sets all structs and functions to default which are needed for this particular theory.
# Not calling any theory should give warnings when ambigious functions are called

## General Theory type and default struct
abstract type AbstractTheory end

struct NoTheory <: AbstractTheory end

# for converting type to symbol and back
const _theoryType = Dict{Symbol, DataType}(:none => NoTheory)
const _theorySymbol = Dict{DataType, Symbol}(NoTheory => :none)
# supported theoris
const _theories = Symbol[]
const _initialized_theories = Set{Symbol}()
const _default_theories = (:none)

# Helper functions to get list of supported theories and current theories
"Returns a list of supported Theorys"
theories() = _theories

"Returns the name of the current Theory"
theory_name() = CURRENT_THEORY.sym

# Returns type of theory if supported
function _theory_instance(sym::Symbol)::AbstractTheory
    haskey(_theoryType, sym) ? _theoryType[sym]() : error("Unsupported theory $sym")
end

# macro for initializing new theories
macro init_theory(s)
    package_str = string(s)
    str = lowercase(package_str)
    sym = Symbol(str)
    T = Symbol(string(s) * "Theory")
    esc(quote
        struct $T <: AbstractTheory end
        export $sym
        $sym(; kw...) = (#=default(; reset = false, kw...);=# theory($T()))
        theory_name(::$T) = Symbol($str)
        push!(_theories, Symbol($str))
        _theoryType[Symbol($str)] = $T
        _theorySymbol[$T] = Symbol($str)
    end)
end

## Struct for currently activated theory
mutable struct CurrentTheory
    sym::Symbol
    theory::AbstractTheory
end
# Constructor
CurrentTheory(sym::Symbol) = CurrentTheory(sym,_theory_instance(sym))

## Default Theory
_fallback_default_theory() = theory(:none)

function _pick_default_theory()
    env_default = get(ENV, "LATSPEC_DEFAULT_THEORY", "")
    if env_default != ""
        sym = Symbol(lowercase(env_default))
        if sym in _theories
            theory(sym)
        else
            @warn("You have set LATSPEC_DEFAULT_THEORY=$env_default but it is not a valid theory.  Choose from:\n\t" *
                 join(sort(_theories), "\n\t"))
            _fallback_default_theory()
        end
    else
        _fallback_default_theory()
    end
end

## Function overloads of theory
"""
Returns the current theory used for spectroscopy.  Initializes package on first call.
"""
function theory()
    if CURRENT_THEORY.sym == :none
        _pick_default_theory()
    end

    CURRENT_THEORY.theory
end

"""
Set the Spectroscopy Theory.
"""
function theory(theory::AbstractTheory)
    sym = theory_name(theory)
    if !(sym in _initialized_theories)
        _initialize_theory(theory)
        push!(_initialized_theories, sym)
    end
    CURRENT_THEORY.sym = sym
    CURRENT_THEORY.theory = theory
    theory
end

function theory(sym::Symbol)
    if sym in _theories
        theory(_theory_instance(sym))
    else
        @warn("`:$sym` is not a supported theory.")
        theory()
    end
end

## Init supported theories
@init_theory none
@init_theory su2higgs
@init_theory su3higgs
@init_theory sp4

# initialize the theory
function _initialize_theory(theory::AbstractTheory)
    sym = theory_name(theory)
    @eval Main begin
        println($sym, " has been initialized")
    end
end


text:
using Test

@testset "AtomsBase.jl" begin
    include("interface.jl")
    include("fast_system.jl")
    include("atom.jl")
    include("properties.jl")
    include("printing.jl")
end


text:
# an AbstractIndex is a thing that can be used to look up ordered things by name, but that
# will also accept a position or set of positions or range or other things and pass them
# through cleanly.
# an Index is the usual implementation.
# a SimpleIndex only works if the things are integer indexes, which is weird.
abstract AbstractIndex

type Index <: AbstractIndex   # an OrderedDict would be nice here...
    lookup::Dict{Symbol, Int}      # name => names array position
    names::Vector{Symbol}
end
function Index(names::Vector{Symbol}; allow_duplicates=true)
    u = make_unique(names, allow_duplicates=allow_duplicates)
    lookup = Dict{Symbol, Int}(zip(u, 1:length(u)))
    Index(lookup, u)
end
Index() = Index(Dict{Symbol, Int}(), Symbol[])
Base.length(x::Index) = length(x.names)
Base.names(x::Index) = copy(x.names)
_names(x::Index) = x.names
Base.copy(x::Index) = Index(copy(x.lookup), copy(x.names))
Base.deepcopy(x::Index) = copy(x) # all eltypes immutable
Base.isequal(x::Index, y::Index) = isequal(x.lookup, y.lookup) && isequal(x.names, y.names)
# Imported in DataFrames.jl for compatibility across Julia 0.4 and 0.5
(==)(x::Index, y::Index) = isequal(x, y)

function names!(x::Index, nms::Vector{Symbol}; allow_duplicates=false)
    if length(nms) != length(x)
        throw(ArgumentError("Length of nms doesn't match length of x."))
    end
    newindex = Index(nms, allow_duplicates=allow_duplicates)
    x.names = newindex.names
    x.lookup = newindex.lookup
    return x
end

function rename!(x::Index, nms)
    for (from, to) in nms
        if haskey(x, to)
            error("Tried renaming $from to $to, when $to already exists in the Index.")
        end
        x.lookup[to] = col = pop!(x.lookup, from)
        x.names[col] = to
    end
    return x
end

rename!(x::Index, from, to) = rename!(x, zip(from, to))
rename!(x::Index, from::Symbol, to::Symbol) = rename!(x, ((from, to),))
rename!(x::Index, f::Function) = rename!(x, [(x,f(x)) for x in x.names])
rename!(f::Function, x::Index) = rename!(x, f)

rename(x::Index, args...) = rename!(copy(x), args...)
rename(f::Function, x::Index) = rename(x, f)

Base.haskey(x::Index, key::Symbol) = haskey(x.lookup, key)
Base.haskey(x::Index, key::Real) = 1 <= key <= length(x.names)
Base.keys(x::Index) = names(x)

# TODO: If this should stay 'unsafe', perhaps make unexported
function Base.push!(x::Index, nm::Symbol)
    x.lookup[nm] = length(x) + 1
    push!(x.names, nm)
    return x
end

function Base.merge!(x::Index, y::Index)
    adds = add_names(x, y)
    i = length(x)
    for add in adds
        i += 1
        x.lookup[add] = i
    end
    append!(x.names, adds)
    return x
end

Base.merge(x::Index, y::Index) = merge!(copy(x), y)

function Base.delete!(x::Index, idx::Integer)
    # reset the lookup's beyond the deleted item
    for i in (idx + 1):length(x.names)
        x.lookup[x.names[i]] = i - 1
    end
    delete!(x.lookup, x.names[idx])
    deleteat!(x.names, idx)
    return x
end

function Base.delete!(x::Index, nm::Symbol)
    if !haskey(x.lookup, nm)
        return x
    end
    idx = x.lookup[nm]
    return delete!(x, idx)
end

function Base.empty!(x::Index)
    empty!(x.lookup)
    empty!(x.names)
    x
end

function Base.insert!(x::Index, idx::Integer, nm::Symbol)
    1 <= idx <= length(x.names)+1 || error(BoundsError())
    for i = idx:length(x.names)
        x.lookup[x.names[i]] = i + 1
    end
    x.lookup[nm] = idx
    insert!(x.names, idx, nm)
    x
end

Base.getindex(x::Index, idx::Symbol) = x.lookup[idx]
Base.getindex(x::AbstractIndex, idx::Real) = @compat Int(idx)
Base.getindex(x::AbstractIndex, idx::AbstractVector{Nullable{Bool}}) =
    getindex(x, convert(Vector{Bool}, idx, false))
Base.getindex{T<:Nullable}(x::AbstractIndex, idx::AbstractVector{T}) =
    getindex(x, dropnull(idx))
Base.getindex(x::AbstractIndex, idx::AbstractVector{Bool}) = find(idx)
Base.getindex(x::AbstractIndex, idx::Range) = [idx;]
Base.getindex{T <: Real}(x::AbstractIndex, idx::AbstractVector{T}) = convert(Vector{Int}, idx)
Base.getindex(x::AbstractIndex, idx::AbstractVector{Symbol}) = [x.lookup[i] for i in idx]

type SimpleIndex <: AbstractIndex
    length::Integer
end
SimpleIndex() = SimpleIndex(0)
Base.length(x::SimpleIndex) = x.length
Base.names(x::SimpleIndex) = nothing
_names(x::SimpleIndex) = nothing

# Helpers

function add_names(ind::Index, add_ind::Index)
    u = names(add_ind)

    seen = Set(_names(ind))
    dups = Int[]

    for i in 1:length(u)
        name = u[i]
        in(name, seen) ? push!(dups, i) : push!(seen, name)
    end
    for i in dups
        nm = u[i]
        k = 1
        while true
            newnm = Symbol("$(nm)_$k")
            if !in(newnm, seen)
                u[i] = newnm
                push!(seen, newnm)
                break
            end
            k += 1
        end
    end

    return u
end


text:
using StatsFuns,Distributions
import ForwardDiff, Zygote

f(x) = poislogpdf(x[1], x[2])
#new Poisson logpdf
poission_lpdf(x) = -x[1]+x[2]*log(x[1])-log(factorial(x[2]))

x = [.5,2.]
f(x) |> display

fd_grad1 =ForwardDiff.gradient(f, x)      # works
fd_grad1 |> display

try
  zd_grade1=Zygote.gradient(f, x)           #fails
catch e
  println(e)
end

poission_lpdf(x)
fd_grad2=ForwardDiff.gradient(poission_lpdf, x)      # works
fd_grad2 |> display

zd_grad2=Zygote.gradient(poission_lpdf, x)           #works
zd_grad2 |> display

fd_grad1 == fd_grad2 == zd_grad2[1]


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

JLLWrappers.@generate_main_file_header("GeographicLibWrapper")
JLLWrappers.@generate_main_file("GeographicLibWrapper", UUID("73ec4600-8e52-54a9-acd2-2fa891c2bd91"))
end  # module GeographicLibWrapper_jll


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: ivs

using AWS.Compat
using AWS.UUIDs
"""
    BatchGetChannel()

Performs GetChannel on multiple ARNs simultaneously.

# Required Parameters
- `arns`: Array of ARNs, one per channel.

"""

batch_get_channel(arns; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/BatchGetChannel", Dict{String, Any}("arns"=>arns); aws_config=aws_config)
batch_get_channel(arns, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/BatchGetChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arns"=>arns), args)); aws_config=aws_config)

"""
    BatchGetStreamKey()

Performs GetStreamKey on multiple ARNs simultaneously.

# Required Parameters
- `arns`: Array of ARNs, one per channel.

"""

batch_get_stream_key(arns; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/BatchGetStreamKey", Dict{String, Any}("arns"=>arns); aws_config=aws_config)
batch_get_stream_key(arns, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/BatchGetStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arns"=>arns), args)); aws_config=aws_config)

"""
    CreateChannel()

Creates a new channel and an associated stream key to start streaming.

# Optional Parameters
- `authorized`: Whether the channel is authorized. Default: false.
- `latencyMode`: Channel latency mode. Default: LOW.
- `name`: Channel name.
- `tags`: See Channeltags.
- `type`: Channel type, which determines the allowable resolution and bitrate. If you exceed the allowable resolution or bitrate, the stream probably will disconnect immediately. Valid values:    STANDARD: Multiple qualities are generated from the original input, to automatically give viewers the best experience for their devices and network conditions. Vertical resolution can be up to 1080 and bitrate can be up to 8.5 Mbps.    BASIC: Amazon IVS delivers the original input to viewers. The viewer’s video-quality choice is limited to the original input. Vertical resolution can be up to 480 and bitrate can be up to 1.5 Mbps.   Default: STANDARD.
"""

create_channel(; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/CreateChannel"; aws_config=aws_config)
create_channel(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/CreateChannel", args; aws_config=aws_config)

"""
    CreateStreamKey()

Creates a stream key, used to initiate a stream, for the specified channel ARN. Note that CreateChannel creates a stream key. If you subsequently use CreateStreamKey on the same channel, it will fail because a stream key already exists and there is a limit of 1 stream key per channel. To reset the stream key on a channel, use DeleteStreamKey and then CreateStreamKey.

# Required Parameters
- `channelArn`: ARN of the channel for which to create the stream key.

# Optional Parameters
- `tags`: See Channeltags.
"""

create_stream_key(channelArn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/CreateStreamKey", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
create_stream_key(channelArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/CreateStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), args)); aws_config=aws_config)

"""
    DeleteChannel()

Deletes the specified channel and its associated stream keys.

# Required Parameters
- `arn`: ARN of the channel to be deleted.

"""

delete_channel(arn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/DeleteChannel", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_channel(arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/DeleteChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), args)); aws_config=aws_config)

"""
    DeletePlaybackKeyPair()

Deletes a specified authorization key pair. This invalidates future viewer tokens generated using the key pair’s privateKey.

# Required Parameters
- `arn`: ARN of the key pair to be deleted.

"""

delete_playback_key_pair(arn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/DeletePlaybackKeyPair", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_playback_key_pair(arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/DeletePlaybackKeyPair", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), args)); aws_config=aws_config)

"""
    DeleteStreamKey()

Deletes the stream key for the specified ARN, so it can no longer be used to stream.

# Required Parameters
- `arn`: ARN of the stream key to be deleted.

"""

delete_stream_key(arn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/DeleteStreamKey", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_stream_key(arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/DeleteStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), args)); aws_config=aws_config)

"""
    GetChannel()

Gets the channel configuration for the specified channel ARN. See also BatchGetChannel.

# Required Parameters
- `arn`: ARN of the channel for which the configuration is to be retrieved.

"""

get_channel(arn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetChannel", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_channel(arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), args)); aws_config=aws_config)

"""
    GetPlaybackKeyPair()

Gets a specified playback authorization key pair and returns the arn and fingerprint. The privateKey held by the caller can be used to generate viewer authorization tokens, to grant viewers access to authorized channels.

# Required Parameters
- `arn`: ARN of the key pair to be returned.

"""

get_playback_key_pair(arn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetPlaybackKeyPair", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_playback_key_pair(arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetPlaybackKeyPair", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), args)); aws_config=aws_config)

"""
    GetStream()

Gets information about the active (live) stream on a specified channel.

# Required Parameters
- `channelArn`: Channel ARN for stream to be accessed.

"""

get_stream(channelArn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetStream", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
get_stream(channelArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetStream", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), args)); aws_config=aws_config)

"""
    GetStreamKey()

Gets stream-key information for a specified ARN.

# Required Parameters
- `arn`: ARN for the stream key to be retrieved.

"""

get_stream_key(arn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetStreamKey", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_stream_key(arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/GetStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), args)); aws_config=aws_config)

"""
    ImportPlaybackKeyPair()

Imports the public portion of a new key pair and returns its arn and fingerprint. The privateKey can then be used to generate viewer authorization tokens, to grant viewers access to authorized channels.

# Required Parameters
- `publicKeyMaterial`: The public portion of a customer-generated key pair.

# Optional Parameters
- `name`: An arbitrary string (a nickname) assigned to a playback key pair that helps the customer identify that resource. The value does not need to be unique.
- `tags`: Any tags provided with the request are added to the playback key pair tags.
"""

import_playback_key_pair(publicKeyMaterial; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ImportPlaybackKeyPair", Dict{String, Any}("publicKeyMaterial"=>publicKeyMaterial); aws_config=aws_config)
import_playback_key_pair(publicKeyMaterial, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ImportPlaybackKeyPair", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("publicKeyMaterial"=>publicKeyMaterial), args)); aws_config=aws_config)

"""
    ListChannels()

Gets summary information about all channels in your account, in the AWS region where the API request is processed. This list can be filtered to match a specified string.

# Optional Parameters
- `filterByName`: Filters the channel list to match the specified name.
- `maxResults`: Maximum number of channels to return.
- `nextToken`: The first channel to retrieve. This is used for pagination; see the nextToken response field.
"""

list_channels(; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListChannels"; aws_config=aws_config)
list_channels(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListChannels", args; aws_config=aws_config)

"""
    ListPlaybackKeyPairs()

Gets summary information about playback key pairs.

# Optional Parameters
- `maxResults`: The first key pair to retrieve. This is used for pagination; see the nextToken response field.
- `nextToken`: Maximum number of key pairs to return.
"""

list_playback_key_pairs(; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListPlaybackKeyPairs"; aws_config=aws_config)
list_playback_key_pairs(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListPlaybackKeyPairs", args; aws_config=aws_config)

"""
    ListStreamKeys()

Gets summary information about stream keys for the specified channel.

# Required Parameters
- `channelArn`: Channel ARN used to filter the list.

# Optional Parameters
- `maxResults`: Maximum number of streamKeys to return.
- `nextToken`: The first stream key to retrieve. This is used for pagination; see the nextToken response field.
"""

list_stream_keys(channelArn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListStreamKeys", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
list_stream_keys(channelArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListStreamKeys", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), args)); aws_config=aws_config)

"""
    ListStreams()

Gets summary information about live streams in your account, in the AWS region where the API request is processed.

# Optional Parameters
- `maxResults`: Maximum number of streams to return.
- `nextToken`: The first stream to retrieve. This is used for pagination; see the nextToken response field.
"""

list_streams(; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListStreams"; aws_config=aws_config)
list_streams(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/ListStreams", args; aws_config=aws_config)

"""
    ListTagsForResource()

Gets information about AWS tags for the specified ARN.

# Required Parameters
- `resourceArn`: The ARN of the resource to be retrieved.

# Optional Parameters
- `maxResults`: Maximum number of tags to return.
- `nextToken`: The first tag to retrieve. This is used for pagination; see the nextToken response field.
"""

list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = ivs("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    PutMetadata()

Inserts metadata into an RTMPS stream for the specified channel. A maximum of 5 requests per second per channel is allowed, each with a maximum 1KB payload.

# Required Parameters
- `channelArn`: ARN of the channel into which metadata is inserted. This channel must have an active stream.
- `metadata`: Metadata to insert into the stream. Maximum: 1 KB per request.

"""

put_metadata(channelArn, metadata; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/PutMetadata", Dict{String, Any}("channelArn"=>channelArn, "metadata"=>metadata); aws_config=aws_config)
put_metadata(channelArn, metadata, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/PutMetadata", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn, "metadata"=>metadata), args)); aws_config=aws_config)

"""
    StopStream()

Disconnects the incoming RTMPS stream for the specified channel. Can be used in conjunction with DeleteStreamKey to prevent further streaming to a channel.  Many streaming client-software libraries automatically reconnect a dropped RTMPS session, so to stop the stream permanently, you may want to first revoke the streamKey attached to the channel.

# Required Parameters
- `channelArn`: ARN of the channel for which the stream is to be stopped.

"""

stop_stream(channelArn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/StopStream", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
stop_stream(channelArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/StopStream", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), args)); aws_config=aws_config)

"""
    TagResource()

Adds or updates tags for the AWS resource with the specified ARN.

# Required Parameters
- `resourceArn`: ARN of the resource for which tags are to be added or updated.
- `tags`: Array of tags to be added or updated.

"""

tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes tags from the resource with the specified ARN.

# Required Parameters
- `resourceArn`: ARN of the resource for which tags are to be removed.
- `tagKeys`: Array of tags to be removed.

"""

untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = ivs("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateChannel()

Updates a channel's configuration. This does not affect an ongoing stream of this channel. You must stop and restart the stream for the changes to take effect.

# Required Parameters
- `arn`: ARN of the channel to be updated.

# Optional Parameters
- `authorized`: Whether the channel is authorized. Default: false.
- `latencyMode`: Channel latency mode. Default: LOW.
- `name`: Channel name.
- `type`: Channel type, which determines the allowable resolution and bitrate. If you exceed the allowable resolution or bitrate, the stream probably will disconnect immediately. Valid values:    STANDARD: Multiple qualities are generated from the original input, to automatically give viewers the best experience for their devices and network conditions. Vertical resolution can be up to 1080 and bitrate can be up to 8.5 Mbps.    BASIC: Amazon IVS delivers the original input to viewers. The viewer’s video-quality choice is limited to the original input. Vertical resolution can be up to 480 and bitrate can be up to 1.5 Mbps.   Default: STANDARD.
"""

update_channel(arn; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/UpdateChannel", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
update_channel(arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = ivs("POST", "/UpdateChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), args)); aws_config=aws_config)


text:
using jlSPICE
using Test

@testset "jlSPICE.jl" begin
    # Write your tests here.
end


text:
@testset "c\"struct ...\"" begin
	@eval c``

	@eval c"""
	struct OpaqueStruct;
	"""
	@eval c"""
	struct OpaqueStruct;
	"""
	@test_throws ErrorException sizeof(c"struct OpaqueStruct")
	@test sizeof(c"struct OpaqueStruct *") == sizeof(Ptr{Cvoid})

	@eval c"""
	struct EmptyStruct {
	} __attribute__((packed));
	"""
	@test sizeof(c"struct EmptyStruct") == 0

	@eval c"""
	struct CcharStruct {
		char c;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct CcharStruct") == sizeof(Cchar)
	@test :c in fieldnames(c"struct CcharStruct")

	@eval c"""
	struct CintStruct {
		int i;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct CintStruct") == sizeof(Cint)
	@test :i in fieldnames(c"struct CintStruct")

	@eval c"""
	struct CstructStruct {
		struct CintStruct i;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct CstructStruct") == sizeof(c"struct CintStruct")
	@test :i in fieldnames(c"struct CstructStruct")

	@eval c"""
	struct CintCcharStruct {
		struct CintStruct i;
		struct CcharStruct c;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct CintCcharStruct") == sizeof(c"struct CintStruct")+sizeof(c"struct CcharStruct")
	@test :i in fieldnames(c"struct CintCcharStruct")
	@test :c in fieldnames(c"struct CintCcharStruct")

	@eval c"""
	struct CcharCintPackedStruct {
		struct CcharStruct c;
		struct CintStruct i;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct CcharCintPackedStruct") == sizeof(c"struct CintStruct")+sizeof(c"struct CcharStruct")
	@test :c in fieldnames(c"struct CcharCintPackedStruct")
	@test :i in fieldnames(c"struct CcharCintPackedStruct")

	@eval c"""
	struct CcharCintAlignedStruct {
		struct CcharStruct c;
		_Alignas(sizeof(struct CintStruct)*2) struct CintStruct i;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct CcharCintAlignedStruct") == max(sizeof(c"struct CintStruct")*4, max(sizeof(c"struct CintStruct"), sizeof(Cint))+max(sizeof(c"struct CcharStruct"), sizeof(Cint)))
	@test :c in fieldnames(c"struct CcharCintAlignedStruct")
	@test :i in fieldnames(c"struct CcharCintAlignedStruct")

	@eval c"""
	struct CunionStruct {
		union {
			unsigned char c[4];
			unsigned short s[2];
			unsigned int i;
		} __attribute__((packed));
	} __attribute__((packed));
	"""
	@test :c in fieldnames(c"struct CunionStruct")
	@test :s in fieldnames(c"struct CunionStruct")
	@test :i in fieldnames(c"struct CunionStruct")

	cus = c"struct CunionStruct"()
	@test cus.i == 0
	@test cus.c[1] == 0
	@test cus.c[4] == 0
	@test cus.s[1] == 0
	@test cus.s[2] == 0

	cus = c"struct CunionStruct"(cus, i = 0xff0000ff)
	@test cus.i == 0xff0000ff
	@test cus.c[1] == 0xff
	@test cus.c[4] == 0xff
	@test cus.s[1] != 0
	@test cus.s[2] != 0

	@eval c"""
	struct PtrToStructStruct;
	"""
	@eval c"""
	struct PtrToPtrToStructStruct {
		struct PtrToStructStruct *p;
	} __attribute__((packed));
	"""
	@eval c"""
	struct PtrToStructStruct {
		struct PtrToStructStruct *p;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct PtrToPtrToStructStruct") == sizeof(Ptr)
	@test :p in fieldnames(c"struct PtrToPtrToStructStruct")
	@test sizeof(c"struct PtrToStructStruct") == sizeof(Ptr)
	@test :p in fieldnames(c"struct PtrToStructStruct")

	@eval c"""
	struct Cint32Bitfield {
		int i:32;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cint32Bitfield") == sizeof(Cint)
	@test :i in fieldnames(c"struct Cint32Bitfield")

	@eval c"""
	struct Cuint32Bitfield {
		unsigned int u:32;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cuint32Bitfield") == sizeof(Cuint)
	@test :u in fieldnames(c"struct Cuint32Bitfield")

	@eval c"""
	struct Cint2Bitfield {
		int i:2;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cint2Bitfield") == sizeof(Cchar)
	@test :i in fieldnames(c"struct Cint2Bitfield")

	bf = c"struct Cint2Bitfield"()
	@test bf.i == 0
	for i in -3:3
		bf = c"struct Cint2Bitfield"(i = i)
		@test bf.i == (-2:1)[mod(i+2, 4)+1]
	end

	@eval c"""
	struct Cuint2Bitfield {
		unsigned int u:2;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cuint2Bitfield") == sizeof(Cchar)
	@test :u in fieldnames(c"struct Cuint2Bitfield")

	@eval c"""
	struct Cuint32Cint32Bitfields {
		unsigned int u:32;
		int i:32;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cuint32Cint32Bitfields") == sizeof(Cuint)+sizeof(Cint)
	@test :u in fieldnames(c"struct Cuint32Cint32Bitfields")
	@test :i in fieldnames(c"struct Cuint32Cint32Bitfields")

	@eval c"""
	struct Cuint16Cint16Bitfields {
		unsigned int u:16;
		int i:16;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cuint16Cint16Bitfields") == sizeof(Cuint)
	@test :u in fieldnames(c"struct Cuint16Cint16Bitfields")
	@test :i in fieldnames(c"struct Cuint16Cint16Bitfields")

	bf = c"struct Cuint16Cint16Bitfields"()
	@test bf.u == 0
	@test bf.i == 0

	bf = c"struct Cuint16Cint16Bitfields"(bf, i = Cint(-1))
	@test bf.u == 0
	@test bf.i == -1

	bf = c"struct Cuint16Cint16Bitfields"(bf, u = Cuint(0x0001))
	@test bf.u == 0x0001
	@test bf.i == -1

	@eval c"""
	struct Cint16Cint16Bitfields {
		int i1:16;
		int i2:16;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cint16Cint16Bitfields") == sizeof(Cuint)

	@eval c"""
	struct Cint16CintBitfields {
		unsigned int u:16;
		int i;
	} __attribute__((packed));
	"""
	@test sizeof(c"struct Cint16CintBitfields") == 2+sizeof(Cint)
	@test :u in fieldnames(c"struct Cint16CintBitfields")
	@test :i in fieldnames(c"struct Cint16CintBitfields")

	# https://github.com/analytech-solutions/CBinding.jl/issues/6
	@eval c"""
	struct AnonymousAggregateField {
		union {
			int x;
			char y;
		} u;
	};
	"""
	aaf = c"struct AnonymousAggregateField"(u = (x = -1,))
	@test aaf.u.x == -1
	aaf = c"struct AnonymousAggregateField"(aaf, u = (x = 0,))
	@test aaf.u.x == 0

	@eval c"""
	struct NestedAnonymousAggregateField {
		union {
			union {
				int z;
			} x;
			char y;
		} u;
	};
	"""
	aaf = c"struct NestedAnonymousAggregateField"(u = (x = (z = -1,),))
	@test aaf.u.x.z == -1
	aaf = c"struct NestedAnonymousAggregateField"(aaf, u = (x = (z = 0,),))
	@test aaf.u.x.z == 0

	# https://github.com/analytech-solutions/CBinding.jl/issues/8
	@eval c"""
	struct BrokenLayout {
		struct {
			int y;
			int z;
		} x;
	};
	"""
	bl = c"struct BrokenLayout"()
	@test bl.x.y == 0 && bl.x.z == 0
	bl = c"struct BrokenLayout"(bl, x = (z = 123,))
	@test bl.x.y == 0 && bl.x.z == 123

	@eval c"""
	struct ArrayOfStruct {
		struct BrokenLayout x[2];
	};
	"""
	@test sizeof(c"struct ArrayOfStruct") == 2*2*sizeof(Cint)
	aos = c"struct ArrayOfStruct"()
	@test typeof(aos) <: c"struct ArrayOfStruct"
	@test typeof(aos.x) <: Carray
	@test typeof(aos.x[2]) <: c"struct BrokenLayout"
	@test typeof(aos.x[2].x) <: Cstruct
	@test typeof(aos.x[2].x.y) <: Cint

	@eval c"""
	enum SmallEnum {
		SMALL_1 = 1,
		SMALL_2 = 2,
		SMALL_3 = 4,
	};
	"""
	@eval c"""
	struct SmallEnumBitfield {
		enum SmallEnum a:8;
		enum SmallEnum b:8;
	};
	"""
	@test sizeof(c"struct SmallEnumBitfield") == sizeof(c"enum SmallEnum")
	x = c"struct SmallEnumBitfield"()
	@test getfield(x, :mem) == (0x00, 0x00, 0x00, 0x00)
	x = c"struct SmallEnumBitfield"(x, a = c"SMALL_1")
	@test getfield(x, :mem) == (convert(UInt8, c"SMALL_1"), 0x00, 0x00, 0x00)
	x = c"struct SmallEnumBitfield"(x, b = c"SMALL_3")
	@test getfield(x, :mem) == (convert(UInt8, c"SMALL_1"), convert(UInt8, c"SMALL_3"), 0x00, 0x00)
end


text:
"""
    AbstractFunction

A function of the form:

```math
y = f(x,u; t)
```
where ``x`` is a state vector of size `n` and ``u`` is a control input of size `m` for some
dynamical system.  The output ``y`` is of size `p`.
Here ``t`` represents any number of optional extra parameters to the function (such as time,
time step, etc.).

The dimensions `n`, `m`, and `p` can
be queried individual by calling `state_dim`, `control_dim`, and `ouput_dim` or
collectively via [`dims`](@ref).

# Evaluation
An `AbstractFunction` can be evaluated in-place by calling any of

    evaluate!(fun, y, z::AbstractKnotPoint)
    evaluate!(fun, y, x, u, p)
    evaluate!(fun, y, x, u)

The top method should be preferred by generic APIs.

The function can also be evaluated out-of-place (generally returning a
`StaticArrays.SVector`) by calling any of

    evaluate(fun, z::AbstractKnotPoint)
    evaluate(fun, x, u, p)
    evaluate(fun, x, u)

Alternatively, the user can dispatch on [`FunctionSignature`](@ref) by calling

    evaluate!(::FunctionSignature, fun, y, z::AbstractKnotPoint)

# Jacobians
The Jacobian of ``f(x,u)`` with respect to both ``x`` and ``u`` can be computed
by calling the following function:

    jacobian!(::FunctionSignature, ::DiffMethod, fun, J, y, z::AbstractKnotPoint)

Where `J` is the Jacobian of size `(p, n + m)` and `y` is provided as an extra storage
vector for evaluating the Jacobian. Some users may whish to store the function output
value in this vector if evaluating the Jacobian and output value together offers some
computational savings, but this behavior is left to the user to implement on their
particular function.

A user-defined Jacobian, called by passing the `UserDefined` [`DiffMethod`](@ref),
can be implemented for the function by defining any of the following methods for
the `AbstractFunction`:

    jacobian!(fun, J, y, z::AbstractKnotPoint)
    jacobian!(fun, J, y, x, u, p)
    jacobian!(fun, J, y, x, u)

# Functions of just the state or control
Alternatively, the function can also be limited to an input of just the state or control
by defining the [`FunctionInputs`](@ref) trait. See trait documentation for more
information.

# Convenience functions
The methods `fill(fun, v)`, `randn(fun)`, `rand(fun)`, and `zeros(fun)` are defined
that provide a tuple with a state and control vector initialized using the corresponding
function. The data type can be specified as the first argument for `randn`, `rand`, and
`zeros`.
"""
abstract type AbstractFunction end

"""
    FunctionSignature

Specifies which function signature to call when evaluating the function, must be
either `StaticReturn` or `InPlace`. The default signature for a function can be
queried via `default_signature(fun)`.
"""
abstract type FunctionSignature end
struct InPlace <: FunctionSignature end
struct StaticReturn <: FunctionSignature end

"""
    DiffMethod

Represents the method used to evaluate the Jacobian of the function. Allows the
user to implement multiple methods and switch between them as needed. The following
methods are provided:

* `ForwardAD`: forward-mode automatic differentiation using
   [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
* `FiniteDifference`: finite difference approximation using
   [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl)
* `UserDefined`: analytical function defined by the user.
* `ImplicitFunctionTheorem`: Only available for a [`DiscretizedDynamics`](@ref) model
   with an an implicit integrator. Parameterized by the `DiffMethod` of the underlying
   continuous dynamics model.

The `UserDefined` method can be provided by defining any of the following methods
on your [`AbstractFunction`](@ref):

    jacobian!(fun, J, y, z::AbstractKnotPoint)
    jacobian!(fun, J, y, x, u, p)
    jacobian!(fun, J, y, x, u)

The other two modes must be explicitly defined on your function, since some functions
may wish to tailor the method to their particular function.  Efficient default
implementations can be defined automatically via the [`@autodiff`](@ref) macro.
"""
abstract type DiffMethod end
struct ForwardAD <: DiffMethod end
struct FiniteDifference <: DiffMethod end
struct UserDefined <: DiffMethod end
struct ImplicitFunctionTheorem{D<:DiffMethod} <: DiffMethod
    function ImplicitFunctionTheorem(diff::D) where D <: DiffMethod
        if D <: ImplicitFunctionTheorem
            throw(ArgumentError("DiffMethod for ImplicitFunctionTheorem cannot also be ImplicitFunctionTheorem."))
        end
        new{D}()
    end
end
default_diffmethod(::AbstractFunction) = UserDefined()
default_signature(::AbstractFunction) = StaticReturn()

"""
    dims(fun::AbstractFunction)

Get the tuple `(n,m,p)` where `n` is the dimension of the state vector, `m` is the
dimension of the control vector, and `p` is the output dimension of the function.
"""
dims(fun::AbstractFunction) = (state_dim(fun), control_dim(fun), output_dim(fun))
state_dim(fun::AbstractFunction) = throw(NotImplementedError("state_dim needs to be implemented for $(typeof(fun))."))
control_dim(fun::AbstractFunction) = throw(NotImplementedError("control_dim needs to be implemented for $(typeof(fun))."))
output_dim(fun::AbstractFunction) = throw(NotImplementedError("output_dim needs to be implemented for $(typeof(fun))"))
errstate_dim(fun::AbstractKnotPoint) = state_dim(fun)
jacobian_width(fun::AbstractFunction) = errstate_dim(fun) + control_dim(fun)

"""
    FunctionInputs

A trait of an [`AbstractFunction`](@ref) that specified the inputs to the function.
By default, the input is assumed to be both the state and the control vector, along with
any extra parameters. This trait allows the user to change this assumption, defining
functions on just the state or control vectors. The number of columns in the Jacobian
will change accordingly, and can be queried via [`input_dim`](@ref). Functions of
just the state or control are not allowed to take in any extra parameters as arguments.

This trait is queried via [`functioninputs(fun)`](@ref).

The following three options are provided:
* `StateControl`: a function of the form ``y = f(x,u)`` (default)
* `StateOnly`: a function of the form ``y = f(x)``
* `ControlOnly`: a function of the form ``y = f(u)``

When defining methods for these functions, you need to disambiguate methods of the form

    evaluate(fun::AbstractFunction, z::AbstractKnotPoint)

from

    evaluate(fun::MyStateOnlyFunction, x)

It's not enough to annotate `x` as an `AbstractVector` because
`AbstractKnotPoint <: AbstractVector`. For most cases, the `RobotDynamics.DataVector`
type should be sufficient to accomplish this. For `StateOnly` or `ControlOnly` functions,
your methods should look like:

    evaluate(fun::MyFunction, x::RobotDynamics.DataVector)
    evaluate!(fun::MyFunction, y, x::RobotDynamics.DataVector)
    jacobian!(fun::MyFunction, J, y, x::RobotDynamics.DataVector)

"""
abstract type FunctionInputs end
struct StateOnly <: FunctionInputs end
struct ControlOnly <: FunctionInputs end
struct StateControl <: FunctionInputs end

"""
    functioninputs(fun::AbstractFunction)

Get the [`FunctionInputs`](@ref) trait on the function, specifying whether the
function takes as input both state and control vectors (default), or just the
state or control vector independently.
"""
functioninputs(fun::AbstractFunction) = StateControl()

"""
    input_dim(fun::AbstractFunction)

Get the dimension of the inputs to the function. Will be equal to `n + m` for
a `StateVector` function, `n` for a `StateOnly` function, and `m` for a `ControlOnly`
function.
"""
input_dim(fun::AbstractFunction) = input_dim(functioninputs(fun), fun)
input_dim(::StateOnly, fun::AbstractFunction) = state_dim(fun)
input_dim(::ControlOnly, fun::AbstractFunction) = control_dim(fun)
input_dim(::StateControl, fun::AbstractFunction) = state_dim(fun) + control_dim(fun)

# Gets the input to the function given a KnotPoint as a single vector
getinput(::StateControl, z::AbstractKnotPoint) = getdata(z)
getinput(::StateOnly, z::AbstractKnotPoint) = state(z)
getinput(::ControlOnly, z::AbstractKnotPoint) = control(z)

# Gets the input to the function given a KnotPoint as separate arguments
getargs(::StateControl, z::AbstractKnotPoint) = state(z), control(z), getparams(z)
getargs(::StateOnly, z::AbstractKnotPoint) = (state(z),)
getargs(::ControlOnly, z::AbstractKnotPoint) = (control(z),)

# @inline getinputs(z::AbstractVector) = z
# @inline setinputs!(dest::AbstractVector{<:Real}, src::AbstractVector{<:Real}) = src
datatype(::Type{<:AbstractFunction}) = Float64

# Top-level command that can be overridden
# Should only be overridden if using hand-written Jacobian methods
evaluate!(fun::AbstractFunction, y, z::AbstractKnotPoint) = evaluate!(functioninputs(fun), fun, y, z)
evaluate(fun::AbstractFunction, z::AbstractKnotPoint) = evaluate(functioninputs(fun), fun, z)

evaluate!(inputtype::FunctionInputs, fun::AbstractFunction, y, z::AbstractKnotPoint) = evaluate!(fun, y, getargs(inputtype, z)...)
evaluate(inputtype::FunctionInputs, fun::AbstractFunction, z::AbstractKnotPoint) = evaluate(fun, getargs(inputtype, z)...)

# Strip the parameter
evaluate!(fun::AbstractFunction, y, x, u, p) = evaluate!(fun, y, x, u)
evaluate(fun::AbstractFunction, x, u, p) = evaluate(fun, x, u)

# Dispatch on function signature
evaluate!(::StaticReturn, fun::AbstractFunction, y, args...) = y .= evaluate(fun, args...)
evaluate!(::InPlace, fun::AbstractFunction, y, args...) = evaluate!(fun, y, args...)

# Minimal call that must be implemented
evaluate!(fun::AbstractFunction, y, x, u) =
    throw(NotImplementedError("User-defined in-place function not implemented for $(typeof(fun))"))
evaluate(fun::AbstractFunction, x, u) =
    throw(NotImplementedError("User-defined static return function not implemented for $(typeof(fun))"))

# Dispatch on FunctionSignature
evaluate!(::StaticReturn, fun::AbstractFunction, y, z::AbstractKnotPoint) =
    y .= evaluate(fun, z)
evaluate!(::InPlace, fun::AbstractFunction, y, z::AbstractKnotPoint) =
    evaluate!(fun, y, z)

# Jacobian
# Most generic method: this should the version called by external APIs
jacobian!(::FunctionSignature, ::UserDefined, fun::AbstractFunction, J, y, z) = jacobian!(fun, J, y, z)

# Highest-level function to implement for User-defined Jacobians
jacobian!(fun::AbstractFunction, J, y, z) = jacobian!(functioninputs(fun), fun, J, y, z)

# Dispatches on the input signature to call one of the following:
#  jacobian!(fun, J, y, x, u, p)  # StateControl
#  jacobian!(fun, J, y, x)        # StateOnly
#  jacobian!(fun, J, y, u)        # ControlOnly
jacobian!(inputtype::FunctionInputs, fun::AbstractFunction, J, y, z::AbstractKnotPoint) = jacobian!(fun, J, y, getargs(inputtype, z)...)

# Strip the parameter
jacobian!(fun::AbstractFunction, J, y, x, u, p) = jacobian!(fun, J, y, x, u)

# Minimal call to be implemented for StateControl inputs
jacobian!(fun::AbstractFunction, J, y, x, u) = throw(NotImplementedError("User-defined Jacobian not implemented for $(typeof(fun))"))

# Jacobian of Jacobian-vector product
∇jacobian!(::FunctionSignature, ::UserDefined, fun::AbstractFunction, H, b, y, z) =
    ∇jacobian!(fun, H, b, y, z)
∇jacobian!(fun::AbstractFunction, H, b, y, z::AbstractKnotPoint) =
    ∇jacobian!(fun, H, b, y, state(z), control(z), getparams(z))
∇jacobian!(fun::AbstractFunction, H, b, y, x, u, p) =
    ∇jacobian!(fun, H, b, y, state(z), control(z))
∇jacobian!(fun::AbstractFunction, H, b, y, x, u) =
    throw(UserDefined("User-defined Jacobian of Jacobian-vector product is undefined for $(typeof(fun))"))

# Some convenience methods
Base.rand(::Type{T}, model::AbstractFunction) where {T} = ((@SVector rand(T, state_dim(model))), (@SVector rand(T, control_dim(model))))
Base.randn(::Type{T}, model::AbstractFunction) where {T} = ((@SVector randn(T, state_dim(model))), (@SVector randn(T, control_dim(model))))
Base.zeros(::Type{T}, model::AbstractFunction) where {T} = ((@SVector zeros(T, state_dim(model))), (@SVector zeros(T, control_dim(model))))
Base.fill(model::AbstractFunction, v) = ((@SVector fill(v, state_dim(model))), (@SVector fill(v, control_dim(model))))

Base.rand(model::AbstractFunction) = Base.rand(Float64, model)
Base.randn(model::AbstractFunction) = Base.randn(Float64, model)
Base.zeros(model::AbstractFunction) = Base.zeros(Float64, model)


text:
import CommutativeRings: tosuper, tosub

@testset "scriptdigits" begin

    @test tosuper(0) == "⁰"
    @test tosuper(0, sign=true) == "⁺⁰"
    @test tosuper(1234) == "¹²³⁴"
    @test tosuper(1234, sign=true) == "⁺¹²³⁴"
    @test tosuper(-1234) == "⁻¹²³⁴"
    @test tosuper(-1234, sign=true) == "⁻¹²³⁴"
    @test tosub(0) == "₀"
    @test tosub(0, sign=true) == "₊₀"
    @test tosub(1234) == "₁₂₃₄"
    @test tosub(1234, sign=true) == "₊₁₂₃₄"
    @test tosub(-1234) == "₋₁₂₃₄"
    @test tosub(-1234, sign=true) == "₋₁₂₃₄"

end

@testset "degrees" begin
    @test deg(0) == -1
    @test deg(1) == 0
    @test deg(ZZ(0)) == -1
    @test deg(ZZ(10)) == 0
end

@testset "testrules $R" for R in (ZZ/25, GF(25))
    io = IOBuffer()
    @test sprint(CommutativeRings.testrules, R) == ""
end

@testset "basetypes($G)" for G in (Int, ZZ{Int}, ZZ/2, (ZZ/3)[:x], GF(2,2))
    @test length(basetypes(G)) == depth(G) + 1
end

@testset "homomorphisms" begin
    R = ZZ/7
    S = GF(7, 2)
    h = Hom{R,S}(S)
    @test h(R(2)) == S(2)
end

@testset "isfield $F" for (F, r) in ((ZZ{Int}, false), (ZZ/6, false), (ZZ/3, true), (QQ{Int}, true), (ZZ{Int}[:x], false))
    @test isfield(F) == r
end

@testset "promotion1" begin
    P = QQ{Int}[:x]
    x = monom(P)
    Q = P / (x^2 - 2)
    @test x + 1//7 isa P
    @test 1//7 + x isa P
    @test Q(x) + 1 // 7 isa Q
    @test 1 // 7 + Q(x) isa Q
end

@testset "promotion2" begin
    P = ZZ{Int}[:x]
    x = monom(P)
    Q = P / ( x^2 - 3)
    PP = Q[:y]
    y = monom(PP)
    @test x * y == Q(x) * y
    @test Q(x) * y isa PP
end

@testset "promotion3" begin
    P = ZZ{Int}[:x]
    x = monom(P)
    R = P[:y]
    y = monom(R)
    S = R[:z]
    z = monom(S)
    @test z + x isa S
    @test (x + y + z)^2 == x^2 + y^2 + z^2 + 2x*y + 2x*z + 2y*z

end

using CommutativeRings: FieldTrait, EuclidianDomainTrait, UniqueFactorizationDomainTrait, IntegralDomainTrait, CommutativeRingTrait

@testset "category traits $T" for (T,C) in [(ZZ, EuclidianDomainTrait),
                                            (ZZ/2, FieldTrait),
                                            (GF(4), FieldTrait),
                                            (QQ{BigInt}, FieldTrait),
                                            (ZZ/6, CommutativeRingTrait),
                                            (ZZ{Int}[:x], UniqueFactorizationDomainTrait),
                                            (ZZ{Int}[:x][:y], UniqueFactorizationDomainTrait),
                                            (ZZ{Int}[:x,:y], UniqueFactorizationDomainTrait),
                                            ((ZZ/3)[:x], EuclidianDomainTrait),
                                            ((ZZ/3)[:x,:y], UniqueFactorizationDomainTrait)
                                          ]
    @test category_trait(T) == C
end

const X = monom(ZZ{Int}[:x])
@testset "category traits $T/($p)" for (T, p, C) in [(ZZ/5, X^2 + 2, FieldTrait),
                                                     (ZZ/5, X^2 + 1, CommutativeRingTrait),
                                                     (ZZ{Int}, X^2 + 1, FieldTrait),
                                                     (ZZ{Int}, X^2 - 1, CommutativeRingTrait)
                                                  ]
    x = monom(T[:x])
    P = p(x)
    @test category_trait(T[:x]/P) == C
end

@testset "category traits Frac{$T}" for (T, C) in [(ZZ{Int}, FieldTrait),
                                                   (ZZ{Int}[:x], FieldTrait),
                                                   ((ZZ/6)[:x], CommutativeRingTrait),
                                                  ]
    @test category_trait(Frac{T}) == C
end


text:
module ParetoSets

export dominates, naive_pareto
export ParetoSet, pareto_ids, pareto_ys, dominated_ids, dominated_ys, all_ys, ispareto

"""
    dominates(y, y´)

Returns true if y dominates y´, false otherwise.
"""
dominates(y::Vector{Float64}, y´::Vector{Float64}) = all(y .≤ y´) && any(y .< y´)

"""
    naive_pareto{T}(ys::AbstractVector{Vector{T}})

Returns the pareto set of a collection of points ys as a vector of indices into ys.  Uses naive O(n^2) algorithm.
"""
function naive_pareto(ys::AbstractVector{Vector{Float64}})
    pareto_ids = Int[]
    for (i,y) in enumerate(ys)
        if !any(dominates(y´,y) for y´ in ys)
            push!(pareto_ids, i)
        end
    end
    return pareto_ids
end

mutable struct ParetoSet
    ys::Vector{Vector{Float64}}  #collection of all points
    pareto_ids::Vector{Int}  #vector of indices into full containing pareto set
end
ParetoSet() = ParetoSet(Vector{Float64}[], Int[])

function Base.push!(p::ParetoSet, y::Vector{Float64})
    push!(p.ys, y)
    yi = length(p.ys)
    ispareto(p, y) && push!(p.pareto_ids, yi)
    ids = find(y´->dominates(y,y´), pareto_ys(p))
    deleteat!(p.pareto_ids, ids)
end

pareto_ids(p::ParetoSet) = p.pareto_ids
dominated_ids(p::ParetoSet) = setdiff(1:length(p.ys), p.pareto_ids)
pareto_ys(p::ParetoSet) = view(p.ys, p.pareto_ids)
dominated_ys(p::ParetoSet) = view(p.ys, dominated_ids(p))
all_ys(p::ParetoSet) = p.ys

Base.getindex(p::ParetoSet, i) = p.ys[i]
function Base.append!(ps::ParetoSet, ys)
    for y in ys
        push!(ps, y)
    end
end

ispareto(p::ParetoSet, y::Vector{Float64}) = !any(dominates(y´, y) for y´ in pareto_ys(p))

end # module


text:
# module for landmark estimation by graph based slam

using LinearAlgebra

mutable struct MapEdge
  x
  z
  m
  omega
  xi

  # init
  # head_t: time when target landmark was observed for the first time
  # head_z: observed values when target landmark was observed for the first time
  function MapEdge(t, z, head_t, head_z, pose_list;
                   sensor_noise_rate=[0.14, 0.05])
    self = new()
    self.x = pose_list[string(t)]
    self.z = z

    # landmark's pose
    delta_obsrv = [
      self.z[1]*cos(self.x[3]+self.z[2]),
      self.z[1]*sin(self.x[3]+self.z[2])
    ]
    self.m = self.x[1:2] + delta_obsrv

    # precision matrix
    Q1 = diagm(0 => [(self.z[1]*sensor_noise_rate[1])^2,
                     sensor_noise_rate[2]^2])
    s1 = sin(self.x[3] + self.z[2])
    c1 = cos(self.x[3] + self.z[2])
    R = [-c1  self.z[1]*s1;
         -s1 -self.z[1]*c1]
    self.omega = R * Q1 * R'

    # coefficient vector
    self.xi = self.omega * self.m

    return self
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module KMeans
  using ..KNeighbors
  using ...Statistics

  function runKMeans(centroids,dataset,k)
    assignment = 1:k
    return KNeighbors.KNearestNeighbors(centroids,dataset,assignment)
  end

  function trainKmeans(training_set,k,seed)
    changed = true;
    #size is k by n
    centroids = initCentroids(training_set,k,seed)
    assignment = 1:k
    y=zeros(Int64,size(training_set)[1],1)
    index = 0;
    while(changed)
      index+=1
      y=KNeighbors.KNearestNeighbors(centroids,training_set,assignment)
      newCentroids = recalculateCentroids(centroids,training_set,y);
      if centroids==newCentroids
        changed=false
      end
      centroids = newCentroids
    end
    return centroids,y;
  end
  #
  function initCentroids(training_set,k,seed)
    indices = []
    srand(seed)
    indices=randperm(size(training_set)[1])[1:k]
    return training_set[indices, :]
  end

  function recalculateCentroids(centroids,training_set, y)
    newCentroids = zeros(Float64,size(centroids))
    for i=1:size(centroids)[1]
      newCentroids[i,:] = Statistics.meanColumn(training_set[find(x->x==i,y),:])
    end
    return newCentroids
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
# #infinitemps.jl
# struct InfiniteMPS{S<:IndexSpace} <: AbstractMPS{S}
#     sitetensors::Vector{AbstractTensorMap{S, 2, 1}}
#     schmidtvalues::Vector{Vector{Float64}}
#     canonical_form::Symbol
# end
#
# InfiniteMPS(sitetensors,schmidtvalues=nothing,canonical_form=:non_canonical) =
#       InfiniteMPS{spacetype(sitetensors[1])}(sitetensors, schmidtvalues, canonical_form)


text:
module TestWeight
using Test,CorrDCA,DelimitedFiles

Z = [1 2 3 3 3;
    1 2 4 3 3;
    3 2 1 3 3;]

W,Meff = compute_weights(Z,0.9)

@test W == [0.25,0.25,0.25, 0.125,0.125]
@test Meff == 4.0
W,Meff = compute_weights(Z,:auto)
@test W == [0.25,0.25,0.25,0.125,0.125]
@test Meff == 4.0

datadir = joinpath(dirname(@__FILE__),"data")

fileal = joinpath(datadir,"pf14short.fasta.gz")
fileW = joinpath(datadir,"W.dat")
Z = read_fasta_alignment(fileal,0.8)
Wtrue = readdlm(fileW) |> vec
W,Meff=compute_weights(Z,:auto)
Mefftrue = sum(Wtrue)
@test Meff == Mefftrue
@test W .* Meff == Wtrue

printstyled("All TestWeight passed!\n",color=:light_green,bold=true)
end # end module

max_stars_repo_name: pagnani/CorrDCA

max_issues_repo_name: pagnani/CorrDCA

max_forks_repo_name: pagnani/CorrDCA


text:
using TypeVars
using Test

struct Zm{M, T<:Integer}
    val::T
end

import Base: +
+(a::Z, b::Z) where Z<:Zm = Z(mod(a.val + b.val, typevar(Z)))

@testset "TypeVars" begin

    m = big"2"^512 -1
    M = Symbol(hash(m))
    Z = Zm{M,BigInt}
    settypevar(Z, m)

    a = Z(m - 5)
    b = Z(m - 10)
    c = a + b
    @test c.val == m - 15
end


text:
__precompile__()
module COSMO

using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg, DataStructures, Requires, Printf, IterTools, GenericLinearAlgebra, MutableArithmetics
using Reexport
using COSMOAccelerators

using Reexport
@reexport using COSMOAccelerators

export assemble!, warm_start!, empty_model!, update!

const DefaultFloat = Float64
const DefaultInt   = LinearAlgebra.BlasInt
const MA = MutableArithmetics

include("./kktsolver.jl")
# optional dependencies
function __init__()
    @require Pardiso="46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" include("kktsolver_pardiso.jl")

    @require IterativeSolvers="42fd0dbc-a981-5370-80f2-aaf504508153" begin
        @require LinearMaps="7a12625a-238d-50fd-b39a-03d52299707e" include("./kktsolver_indirect.jl")
    end
end

function version()
    v"0.8.1"
end

include("./algebra.jl")
include("./projections.jl")
#include("./acceleration.jl")
include("./trees.jl")
include("./clique_graph.jl")
include("./clique_merging.jl")
include("./settings.jl")
include("./types.jl")
include("./accelerator_interface.jl")

include("./constraint.jl")
include("./parameters.jl")
include("./residuals.jl")
include("./scaling.jl")
include("./infeasibility.jl")
include("./transformations.jl")
include("./chordal_decomposition.jl")
include("./printing.jl")
include("./setup.jl")
include("./solver.jl")
include("./interface.jl")
include("./MOI_wrapper.jl")

# include("precompile.jl")
# _precompile_()

end #end module


text:
"""
    ut1tai(ut11, ut12, dta)

Time scale transformation:  Universal Time, UT1, to International
Atomic Time, TAI.

### Given ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date
- `dta`: UT1-TAI in seconds

### Returned ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date

### Notes ###

1. ut11+ut12 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where ut11 is the Julian
   Day Number and ut12 is the fraction of a day.  The returned
   tai1,tai2 follow suit.

2. The argument dta, i.e. UT1-TAI, is an observed quantity, and is
   available from IERS tabulations.

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
ut1tai

"""
    ut1tt(ut11, ut12, dt)

Time scale transformation:  Universal Time, UT1, to Terrestrial
Time, TT.

### Given ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date
- `dt`: TT-UT1 in seconds

### Returned ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Notes ###

1. ut11+ut12 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where ut11 is the Julian
   Day Number and ut12 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

2. The argument dt is classical Delta T.

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
ut1tt

"""
    ut1utc(ut11, ut12, dut1)

Time scale transformation:  Universal Time, UT1, to Coordinated
Universal Time, UTC.

### Given ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date (Note 1)
- `dut1`: Delta UT1: UT1-UTC in seconds (Note 2)

### Returned ###

- `utc1`, `utc2`: UTC as a 2-part quasi Julian Date (Notes 3,4)

### Notes ###

1. ut11+ut12 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where ut11 is the Julian
   Day Number and ut12 is the fraction of a day.  The returned utc1
   and utc2 form an analogous pair, except that a special convention
   is used, to deal with the problem of leap seconds - see Note 3.

2. Delta UT1 can be obtained from tabulations provided by the
   International Earth Rotation and Reference Systems Service.  The
   value changes abruptly by 1s at a leap second;  however, close to
   a leap second the algorithm used here is tolerant of the "wrong"
   choice of value being made.

3. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the returned quasi JD day UTC1+UTC2 represents
   UTC days whether the length is 86399, 86400 or 86401 SI seconds.

4. The function eraD2dtf can be used to transform the UTC quasi-JD
   into calendar date and clock time, including UTC leap second
   handling.

5. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See eraDat for further details.

### Called ###

- `eraJd2cal`: JD to Gregorian calendar
- `eraDat`: delta(AT) = TAI-UTC
- `eraCal2jd`: Gregorian calendar to JD

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
ut1utc

"""
    utcut1(utc1, utc2, dut1)

Time scale transformation:  Coordinated Universal Time, UTC, to
Universal Time, UT1.

### Given ###

- `utc1`, `utc2`: UTC as a 2-part quasi Julian Date (Notes 1-4)
- `dut1`: Delta UT1 = UT1-UTC in seconds (Note 5)

### Returned ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date (Note 6)

### Notes ###

1. utc1+utc2 is quasi Julian Date (see Note 2), apportioned in any
   convenient way between the two arguments, for example where utc1
   is the Julian Day Number and utc2 is the fraction of a day.

2. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the JD day represents UTC days whether the
   length is 86399, 86400 or 86401 SI seconds.

3. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See eraDat for further details.

4. The function eraDtf2d converts from calendar date and time of
   day into 2-part Julian Date, and in the case of UTC implements
   the leap-second-ambiguity convention described above.

5. Delta UT1 can be obtained from tabulations provided by the
   International Earth Rotation and Reference Systems Service.
   It is the caller's responsibility to supply a dut1 argument
   containing the UT1-UTC value that matches the given UTC.

6. The returned ut11,ut12 are such that their sum is the UT1 Julian
   Date.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

### Called ###

- `eraJd2cal`: JD to Gregorian calendar
- `eraDat`: delta(AT) = TAI-UTC
- `eraUtctai`: UTC to TAI
- `eraTaiut1`: TAI to UT1

"""
utcut1

for name in ("ut1tai",
             "ut1tt",
             "ut1utc",
             "utcut1")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(a, b, c)
            r1 = Ref(0.0)
            r2 = Ref(0.0)
            i = ccall(($fc, liberfa), Cint,
                      (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                      a, b, c, r1, r2)
            @assert i == 0
            r1[], r2[]
        end
    end
end

"""
    utctai(utc1, utc2)

Time scale transformation:  Coordinated Universal Time, UTC, to
International Atomic Time, TAI.

### Given ###

- `utc1`, `utc2`: UTC as a 2-part quasi Julian Date (Notes 1-4)

### Returned ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date (Note 5)

### Notes ###

1. utc1+utc2 is quasi Julian Date (see Note 2), apportioned in any
   convenient way between the two arguments, for example where utc1
   is the Julian Day Number and utc2 is the fraction of a day.

2. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the JD day represents UTC days whether the
   length is 86399, 86400 or 86401 SI seconds.  In the 1960-1972 era
   there were smaller jumps (in either direction) each time the
   linear UTC(TAI) expression was changed, and these "mini-leaps"
   are also included in the ERFA convention.

3. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See eraDat for further details.

4. The function eraDtf2d converts from calendar date and time of day
   into 2-part Julian Date, and in the case of UTC implements the
   leap-second-ambiguity convention described above.

5. The returned TAI1,TAI2 are such that their sum is the TAI Julian
   Date.

### Called ###

- `eraJd2cal`: JD to Gregorian calendar
- `eraDat`: delta(AT) = TAI-UTC
- `eraCal2jd`: Gregorian calendar to JD

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
function utctai(a, b)
    r1 = Ref(0.0)
    r2 = Ref(0.0)
    i = ccall((:eraUtctai, liberfa), Cint,
                (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                a, b, r1, r2)
    if i == 1
        @warn "dubious year (Note 3)"
    elseif i == -1
        throw(ERFAException("unacceptable date"))
    end
    r1[], r2[]
end


text:
@defcomp impactcooling begin
    regions = Index()

    cooling = Variable(index=[time,regions])

    cebm = Parameter(index=[regions])
    gdp90 = Parameter(index=[regions])

    population = Parameter(index=[time,regions])
    pop90 = Parameter(index=[regions])

    income = Parameter(index=[time,regions])
    ceel = Parameter(default = 0.8)

    temp = Parameter(index=[time,regions])
    cenl = Parameter(default = 1.5)

    cumaeei = Parameter(index=[time,regions])

    function run_timestep(p, v, d, t)

        if !is_first(t)
            for r in d.regions
                ypc = p.income[t, r] / p.population[t, r] * 1000.0
                ypc90 = p.gdp90[r] / p.pop90[r] * 1000.0

                v.cooling[t, r] = p.cebm[r] * p.cumaeei[t, r] * p.gdp90[r] * (p.temp[t, r] / 1.0)^p.cenl * (ypc / ypc90)^p.ceel * p.population[t, r] / p.pop90[r]
            end
        end
    end
end


text:
@testset "Magnitude of completeness" begin
    b = 1
    M0 = 0.5
    mags = rand(GRclassic(b,M0), 100000)
    filter!(x->rand() < cdf(Normal(0.8,0.2),x), mags)
    ΔM = 0.1

    x, y = histmag(mags,ΔM)
    @test round(maxcurv(mags, ΔM), digits=2) == 0.95
    @test maxcurv(mags, ΔM) == maxcurv(x,y)
    @test_throws ErrorException("The two input arrays must have the same number of elements") maxcurv([0.1],[1.2,1.6])

    @test round(mbass(mags, ΔM), digits=2) == 1.25
    @test mbass(mags, ΔM) == mbass(x,y)
    @test_throws ErrorException("The two input arrays must have the same number of elements") mbass([0.1],[1.2,1.6])
    @test isnan( mbass(collect(0.1:0.1:1.1),fill(11,11)) )
    @test isnan( mbass([0.1,0.2],[1.2,1.6]) )

    @test round(mbs(mags,ΔM,winsize=5, minevents=50), digits=2) == 1.55
    @test isnan( mbs([0.1,0.2],0.1) )

end