text:
using DataStructures: PriorityQueue, enqueue!, dequeue!, dequeue_pair!,
IntDisjointSets, in_same_set, union!
using CUDA

export learn_chow_liu_tree, learn_chow_liu_trees

"Learn Chow Liu Tree(s). It will run on CPU/GPU based on where `train_x` is.

Arguments:
- `train_x`: training data. If want gpu, move to gpu before calling this `CuArray(train_x)`.

Keyword arguments:
- `num_trees=1`: number of trees you want to learned.
- `dropout_prob=0.0`: drop edges with probability `dropout_prob` when learning maximum spanning tree.
- `weights=nothing`: weights of samples. Weights are all 1 if `nothing`.
- `pseudocount=0.0`: add a total of pseudo count spread out overall all categories.
- `Float=Float64`: precision. `Float32` is faster if `train_x` a large.
"
function learn_chow_liu_trees(train_x;
        num_trees=1, dropout_prob=0.0, weights=nothing,
        pseudocount=1.0, Float=Float64)
    MI = pairwise_MI(train_x; weights, pseudocount, Float)
    if MI isa CuArray
        MI_cpu = Array(MI) # TODO: GPU MST
        CUDA.unsafe_free!(MI)
    else
        MI_cpu = MI
    end

    topk_MST(-MI_cpu; num_trees, dropout_prob)
end

function learn_chow_liu_tree(train_x;
    dropout_prob=0.0, weights=nothing,
    pseudocount=1.0, Float=Float64)
    learn_chow_liu_trees(train_x;
        num_trees=1, dropout_prob, weights, pseudocount, Float)[1]
end

"Top k minimum spanning trees for a complete graph with `weights` as weights
http://www.nda.ac.jp/~yamada/paper/enum-mst.pdf
"
function topk_MST(weights::Matrix;
        num_trees::Integer=1, dropout_prob::Float64=0.0)

    # Priority queue that maintain candidate MSTs
    T = Vector{Tuple{Int, Int}}
    candidates = PriorityQueue{Tuple{T, T, T}, Float32}()

    included_edges::T = T()
    excluded_edges::T = T()
    reuse = similar(weights)
    topk_msts::Vector{T} = Vector{T}()

    # Initialize `candidate` with the global MST
    mst_edges, total_weight = MST(weights, included_edges, excluded_edges;
                                  reuse, dropout_prob=0.0)
    enqueue!(candidates, (mst_edges, included_edges, excluded_edges), total_weight)

    if Threads.nthreads() == 1

        # Sequential code
        for idx = 1 : num_trees
            if isempty(candidates)
                break
            end

            (mst_edges, included_edges, excluded_edges), total_weight = dequeue_pair!(candidates)

            # Record the current ST into `topk_msts`
            push!(topk_msts, mst_edges)

            if idx == num_trees
                break
            end

            edge_added = false
            for edge_idx = 1 : length(mst_edges)
                if mst_edges[edge_idx] in included_edges
                    continue
                end

                if edge_added
                    push!(included_edges, pop!(excluded_edges))
                end
                push!(excluded_edges, mst_edges[edge_idx])
                edge_added = true

                candidate_mst, total_weight = MST(weights, included_edges, excluded_edges;
                                                reuse, dropout_prob)
                if candidate_mst !== nothing
                    # A shallow copy of the vectors `included_edges` and `excluded_edges` is sufficient
                    enqueue!(candidates, (candidate_mst, copy(included_edges), copy(excluded_edges)), total_weight)
                end
            end
        end

    else

        # Parallel code
        reuse = map(1:Threads.nthreads()) do idx
            similar(weights)
        end

        weights = map(1:Threads.nthreads()) do idx
            deepcopy(weights)
        end

        l = ReentrantLock()

        for idx = 1 : num_trees
            if isempty(candidates)
                break
            end

            (mst_edges, included_edges, excluded_edges), total_weight = dequeue_pair!(candidates)

            # Record the current ST into `topk_msts`
            push!(topk_msts, mst_edges)

            if idx == num_trees
                break
            end

            Threads.@threads for edge_idx = 1 : length(mst_edges)
                curr_included_edges = copy(included_edges)
                curr_excluded_edges = copy(excluded_edges)
                for edge in mst_edges[1:edge_idx-1]
                    if !(edge in included_edges)
                        push!(curr_included_edges, edge)
                    end
                end
                if !(mst_edges[edge_idx] in excluded_edges)
                    push!(curr_excluded_edges, mst_edges[edge_idx])
                end

                id = Threads.threadid()
                candidate_mst, total_weight = MST(weights[id],
                                                curr_included_edges,
                                                curr_excluded_edges;
                                                reuse = reuse[id],
                                                dropout_prob)

                lock(l)
                if candidate_mst !== nothing
                    # A shallow copy of the vectors `included_edges` and `excluded_edges` is sufficient
                    enqueue!(candidates, (candidate_mst, copy(included_edges), copy(excluded_edges)), total_weight)
                end
                unlock(l)
            end
        end

    end

    # Post-process the top-K Spanning Trees
    topk_msts
end

"Compute the Minimum Spanning Tree (MST) of graph g with weights `weights`, with
 constraints such that `included_edges` should be included while `excluded_edges`
 should be excluded."
function MST(weights::Matrix,
             included_edges::Vector{Tuple{Int, Int}},
             excluded_edges::Vector{Tuple{Int, Int}};
             reuse::Matrix, dropout_prob=0.0)
    T = eltype(weights)
    @inbounds @views reuse[:, :] .= weights[:, :]

    # Dropout
    if dropout_prob > 1e-8
        dropped_mask = rand(Bernoulli(dropout_prob), size(reuse, 1), size(reuse, 2))
        @inbounds @views reuse[dropped_mask] .= Inf
    end

    # Add constraints
    map(included_edges) do edge
        reuse[edge[1], edge[2]] = -Inf
        reuse[edge[2], edge[1]] = -Inf
        nothing # Return nothing to save some effort
    end
    map(excluded_edges) do edge
        reuse[edge[1], edge[2]] = Inf
        reuse[edge[2], edge[1]] = Inf
        nothing # Return nothing to save some effort
    end

    mst_edges = kruskal_mst_complete(reuse)

    # Senity check
    valid_tree::Bool = true

    edges = Set(mst_edges)
    map(included_edges) do edge
        if !(edge in edges)
            valid_tree = false
        end
        nothing
    end
    map(excluded_edges) do edge
        if (edge in edges)
            valid_tree = false
        end
        nothing
    end

    if valid_tree
        # Compute the tree weight
        total_weight::T = 0.0
        map(mst_edges) do edge
            total_weight += weights[edge[1], edge[2]]
            nothing
        end
        mst_edges, total_weight
    else
        nothing, nothing
    end
end

function kruskal_mst_complete(distmx::Matrix; minimize=true)
    T = eltype(distmx)

    nv = size(distmx, 1)
    ne = Int(nv * (nv - 1) / 2)
    connected_vs = IntDisjointSets(nv)

    mst = Vector{Tuple{Int64, Int64}}()
    sizehint!(mst, nv - 1)

    weights = Vector{T}()
    sizehint!(weights, ne)
    edge_list = []
    for i in 1 : nv
        for j in i + 1 : nv
            push!(edge_list, (i, j))
        end
    end

    for e in edge_list
        push!(weights, distmx[e[1], e[2]])
    end

    for e in edge_list[sortperm(weights; rev=!minimize)]
        if !in_same_set(connected_vs, e[1], e[2])
            union!(connected_vs, e[1], e[2])
            push!(mst, e)
            (length(mst) >= nv - 1) && break
        end
    end
    return mst
end


text:
"""
Call `f` with each AMP spell in `records`. Useful for accumulators.
"""
function traverse_AMP_spells(f, records)
    for data in values(records)
        foreach(f, data.AMP_spells)
    end
end

ext: jl


text:
"""
    ConfigTypes

Module containing ClimateMachine configuration types.
"""
module ConfigTypes

export ClimateMachineConfigType,
    AtmosConfigType,
    AtmosLESConfigType,
    AtmosGCMConfigType,
    OceanConfigType,
    OceanBoxGCMConfigType,
    OceanSplitExplicitConfigType,
    SingleStackConfigType,
    MultiColumnLandConfigType

abstract type ClimateMachineConfigType end
abstract type AtmosConfigType <: ClimateMachineConfigType end
struct AtmosLESConfigType <: AtmosConfigType end
struct AtmosGCMConfigType <: AtmosConfigType end
abstract type OceanConfigType <: ClimateMachineConfigType end
struct OceanBoxGCMConfigType <: OceanConfigType end
struct OceanSplitExplicitConfigType <: OceanConfigType end
struct SingleStackConfigType <: ClimateMachineConfigType end
struct MultiColumnLandConfigType <: ClimateMachineConfigType end
end


text:
# PROJECT_CONSENSUS    Bring a set of points into consensus.
#
#   project_consensus(v1,v2,...,vn) returns the elementwise
#   average of v1, ..., vn, i.e., the projection of the vi
#   onto the consensus set.

function project_consensus(varargin)
    N = length(varargin)
    x = zeros(size(varargin[1]))

    for i = 1:N
        x = x + varargin[i]
    end
    x = x ./ N

    return x
end


text:
include("zygote_patch.jl")

import YaoExtensions, Random

c = YaoExtensions.variational_circuit(5)
dispatch!(c, :random)

function loss(reg::AbstractRegister, circuit::AbstractBlock{N}) where N
    #copy(reg) |> circuit
    reg = apply!(copy(reg), circuit)
    st = state(reg)
    sum(real(st.*st))
end

reg0 = zero_state(5)
paramsδ = gradient(c->loss(reg0, c), c)[1]
regδ = gradient(reg->loss(reg, c), reg0)[1]


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

module FlockingExample

export run_flocking_example!

using Agents
using LinearAlgebra
using InteractiveDynamics
using CairoMakie
using Makie.GeometryBasics # without this UndefVarError: Point2f0

mutable struct Bird <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    speed::Float64
    cohere_factor::Float64
    separation::Float64
    separate_factor::Float64
    match_factor::Float64
    visual_distance::Float64
end

function random_angle(rng)
    rand_deg = rand(rng,0:360)
    rad = deg2rad(rand_deg)
    return (cos(rad), sin(rad))
end

function initialize_model(;
    n_birds = 100,
    speed = 1.0,
    cohere_factor = 0.25,
    separation = 4.0,
    separate_factor = 0.25,
    match_factor = 0.01,
    visual_distance = 5.0,
    extent = (100, 100),
    spacing = visual_distance / 1.5,
)
    space2d = ContinuousSpace(extent, spacing)
    model = ABM(Bird, space2d, scheduler = Schedulers.randomly)
    for _ in 1:n_birds
        # randomly assign the heading for each bird
        # the * 2 .-1 ensures that we are covering all 360deg
        # otherwise rand will only generate +ve numbers
        # which effectively locks us within 90*
        # vel = Tuple(rand(model.rng, 2) * 2 .- 1)
        vel = random_angle(model.rng)
        add_agent!(
            model,
            vel,
            speed,
            cohere_factor,
            separation,
            separate_factor,
            match_factor,
            visual_distance,
        )
    end
    return model
end

"""
    heading_to_new_pos(current_pos, new_pos)
The heading that agent should turn to go towards the new_pos
Return `heading`
"""
heading_to_new_pos(current_pos, new_pos) = new_pos .- current_pos

function agent_step!(bird, model)
    # Obtain the ids of neighbors within the bird's visual distance
    neighbor_ids = nearby_ids(bird, model, bird.visual_distance)
    N = 0
    match = separate = cohere = (0.0, 0.0)
    # Calculate behaviour properties based on neighbors
    for id in neighbor_ids
        N += 1
        neighbor = model[id]
        # the heading that the bird has to turn, if it were to go towards the neighbor bird
        heading = heading_to_new_pos(bird.pos, neighbor.pos)
        # `cohere` computes the average position of neighboring birds
        cohere = cohere .+ heading
        if edistance(bird.pos, neighbor.pos, model) < bird.separation
            # `separate` repels the bird away from neighboring birds
            # ' - heading' indicates head in the opposite direction of the neighboring bird
            separate = separate .- heading
        end
        # `match` computes the average trajectory of neighboring birds
        match = match .+ model[id].vel
    end
    N = max(N, 1)
    # Normalise results based on model input and neighbor count
    cohere = cohere ./ N .* bird.cohere_factor
    separate = separate ./ N .* bird.separate_factor
    match = match ./ N .* bird.match_factor
    # Compute velocity (heading) based on rules defined above
    bird.vel = (bird.vel .+ cohere .+ separate .+ match) ./ 2
    bird.vel = bird.vel ./ norm(bird.vel)
    # Move bird according to new velocity and speed
    move_agent!(bird, model, bird.speed)
end

const bird_polygon = Polygon(Point2f0[(-0.5, -0.5), (1, 0), (-0.5, 0.5)])

function bird_marker(b::Bird)
    φ = atan(b.vel[2], b.vel[1]) #+ π/2 + π
    scale(rotate2D(bird_polygon, φ), 2)
end

function run_flocking_example!()
    model = initialize_model()
    # And let's also do a nice little video for it:
    abm_video(
        "flocking.mp4", model, agent_step!;
        am = bird_marker,
        framerate = 20, frames = 200,
        title = "Flocking"
    )
end
#run_flocking_example!()
end


text:
module RulesSubtractionMarginalsTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_marginalrules

@testset "marginalrules:typeof(-)" begin

    @testset ":in1_in2 (m_out::UnivariateNormalDistributionsFamily, m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = NormalMeanVariance(3.0, 4.0), m_in1 = NormalMeanVariance(2.0, 2.0), m_in2 = PointMass(2.0)),
                output = (in1 = NormalWeightedMeanPrecision(9/4, 3/4), in2 = PointMass(2.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(3, 4), m_in1 = NormalMeanPrecision(1, 4), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(20.0, 8.0), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(3.0, 4.0), m_in1 = NormalWeightedMeanPrecision(1.0, 4.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(8.0, 8.0), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(2.0, 4.0), m_in1 = NormalMeanVariance(1.0, 2.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(25/2, 9/2), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalMeanVariance(2.0, 4.0), m_in1 = NormalMeanPrecision(1.0, 2.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(11/4, 9/4), in2 = PointMass(1.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(3.0, 4.0), m_in1 = NormalWeightedMeanPrecision(1.0, 2.0), m_in2 = PointMass(2.0)),
                output = (in1 = NormalWeightedMeanPrecision(21.0, 6.0), in2 = PointMass(2.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 2.0), m_in1 = NormalMeanPrecision(1.0, 2.0), m_in2 = PointMass(-1.0)),
                output = (in1 = NormalWeightedMeanPrecision(2.0, 4.0), in2 = PointMass(-1.0))
            ),
            (
                input = (m_out = NormalMeanVariance(3.0, 3.0), m_in1 = NormalWeightedMeanPrecision(2.0, 1.0), m_in2 = PointMass(2.0)),
                output = (in1 = NormalWeightedMeanPrecision(11/3, 4/3), in2 = PointMass(2.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 4.0), m_in1 = NormalMeanVariance(2.0, 2.0), m_in2 = PointMass(1.0)),
                output = (in1 = NormalWeightedMeanPrecision(7.0, 9/2), in2 = PointMass(1.0))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass, m_in2::UnivariateNormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = NormalMeanVariance(3.0, 4.0), m_in1 = PointMass(2.0), m_in2 = NormalMeanVariance(2.0, 2.0)),
                output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(5/4, 3/4))
            ),
            (
                input = (m_out = NormalMeanPrecision(3, 4), m_in1 = PointMass(1.0), m_in2 = NormalMeanPrecision(1, 4)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(12.0, 8.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(3.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalWeightedMeanPrecision(1.0, 4.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(0.0, 8.0))
            ),
            (
                input = (m_out = NormalMeanPrecision(2.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(1.0, 2.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(9/2, 9/2))
            ),
            (
                input = (m_out = NormalMeanVariance(2.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanPrecision(1.0, 2.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(9/4, 9/4))
            ),
            (
                input = (m_out = NormalMeanPrecision(3.0, 4.0), m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(1.0, 2.0)),
                output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(5.0, 6.0))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 2.0), m_in1 = PointMass(-1.0), m_in2 = NormalMeanPrecision(1.0, 2.0)),
                output = (in1 = PointMass(-1.0), in2 = NormalWeightedMeanPrecision(6.0, 4.0))
            ),
            (
                input = (m_out = NormalMeanVariance(3.0, 3.0), m_in1 = PointMass(2.0), m_in2 = NormalWeightedMeanPrecision(2.0, 1.0)),
                output = (in1 = PointMass(2.0), in2 = NormalWeightedMeanPrecision(7/3, 4/3))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(2.0, 4.0), m_in1 = PointMass(1.0), m_in2 = NormalMeanVariance(2.0, 2.0)),
                output = (in1 = PointMass(1.0), in2 = NormalWeightedMeanPrecision(-1.0, 9/2))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanCovariance, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanCovariance([2.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0]), m_in2 = PointMass([1.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([5/4, 19/8], [1.0 -1/4; -1/4 7/8]), in2 = PointMass([1.0, 1.0]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0]), m_in2 = PointMass([1.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([5.0, 9/2], [5/2 1.0; 1.0 5/2]), in2 = PointMass([1.0, 1.0]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in2 = PointMass([2.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([13/6, 13/6], [7/3 -1/6; -1/6 7/3]), in2 = PointMass([2.0, 1.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanPrecision, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0]), m_in2 = PointMass([1.0, -1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([13/4, 13/8], [3/2 -1/4; -1/4 11/8]), in2 = PointMass([1.0, -1.0]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 3.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0]), m_in2 = PointMass([2.0, 2.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([15.0, 15.0], [5.0 2.0; 2.0 5.0]), in2 = PointMass([2.0, 2.0]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0]), m_in2 = PointMass([-2.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([11.0, 16.0], [7.0  3.0; 3.0 5.0]), in2 = PointMass([-2.0, 1.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalWeightedMeanPrecision, m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = MvNormalWeightedMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0]), m_in2 = PointMass([1.0, -1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([13/4, 13/8], [3/2 -1/4; -1/4 11/8]), in2 = PointMass([1.0, -1.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0]), m_in2 = PointMass([2.0, 2.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([12.0, 10.0], [5.0 2.0; 2.0 4.0]), in2 = PointMass([2.0, 2.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 2.0], [4.0 1.0; 1.0 4.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 0.0; 0.0 2.0]), m_in2 = PointMass([-1.0, 1.0])),
                output = (in1 = MvNormalWeightedMeanPrecision([-1.0, 6.0], [7.0  1.0; 1.0 6.0]), in2 = PointMass([-1.0, 1.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanCovariance, m_in1::PointMass, m_in2::PointMassMultivariateNormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanCovariance([2.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = PointMass([1.0, 1.0]), m_in2 = MvNormalMeanCovariance([1.0, 3.0], [3.0 2.0; 2.0 4.0])),
                output = (in1 = PointMass([1.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([1/4, 11/8], [1.0 -1/4; -1/4 7/8]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [2.0 0.0; 0.0 2.0]), m_in1 = PointMass([1.0, 1.0]), m_in2 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0])),
                output = (in1 = PointMass([1.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([4.0, 7/2], [5/2 1.0; 1.0 5/2]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = PointMass([2.0, 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 2.0], [2.0 0.0; 0.0 2.0])),
                output = (in1 = PointMass([2.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([7/6, 13/6], [7/3 -1/6; -1/6 7/3]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalMeanPrecision, m_in1::, m_in2::MultivariateNormalDistributionsFamilyPointMass)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = MvNormalMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = PointMass([1.0, -1.0]), m_in2 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0])),
                output = (in1 = PointMass([1.0, -1.0]), in2 = MvNormalWeightedMeanPrecision([5/4, 29/8], [3/2 -1/4; -1/4 11/8]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 3.0]), m_in1 = PointMass([2.0, 2.0]), m_in2 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0])),
                output = (in1 = PointMass([2.0, 2.0]), in2 = MvNormalWeightedMeanPrecision([-1.0, -1.0], [5.0 2.0; 2.0 5.0]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([3.0, 2.0], [4.0 2.0; 2.0 4.0]), m_in1 = PointMass([-2.0, 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0])),
                output = (in1 = PointMass([-2.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([23.0, 16.0], [7.0  3.0; 3.0 5.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::MvNormalWeightedMeanPrecision, m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = MvNormalWeightedMeanPrecision([2.0, 2.0], [1.0 0.0; 0.0 1.0]), m_in1 = PointMass([1.0, -1.0]), m_in2 = MvNormalMeanCovariance([2.0, 3.0], [3.0 2.0; 2.0 4.0])),
                output = (in1 = PointMass([1.0, -1.0]), in2 = MvNormalWeightedMeanPrecision([5/4, 29/8], [3/2 -1/4; -1/4 11/8]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in1 = PointMass([2.0, 2.0]), m_in2 = MvNormalMeanPrecision([1.0, 1.0], [2.0 1.0; 1.0 2.0])),
                output = (in1 = PointMass([2.0, 2.0]), in2 = MvNormalWeightedMeanPrecision([-4.0, -2.0], [5.0 2.0; 2.0 4.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 2.0], [4.0 1.0; 1.0 4.0]), m_in1 = PointMass([-1.0, 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 0.0; 0.0 2.0])),
                output = (in1 = PointMass([-1.0, 1.0]), in2 = MvNormalWeightedMeanPrecision([5.0, 0.0], [7.0  1.0; 1.0 6.0]))
            ),

        ]

    end

    @testset ":in1_in2 (m_out::NormalDistributionsFamily, m_in1::NormalDistributionsFamily, m_in2::NormalDistributionsFamily)" begin

        @test_marginalrules [ with_float_conversions = true ] typeof(-)(:in1_in2) [
            (
                input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = NormalMeanVariance(3.0, 4.0), m_in2 = NormalMeanVariance(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([5/4, 1/3], [3/4 -1/2; -1/2 2/3]))
            ),
            (
                input = (m_out = NormalMeanPrecision(1.0, 2.0), m_in1 = NormalMeanPrecision(3.0, 4.0), m_in2 = NormalMeanPrecision(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([14.0, 28.0], [6.0 -2.0; -2.0 8.0]))
            ),
            (
                input = (m_out = NormalWeightedMeanPrecision(1.0, 2.0), m_in1 = NormalWeightedMeanPrecision(3.0, 4.0), m_in2 = NormalWeightedMeanPrecision(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([4.0, 4.0], [6.0 -2.0; -2.0 8.0]))
            ),
            (
                input = (m_out = NormalMeanVariance(1.0, 2.0), m_in1 = NormalMeanPrecision(3.0, 4.0), m_in2 = NormalWeightedMeanPrecision(5.0, 6.0)),
                output = (MvNormalWeightedMeanPrecision([25/2, 9/2], [9/2 -1/2; -1/2 13/2]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([1.0, 2.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanCovariance([2.0, 3.0], [3.0 1.0; 1.0 2.0]), m_in2 = MvNormalMeanCovariance([1.0, 2.0], [3.0 1.0; 1.0 2.0])),
                output = (MvNormalWeightedMeanPrecision([1/5, 12/5, 0.0, 0.0], [[0.8 -0.4 -0.4 0.2; -0.4 1.2 0.2 -0.6]; [-0.4 0.2 0.8 -0.4; 0.2 -0.6 -0.4 1.2]]))
            ),
            (
                input = (m_out = MvNormalMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 2.0]), m_in2 = MvNormalMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 2.0])),
                output = (MvNormalWeightedMeanPrecision([10.0, 10.0, 0.0, 0.0], [6.0 2.0 -3.0 -1.0; 2.0 4.0 -1.0 -2.0; -3.0 -1.0 6.0 2.0; -1.0 -2.0 2.0 4.0]))
            ),
            (
                input = (m_out = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0]), m_in1 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 1.0; 1.0 1.0])),
                output = (MvNormalWeightedMeanPrecision([2.0, 4.0, 0.0, 0.0], [6.0 2.0 -3.0 -1.0; 2.0 2.0 -1.0 -1.0; -3.0 -1.0 6.0 2.0; -1.0 -1.0 2.0 2.0]))
            ),
            (
                input = (m_out = MvNormalMeanCovariance([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in1 = MvNormalMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0]), m_in2 = MvNormalWeightedMeanPrecision([1.0, 1.0], [3.0 1.0; 1.0 2.0])),
                output = (MvNormalWeightedMeanPrecision([21/5, 17/5, 4/5, 3/5], [17/5 4/5 -2/5 1/5; 4/5 13/5 1/5 -3/5; -2/5 1/5 17/5 4/5; 1/5 -3/5 4/5 13/5]))
            ),

        ]

    end

end
end


text:
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end
parsed_args = parse_args(ARGS, s)

test_suffix = parsed_args["test"] ? "_test" : ""
octopi = [[parse(Int8, d) for d in l] for l in readlines("input/day_11$test_suffix.txt")]
my = length(octopi)
mx = length(octopi[1])

neighborhood = [(-1, -1), (0, -1), (+1, -1), (-1, 0), (1, 0), (-1, 1), (0, 1), (1, 1)]
function step()
    flashers = []
    function increment(x, y)
        if x <= 0 || x > mx || y <= 0 || y > my
            return
        end
        octopi[x][y] += 1
        if octopi[x][y] == 10
            push!(flashers, (x = x, y = y))
        end
    end
    function reset(x, y)
        if octopi[x][y] > 9
            octopi[x][y] = 0
        end
    end

    function all_squares(f)
        for y = 1:my
            for x = 1:mx
                f(x, y)
            end
        end
    end

    all_squares(increment)
    count = 0
    while length(flashers) != 0
        flash = popfirst!(flashers)
        foreach(n -> increment(flash.x + n[1], flash.y + n[2]), neighborhood)
        count += 1
    end
    all_squares(reset)
    count
end

function part1()
    sum(map(n -> step(), 1:100))
end

function part2()
    i = 1
    while true
        count = step()
        if count == mx * my
            return i
        end
        i = i + 1
    end
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
plane = ComplexPlane(Quaternion(rand(), ℝ³(rand(3)))) # Construct a unit Quaternion
w, z = plane.z₁, plane.z₂
B = Complex(10rand(), 10rand())
clifford = Clifford(w, z, B)

@test typeof(clifford.total) <: S³
@test typeof(clifford.base) <: S²
@test isapprox(clifford.total, plane)

cl = ComplexLine(100rand() + im * 100rand())
q = σmap(cl)
h = τmap(cl)

@test typeof(q) <: S³
@test typeof(q) <: S³

q = Quaternion(rand() * 2pi, normalize(ℝ³(rand(3))))
c = U1(rand() * 2pi - pi)
h = S¹action(q, c)

@test typeof(h) <: S³
@test isapprox(norm(q), norm(h))

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
function steig!{T<:AbstractFloat}(d::Array{T,1}, e::Array{T,1}, z::Array{T,1}, maxits::Integer)
    #
    # Finds the eigenvalues and first components of the normalised
    # eigenvectors of a symmetric tridiagonal matrix by the implicit
    # QL method.
    #
    # d[i]   On entry, holds the ith diagonal entry of the matrix.
    #        On exit, holds the ith eigenvalue.
    #
    # e[i]   On entry, holds the [i+1,i] entry of the matrix for
    #        i = 1, 2, ..., n-1.  (The value of e[n] is not used.)
    #        On exit, e is overwritten.
    #
    # z[i]   On exit, holds the first component of the ith normalised
    #        eigenvector associated with d[i].
    #
    # maxits The maximum number of QL iterations.
    #
    # Martin and Wilkinson, Numer. Math. 12: 377-383 (1968).
    # Dubrulle, Numer. Math. 15: 450 (1970).
    # Handbook for Automatic Computation, Vol ii, Linear Algebra,
    #        pp. 241-248, 1971.
    #
    # This is a modified version of the Eispack routine imtql2.
    #
    n = length(z)
    z[1] = 1
    z[2:n] = 0
    e[n] = 0

    if n == 1 # Nothing to do for a 1x1 matrix.
        return
    end
    for l = 1:n
        for j = 1:maxits
            # Look for small off-diagonal elements.
            m = n
            for i = l:n-1
                if abs(e[i]) <= eps(T) * ( abs(d[i]) + abs(d[i+1]) )
                    m = i
                    break
                end
            end
            p = d[l]
            if m == l
                continue
            end
            if j == maxits
                msg = @sprintf("No convergence after %d iterations", j)
                msg *= " (try increasing maxits)"
                error(msg)
            end
            # Form shift
            g = ( d[l+1] - p ) / ( 2 * e[l] )
            r = hypot(g, one(T))
            g = d[m] - p + e[l] / ( g + copysign(r, g) )
            s = one(T)
            c = one(T)
            p = zero(T)
            for i = m-1:-1:l
                f = s * e[i]
                b = c * e[i]
                if abs(f) <  abs(g)
                    s = f / g
                    r = hypot(s, one(T))
                    e[i+1] = g * r
                    c = one(T) / r
                    s *= c
                else
                    c = g / f
                    r = hypot(c, one(T))
                    e[i+1] = f * r
                    s = one(T) / r
                    c *= s
                end
                g = d[i+1] - p
                r = ( d[i] - g ) * s + 2 * c * b
                p = s * r
                d[i+1] = g + p
                g = c * r - b
                # Form first component of vector.
                f = z[i+1]
                z[i+1] = s * z[i] + c * f
                z[i]   = c * z[i] - s * f
            end # loop over i
            d[l] -= p
            e[l] = g
            e[m] = zero(T)
        end # loop over j
    end # loop over l
end


text:
module CVVideoIO

export VideoCapture, isOpened, release

using LibOpenCV
using CVCore
using Cxx

libopencv_videoio = LibOpenCV.find_library_e("libopencv_videoio")
try
    Libdl.dlopen(libopencv_videoio, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

cxx"""
#include <opencv2/videoio.hpp>
"""

include("const.jl")

const cvVideoCapture = cxxt"cv::VideoCapture"
cvVideoCapture(idx::Int) = icxx"return cv::VideoCapture($idx);"

const VideoCapture = cvVideoCapture

function Base.read!(cap::VideoCapture, img::Mat{UInt8})
    icxx"$cap.read($(img.handle));"
end

function Base.read(cap::VideoCapture)
    mat = Mat{UInt8}()
    ok = read!(cap, mat)
    ok, Mat(mat)
end

for f in [
        :grab,
        :isOpened,
        :release
        ]
    body = Expr(:macrocall, Symbol("@icxx_str"), "\$cap.$f();")
    @eval $f(cap::VideoCapture) = $body
end

Base.close(cap::VideoCapture) = release(cap)

set(cap::VideoCapture, propId, value) = icxx"$cap.set($propId, $value);"
get(cap::VideoCapture, propId) = icxx"$cap.get($propId);"

end # module


text:
"""
    SubDataFrame{<:AbstractDataFrame, <:AbstractIndex, <:AbstractVector{Int}} <: AbstractDataFrame

A view of an `AbstractDataFrame`. It is returned by a call to the `view` function
on an `AbstractDataFrame` if a collections of rows and columns are specified.

A `SubDataFrame` is an `AbstractDataFrame`, so expect that most
DataFrame functions should work. Such methods include `describe`,
`summary`, `nrow`, `size`, `by`, `stack`, and `join`.

If the selection of columns in a parent data frame is passed as `:` (a colon)
then `SubDataFrame` will always have all columns from the parent,
even if they are added or removed after its creation.

# Examples
```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8)
8×3 DataFrame
 Row │ a      b      c
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      2      1
   2 │     2      1      2
   3 │     3      2      3
   4 │     4      1      4
   5 │     1      2      5
   6 │     2      1      6
   7 │     3      2      7
   8 │     4      1      8

julia> sdf1 = view(df, :, 2:3) # column subsetting
8×2 SubDataFrame
 Row │ b      c
     │ Int64  Int64
─────┼──────────────
   1 │     2      1
   2 │     1      2
   3 │     2      3
   4 │     1      4
   5 │     2      5
   6 │     1      6
   7 │     2      7
   8 │     1      8

julia> sdf2 = @view df[end:-1:1, [1, 3]]  # row and column subsetting
8×2 SubDataFrame
 Row │ a      c
     │ Int64  Int64
─────┼──────────────
   1 │     4      8
   2 │     3      7
   3 │     2      6
   4 │     1      5
   5 │     4      4
   6 │     3      3
   7 │     2      2
   8 │     1      1

julia> sdf3 = groupby(df, :a)[1]  # indexing a GroupedDataFrame returns a SubDataFrame
2×3 SubDataFrame
 Row │ a      b      c
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      2      1
   2 │     1      2      5
```
"""
struct SubDataFrame{D<:AbstractDataFrame, S<:AbstractIndex, T<:AbstractVector{Int}} <: AbstractDataFrame
    parent::D
    colindex::S
    rows::T # maps from subdf row indexes to parent row indexes
end

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector{Int}, cols)
    @boundscheck if !checkindex(Bool, axes(parent, 1), rows)
        throw(BoundsError(parent, (rows, cols)))
    end
    SubDataFrame(parent, SubIndex(index(parent), cols), rows)
end
Base.@propagate_inbounds SubDataFrame(parent::DataFrame, ::Colon, cols) =
    SubDataFrame(parent, axes(parent, 1), cols)
@inline SubDataFrame(parent::DataFrame, row::Integer, cols) =
    throw(ArgumentError("invalid row index: $row of type $(typeof(row))"))

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector{<:Integer}, cols)
    if any(x -> x isa Bool, rows)
        throw(ArgumentError("invalid row index of type `Bool`"))
    end
    return SubDataFrame(parent, convert(Vector{Int}, rows), cols)
end

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector{Bool}, cols)
    if length(rows) != nrow(parent)
        throw(ArgumentError("invalid length of `AbstractVector{Bool}` row index " *
                            "(got $(length(rows)), expected $(nrow(parent)))"))
    end
    return SubDataFrame(parent, _findall(rows), cols)
end

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector, cols)
    if !all(x -> (x isa Integer) && !(x isa Bool), rows)
        throw(ArgumentError("only `Integer` indices are accepted in `rows`"))
    end
    return SubDataFrame(parent, convert(Vector{Int}, rows), cols)
end

Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind, cols) =
    SubDataFrame(parent(sdf), rows(sdf)[rowind], parentcols(index(sdf), cols))
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind::Bool, cols) =
    throw(ArgumentError("invalid row index of type Bool"))
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind, ::Colon) =
    if index(sdf) isa Index # sdf was created using : as row selector
        SubDataFrame(parent(sdf), rows(sdf)[rowind], :)
    else
        SubDataFrame(parent(sdf), rows(sdf)[rowind], parentcols(index(sdf), :))
    end
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind::Bool, ::Colon) =
    throw(ArgumentError("invalid row index of type Bool"))
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, ::Colon, cols) =
    SubDataFrame(parent(sdf), rows(sdf), parentcols(index(sdf), cols))
@inline SubDataFrame(sdf::SubDataFrame, ::Colon, ::Colon) = sdf

rows(sdf::SubDataFrame) = getfield(sdf, :rows)
Base.parent(sdf::SubDataFrame) = getfield(sdf, :parent)
Base.parentindices(sdf::SubDataFrame) = (rows(sdf), parentcols(index(sdf)))

Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds, colind::ColumnIndex) =
    view(adf[!, colind], rowinds)
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, ::typeof(!), colind::ColumnIndex) =
    view(adf[!, colind], :)
@inline Base.view(adf::AbstractDataFrame, rowinds, colind::Bool) =
    throw(ArgumentError("invalid column index $colind of type `Bool`"))
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds,
                                   colinds::MultiColumnIndex) =
    SubDataFrame(adf, rowinds, colinds)
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds::typeof(!),
                                   colinds::MultiColumnIndex) =
    SubDataFrame(adf, :, colinds)
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds::Not,
                                   colinds::MultiColumnIndex) =
    SubDataFrame(adf, axes(adf, 1)[rowinds], colinds)

##############################################################################
##
## AbstractDataFrame interface
##
##############################################################################

index(sdf::SubDataFrame) = getfield(sdf, :colindex)

nrow(sdf::SubDataFrame) = ncol(sdf) > 0 ? length(rows(sdf))::Int : 0
ncol(sdf::SubDataFrame) = length(index(sdf))

Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowind::Integer, colind::ColumnIndex) =
    parent(sdf)[rows(sdf)[rowind], parentcols(index(sdf), colind)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowind::Bool, colind::ColumnIndex) =
    throw(ArgumentError("invalid row index of type Bool"))

Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowinds::Union{AbstractVector, Not},
                                       colind::ColumnIndex) =
    parent(sdf)[rows(sdf)[rowinds], parentcols(index(sdf), colind)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, ::Colon, colind::ColumnIndex) =
    parent(sdf)[rows(sdf), parentcols(index(sdf), colind)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, ::typeof(!), colind::ColumnIndex) =
    view(parent(sdf), rows(sdf), parentcols(index(sdf), colind))

Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowinds::Union{AbstractVector, Not},
                                       colinds::MultiColumnIndex) =
    parent(sdf)[rows(sdf)[rowinds], parentcols(index(sdf), colinds)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, ::Colon,
                                       colinds::MultiColumnIndex) =
    parent(sdf)[rows(sdf), parentcols(index(sdf), colinds)]
Base.@propagate_inbounds Base.getindex(df::SubDataFrame, row_ind::typeof(!),
                                       col_inds::MultiColumnIndex) =
    select(df, col_inds, copycols=false)

Base.@propagate_inbounds function Base.setindex!(sdf::SubDataFrame, val::Any, idx::CartesianIndex{2})
    setindex!(sdf, val, idx[1], idx[2])
end
Base.@propagate_inbounds function Base.setindex!(sdf::SubDataFrame, val::Any, ::Colon, colinds::Any)
    parent(sdf)[rows(sdf), parentcols(index(sdf), colinds)] = val
    return sdf
end
Base.@propagate_inbounds function Base.setindex!(sdf::SubDataFrame, val::Any, ::typeof(!), colinds::Any)
    throw(ArgumentError("setting index of SubDataFrame using ! as row selector is not allowed"))
end
Base.@propagate_inbounds function Base.setindex!(sdf::SubDataFrame, val::Any, rowinds::Any, colinds::Any)
    parent(sdf)[rows(sdf)[rowinds], parentcols(index(sdf), colinds)] = val
    return sdf
end
Base.@propagate_inbounds Base.setindex!(sdf::SubDataFrame, val::Any, rowinds::Bool, colinds::Any) =
    throw(ArgumentError("invalid row index of type Bool"))

Base.setproperty!(::SubDataFrame, ::Symbol, ::Any) =
    throw(ArgumentError("Replacing or adding of columns of a SubDataFrame is not allowed. " *
                        "Instead use `df[:, col_ind] = v` or `df[:, col_ind] .= v` " *
                        "to perform an in-place assignment."))
Base.setproperty!(::SubDataFrame, ::AbstractString, ::Any) =
    throw(ArgumentError("Replacing or adding of columns of a SubDataFrame is not allowed. " *
                        "Instead use `df[:, col_ind] = v` or `df[:, col_ind] .= v` " *
                        "to perform an in-place assignment."))

##############################################################################
##
## Miscellaneous
##
##############################################################################

Base.copy(sdf::SubDataFrame) = parent(sdf)[rows(sdf), parentcols(index(sdf), :)]

Base.delete!(df::SubDataFrame, ind) =
    throw(ArgumentError("SubDataFrame does not support deleting rows"))

function DataFrame(sdf::SubDataFrame; copycols::Bool=true)
    if copycols
        sdf[:, :]
    else
        DataFrame(collect(eachcol(sdf)), _names(sdf), copycols=false)
    end
end

Base.convert(::Type{DataFrame}, sdf::SubDataFrame) = DataFrame(sdf)


text:
using GridWorld
using Test

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


text:
"""
    Normed{T <: Unsigned, f} <: FixedPoint{T, f}

`Normed{T,f}` maps `Unsigned` integers from `0` to `2^f-1` to the range
[0.0, 1.0]. For example, `Normed{UInt8,8}` maps `0x00` to `0.0` and `0xff` to
`1.0`.

There are the typealiases for `Normed` in the `NXfY` notation, where `Y` is
the number of fractional bits (i.e. `f`), and `X+Y` equals the number of
underlying bits used. For example, `N0f8` is aliased to `Normed{UInt8,8}` and
`N4f12` is aliased to `Normed{UInt16,12}`.
"""
struct Normed{T <: Unsigned, f} <: FixedPoint{T, f}
    i::T

    function Normed{T, f}(i::Integer, _) where {T, f}
        1 <= f <= bitwidth(T) || throw(DomainError(f, "f must be between 1 and $(bitwidth(T)) (i.e. the number of bits of `T=$T`)"))
        new{T, f}(i % T)
    end
end

typechar(::Type{X}) where {X <: Normed} = 'N'

for T in (UInt8, UInt16, UInt32, UInt64)
    io = IOBuffer()
    for f in 1:bitwidth(T)
        sym = Symbol(String(take!(showtype(io, Normed{T,f}))))
        @eval begin
            const $sym = Normed{$T,$f}
            export $sym
        end
    end
end

function rawone(::Type{Normed{T,f}}) where {T <: Unsigned, f}
    typemax(T) >> (bitwidth(T) - f)
end

# constructor-style conversions
function _convert(::Type{N}, x::Normed{T2,f}) where {T, T2, f, N <: Normed{T,f}}
    reinterpret(N, convert(T, x.i)) # TODO: input range checking
end

function _convert(::Type{N}, x::Normed{T2,f2}) where {T, T2, f, f2, N <: Normed{T,f}}
    y = round((rawone(N)/rawone(x))*reinterpret(x))
    (0 <= y) & (y <= typemax(T)) || throw_converterror(N, x)
    reinterpret(N, _unsafe_trunc(T, y))
end

function _convert(::Type{N}, x::Normed{UInt8,8}) where {N <: Normed{UInt16,16}} # TODO: generalization
    reinterpret(N0f16, convert(UInt16, 0x0101*reinterpret(x)))
end

function _convert(::Type{N}, x::Real) where {T, f, N <: Normed{T,f}}
    if T == UInt128 # for UInt128, we can't widen
        # the upper limit is not exact
        (0 <= x) & (x <= (typemax(T)/rawone(N))) || throw_converterror(N, x)
        y = round(rawone(N)*x)
    else
        y = round(widen1(rawone(N))*x)
        (0 <= y) & (y <= typemax(T)) || throw_converterror(N, x)
    end
    reinterpret(N, _unsafe_trunc(T, y))
end
# Prevent overflow (https://discourse.julialang.org/t/saving-greater-than-8-bit-images/6057)
function _convert(::Type{N}, x::Float16) where {T, f, N <: Normed{T,f}}
    if Float16(typemax(T)/rawone(N)) > Float32(typemax(T)/rawone(N))
        x == Float16(typemax(T)/rawone(N)) && return typemax(N)
    end
    return _convert(N, Float32(x))
end
function _convert(::Type{N}, x::Tf) where {T, f, N <: Normed{T,f}, Tf <: Union{Float32, Float64}}
    if T === UInt128 && f == 53
        0 <= x <= Tf(3.777893186295717e22) || throw_converterror(N, x)
    else
        0 <= x <= Tf((typemax(T)-rawone(N))/rawone(N)+1) || throw_converterror(N, x)
    end

    f == 1 && x == Tf(typemax(N)) && return typemax(N)
    if f <= (significand_bits(Tf) + 1) && bitwidth(T) < significand_bits(Tf)
        return reinterpret(N, unsafe_trunc(T, round(rawone(N) * x)))
    end
    # cf. the implementation of `frexp`
    Tw = f < bitwidth(T) ? T : widen1(T)
    bits = bitwidth(Tw) - 1
    xu = reinterpret(Unsigned, x)
    k = Int(xu >> significand_bits(Tf))
    k == 0 && return zero(N) # neglect subnormal numbers
    significand = xu | (oneunit(xu) << significand_bits(Tf))
    yh = unsafe_trunc(Tw, significand) << (bits - significand_bits(Tf))
    ex = exponent_bias(Tf) - k + bits - f
    yi = bits >= f ? yh - (yh >> f) : yh
    if ex <= 0
        ex == 0 && return reinterpret(N, unsafe_trunc(T, yi))
        ex != -1 || signbit(signed(yi)) && return typemax(N)
        return reinterpret(N, unsafe_trunc(T, yi + yi))
    end
    ex > bits && return reinterpret(N, ex == bits + 1 ? oneunit(T) : zero(T))
    yi += oneunit(Tw)<<((ex - 1) & bits) # RoundNearestTiesUp
    return reinterpret(N, unsafe_trunc(T, yi >> (ex & bits)))
end

function _convert(::Type{N}, x::Rational) where {T, f, N <: Normed{T,f}}
    if 0 <= x <= Rational(typemax(N))
        reinterpret(N, round(T, convert(floattype(T), x) * rawone(N)))
    else
        throw_converterror(N, x)
    end
end

_rem(x::N, ::Type{N}) where {N <: Normed} = x
_rem(x::Normed, ::Type{N}) where {T, N <: Normed{T}} =
    reinterpret(N, _unsafe_trunc(T, round((rawone(N)/rawone(x))*reinterpret(x))))
function _rem(x::Real, ::Type{N}) where {T, N <: Normed{T}}
    bitwidth(T) < 32 || isfinite(x) || return zero(N)
    reinterpret(N, _unsafe_trunc(T, round(rawone(N) * x)))
end
_rem(x::Float16, ::Type{X}) where {X <: Normed} = _rem(Float32(x), X)  # avoid overflow
# Float32 and Float64 cannot exactly represent `rawone(N)` with `f` greater than
# the number of their significand bits, resulting in rounding errors (issue #150).
# So, we use another strategy for the large `f`s explained in:
# https://github.com/JuliaMath/FixedPointNumbers.jl/pull/166#issuecomment-574135643
function _rem(x::Float32, ::Type{N}) where {f, N <: Normed{UInt32,f}}
    isfinite(x) || return zero(N)
    f <= 24 && return reinterpret(N, _unsafe_trunc(UInt32, round(rawone(N) * x)))
    r = _unsafe_trunc(UInt32, round(x * @f32(0x1p24)))
    reinterpret(N, r << UInt8(f - 24) - unsigned(signed(r) >> 0x18))
end
function _rem(x::Float64, ::Type{N}) where {f, N <: Normed{UInt64,f}}
    isfinite(x) || return zero(N)
    f <= 53 && return reinterpret(N, _unsafe_trunc(UInt64, round(rawone(N) * x)))
    r = _unsafe_trunc(UInt64, round(x * 0x1p53))
    reinterpret(N, r << UInt8(f - 53) - unsigned(signed(r) >> 0x35))
end

function (::Type{T})(x::Normed) where {T <: AbstractFloat}
    # The following optimization for constant division may cause rounding errors.
    # y = reinterpret(x)*(one(rawtype(x))/convert(T, rawone(x)))
    # Therefore, we use a simple form here.
    # If you prefer speed over accuracy, consider using `scaledual` instead.
    y = reinterpret(x) / convert(promote_type(T, floattype(x)), rawone(x))
    convert(T, y)  # needed for types like Float16 which promote arithmetic to Float32
end

function Base.Float16(x::Normed{Ti,f}) where {Ti <: Union{UInt8, UInt16, UInt32}, f}
    f == 1 ? Float16(x.i) : Float16(Float32(x))
end
function Base.Float16(x::Normed{Ti,f}) where {Ti <: Union{UInt64, UInt128}, f}
    f == 1 ? Float16(x.i) : Float16(Float64(x))
end

function Base.Float32(x::Normed{UInt8,f}) where f
    f == 1 && return Float32(x.i)
    f == 2 && return Float32(Int32(x.i) * 0x101) * @f32(0x550055p-32)
    f == 3 && return Float32(Int32(x.i) * 0x00b) * @f32(0xd4c77bp-30)
    f == 4 && return Float32(Int32(x.i) * 0x101) * @f32(0x110011p-32)
    f == 5 && return Float32(Int32(x.i) * 0x003) * @f32(0xb02c0bp-30)
    f == 6 && return Float32(Int32(x.i) * 0x049) * @f32(0xe40039p-36)
    f == 7 && return Float32(Int32(x.i) * 0x01f) * @f32(0x852b5fp-35)
    f == 8 && return Float32(Int32(x.i) * 0x155) * @f32(0xc0f0fdp-40)
    0.0f0
end
function Base.Float32(x::Normed{UInt16,f}) where f
    f32 = Float32(x.i)
    f ==  1 && return f32
    f ==  2 && return f32 * @f32(0x55p-8)  + f32 * @f32(0x555555p-32)
    f ==  3 && return f32 * @f32(0x49p-9)  + f32 * @f32(0x249249p-33)
    f ==  4 && return f32 * @f32(0x11p-8)  + f32 * @f32(0x111111p-32)
    f ==  5 && return f32 * @f32(0x21p-10) + f32 * @f32(0x108421p-35)
    f ==  6 && return f32 * @f32(0x41p-12) + f32 * @f32(0x041041p-36)
    f ==  7 && return f32 * @f32(0x81p-14) + f32 * @f32(0x204081p-42)
    f == 16 && return f32 * @f32(0x01p-16) + f32 * @f32(0x010001p-48)
    Float32(x.i / rawone(x))
end
function Base.Float32(x::Normed{UInt32,f}) where f
    f == 1 && return Float32(x.i)
    i32 = unsafe_trunc(Int32, x.i)
    if f == 32
        rh, rl = Float32(i32>>>0x10), Float32((i32&0xFFFF)<<0x8 | i32>>>0x18)
        return muladd(rh, @f32(0x1p-16), rl * @f32(0x1p-40))
    elseif f >= 25
        rh, rl = Float32(i32>>>0x10), Float32((i32&0xFFFF)<<0xE + i32>>>UInt8(f-14))
        return muladd(rh, Float32(@exp2(16-f)), rl * Float32(@exp2(-14-f)))
    end
    # FIXME: avoid the branch in native x86_64 (non-SIMD) codes
    m = ifelse(i32 < 0, 0x1p32 * inv_rawone(x), 0.0)
    Float32(muladd(Float64(i32), inv_rawone(x), m))
end
function Base.Float32(x::Normed{Ti,f}) where {Ti <: Union{UInt64, UInt128}, f}
    f == 1 ? Float32(x.i) : Float32(Float64(x))
end

function Base.Float64(x::Normed{Ti,f}) where {Ti <: Union{UInt8, UInt16}, f}
    Float64(Normed{UInt32,f}(x))
end
function Base.Float64(x::Normed{UInt32,f}) where f
    f64 = Float64(x.i)
    f ==  1 && return f64
    f ==  2 && return (f64 * 0x040001p0) * 0x15555000015555p-72
    f ==  3 && return (f64 * 0x108421p0) * 0x11b6db76924929p-75
    f ==  4 && return (f64 * 0x010101p0) * 0x11000011000011p-72
    f ==  5 && return (f64 * 0x108421p0) * 0x04000002000001p-75
    f ==  6 && return (f64 * 0x09dfb1p0) * 0x1a56b8e38e6d91p-78
    f ==  7 && return (f64 * 0x000899p0) * 0x0f01480001e029p-70
    f ==  8 && return (f64 * 0x0a5a5bp0) * 0x18d300000018d3p-80
    f ==  9 && return (f64 * 0x001001p0) * 0x080381c8e3f201p-72
    f == 10 && return (f64 * 0x100001p0) * 0x04010000000401p-80
    f == 11 && return (f64 * 0x000009p0) * 0x0e3aaae3955639p-66
    f == 12 && return (f64 * 0x0a8055p0) * 0x186246e46e4cfdp-84
    f == 13 && return (f64 * 0x002001p0) * 0x10000004000001p-78
    f == 14 && return (f64 * 0x03400dp0) * 0x13b13b14ec4ec5p-84
    f == 15 && return (f64 * 0x000259p0) * 0x06d0c5a4f3a5e9p-75
    f == 16 && return (f64 * 0x011111p0) * 0x00f000ff00fff1p-80
    f == 18 && return (f64 * 0x0b06d1p0) * 0x17377445dd1231p-90
    f == 19 && return (f64 * 0x080001p0) * 0x00004000000001p-76
    f == 20 && return (f64 * 0x000101p0) * 0x0ff010ef10ff01p-80
    f == 21 && return (f64 * 0x004001p0) * 0x01fff8101fc001p-84
    f == 22 && return (f64 * 0x002945p0) * 0x18d0000000018dp-88
    f == 23 && return (f64 * 0x044819p0) * 0x07794a23729429p-92
    f == 27 && return (f64 * 0x000a21p0) * 0x0006518c7df9e1p-81
    f == 28 && return (f64 * 0x00000dp0) * 0x13b13b14ec4ec5p-84
    f == 30 && return (f64 * 0x001041p0) * 0x00fc003f03ffc1p-90
    f == 32 && return (f64 * 0x010101p0) * 0x00ff0000ffff01p-96
    f64 / rawone(x)
end
function Base.Float64(x::Normed{UInt64,f}) where f
    f == 1 && return Float64(x.i)
    if f >= 53
        rh = Float64(unsafe_trunc(Int64, x.i>>0x10)) * @exp2(16-f) # upper 48 bits
        rl = Float64(unsafe_trunc(Int32, x.i&0xFFFF)) * @exp2(-f)  # lower 16 bits
        return rh + muladd(rh, @exp2(-f), rl)
    end
    x.i / rawone(x)
end
function Base.Float64(x::Normed{UInt128,f}) where f
    f == 1 && return Float64(x.i)
    ih, il = unsafe_trunc(Int64, x.i>>0x40), unsafe_trunc(Int64, x.i)
    rh = Float64(ih>>>0x10) * @exp2(f <= 53 ? 80 : 80 - f) # upper 48 bits
    km = @exp2(f <= 53 ? 48 : 48 - f) # for middle 32 bits
    rm = Float64(unsafe_trunc(Int32, ih&0xFFFF)) * (0x1p16 * km) +
         Float64(unsafe_trunc(Int32, il>>>0x30)) * km
    rl = Float64(il&0xFFFFFFFFFFFF) * @exp2(f <= 53 ? 0 : -f) # lower 48 bits
    if f <= 53
        return (rh + (rm + rl)) / unsafe_trunc(Int64, rawone(x))
    elseif f < 76
        return rh + (rm + muladd(rh, @exp2(-f), rl))
    else
        return rh + (rm + rl)
    end
end

Base.BigFloat(x::Normed) = reinterpret(x) / BigFloat(rawone(x))

Base.Rational(x::Normed) = reinterpret(x)//rawone(x)

# Division by `2^f-1` with RoundNearest. The result would be in the lower half bits.
div_2fm1(x::T, ::Val{f}) where {T, f} = (x + (T(1)<<(f - 1) - 0x1)) ÷ (T(1) << f - 0x1)
div_2fm1(x::T, ::Val{1}) where T = x
div_2fm1(x::UInt16,  ::Val{8})  = (((x + 0x80) >> 0x8) + x + 0x80) >> 0x8
div_2fm1(x::UInt32,  ::Val{16}) = (((x + 0x8000) >> 0x10) + x + 0x8000) >> 0x10
div_2fm1(x::UInt64,  ::Val{32}) = (((x + 0x80000000) >> 0x20) + x + 0x80000000) >> 0x20
div_2fm1(x::UInt128, ::Val{64}) = (((x + 0x8000000000000000) >> 0x40) + x + 0x8000000000000000) >> 0x40

# wrapping arithmetic
function wrapping_mul(x::N, y::N) where {T <: Union{UInt8,UInt16,UInt32,UInt64}, f, N <: Normed{T,f}}
    z = widemul(x.i, y.i)
    N(div_2fm1(z, Val(Int(f))) % T, 0)
end

# saturating arithmetic
function saturating_mul(x::N, y::N) where {T <: Union{UInt8,UInt16,UInt32,UInt64}, f, N <: Normed{T,f}}
    f == bitwidth(T) && return wrapping_mul(x, y)
    z = min(widemul(x.i, y.i), widemul(typemax(N).i, rawone(N)))
    N(div_2fm1(z, Val(Int(f))) % T, 0)
end

# checked arithmetic
function checked_mul(x::N, y::N) where {N <: Normed}
    z = float(x) * float(y)
    z < typemax(N) + eps(N)/2 || throw_overflowerror(:*, x, y)
    z % N
end
function checked_mul(x::N, y::N) where {T <: Union{UInt8,UInt16,UInt32,UInt64}, f, N <: Normed{T,f}}
    f == bitwidth(T) && return wrapping_mul(x, y)
    z = widemul(x.i, y.i)
    m = widemul(typemax(N).i, rawone(N)) + (rawone(N) >> 0x1)
    z < m || throw_overflowerror(:*, x, y)
    N(div_2fm1(z, Val(Int(f))) % T, 0)
end

# TODO: decide the default arithmetic for `Normed` mul
# Override the default arithmetic with `checked` for backward compatibility
*(x::N, y::N) where {N <: Normed} = checked_mul(x, y)

# Functions
trunc(x::N) where {N <: Normed} = floor(x)
floor(x::N) where {N <: Normed} = reinterpret(N, x.i - x.i % rawone(N))
function ceil(x::Normed{T,f}) where {T, f}
    f == 1 && return x
    if typemax(T) % rawone(x) != 0
        upper = typemax(T) - typemax(T) % rawone(x)
        x.i > upper && throw_converterror(Normed{T,f}, ceil(T, typemax(x)))
    end
    r = x.i % rawone(x)
    reinterpret(Normed{T,f}, x.i - r + (r > 0 ? rawone(x) : zero(T)))
end
function round(x::Normed{T,f}) where {T, f}
    r = x.i % rawone(x)
    q = rawone(x) - r
    reinterpret(Normed{T,f}, r > q ? x.i + q : x.i - r)
end

trunc(::Type{Ti}, x::Normed) where {Ti <: Integer} = floor(Ti, x)
function floor(::Type{Ti}, x::Normed) where {Ti <: Integer}
    convert(Ti, reinterpret(x) ÷ rawone(x))
end
function ceil(::Type{Ti}, x::Normed) where {Ti <: Integer}
    d, r = divrem(x.i, rawone(x))
    convert(Ti, r > 0 ? d + oneunit(rawtype(x)) : d)
end
function round(::Type{Ti}, x::Normed) where {Ti <: Integer}
    d, r = divrem(x.i, rawone(x))
    convert(Ti, r > (rawone(x) >> 0x1) ? d + oneunit(rawtype(x)) : d)
end

function decompose(x::Normed)
    g = gcd(reinterpret(x), rawone(x))
    div(reinterpret(x),g), 0, div(rawone(x),g)
end

# Range lengths
length(r::AbstractUnitRange{N}) where {N <: Normed{<:UShorterThanInt}} =
    floor(Int, last(r)) - floor(Int, first(r)) + 1
length(r::AbstractUnitRange{N}) where {N <: Normed{T}} where {T<:Unsigned} =
    r.start > r.stop ? T(0) : checked_add(floor(T, last(r)) - floor(T, first(r)), oneunit(T))


text:
#!/usr/bin/env julia
# -*- coding: UTF-8 -*-
# __julia-version__ = 1.7.2
# __author__        = "Max Kannenberg, Martin Scheidt"
# __copyright__     = "2020-2022"
# __license__       = "ISC"
__precompile__(true)

module TrainRuns

## loading standard library packages
using UUIDs, Dates, Statistics
## loading external packages
using YAML, JSONSchema, DataFrames

export
## Interface
trainrun, Train, Path, Settings

## global variables
global g      = 9.80665  # acceleration due to gravity (in m/s^2)
global μ      = 0.2      # friction as constant, todo: implement as function
global Δv_air = 15.0/3.6 # coefficient for velocitiy difference between train and outdoor air (in m/s)

## include package files
include("types.jl")
include("constructors.jl")
include("formulary.jl")
include("calc.jl")
include("characteristics.jl")
include("behavior.jl")
include("output.jl")

## main function
"""
    trainrun(train::Train, path::Path, settings::Settings)

Calculate the running time of a `train` on a `path`.
The `settings` provides the choice of models for the calculation.
`settings` can be omitted. If so, a default is used.
The running time will be return in seconds.

# Examples
```julia-repl
julia> trainrun(train, path)
xxx.xx # in seconds
```
"""
function trainrun(train::Train, path::Path, settings=Settings()::Settings)
    # prepare the input data
    movingSection = determineCharacteristics(path, train, settings)
    settings.outputDetail == :everything && println("The moving section has been prepared.")

    # calculate the train run for oparation mode "minimum running time"
    (movingSection, drivingCourse) = calculateMinimumRunningTime!(movingSection, settings, train)
    settings.outputDetail == :everything && println("The driving course for the shortest running time has been calculated.")

    # accumulate data and create an output dictionary
    output = createOutput(train, settings, path, movingSection, drivingCourse)

    return output
end # function trainrun

end # module TrainRuns


text:
using Hydrogen
using AtomicLevels
using Test

# The Lyman and Balmer series radial dipole moments are given by
# Eq. (63.4)

lyman²(n) = 2^8*n^7*(n-1)^(2n-5)/((n+1)^(2n+5))

balmer_s²(n) = 2^17*n^7*(n^2-1)*(n-2)^(2n-6)/((n+2)^(2n+6))

function balmer_p²(n,ℓ)
    if ℓ == 0
        2^15*n^9*(n-2)^(2n-6)/(3*(n+2)^(2n+6))
    elseif ℓ == 2
        2^19*n^9*(n^2-1)*(n-2)^(2n-7)/(3*(n+2)^(2n+7))
    end
end

function lyman_balmer²(o, n, ℓ)
    # Eq. (BS63.5)
    n == o.n && return (3/2*n*√(n^2-max(ℓ,o.ℓ)^2))^2
    n = float(n)
    if o == o"1s"
        lyman²(n)
    elseif o == o"2s"
        balmer_s²(n)
    elseif o == o"2p"
        balmer_p²(n, ℓ)
    end
end

@testset "Hydrogen.jl" begin
    @testset "Factorial ratios" begin
        for a = 1:6
            for b = 1:6
                @test Hydrogen.factorial_ratio(a,b) == factorial(a)/factorial(b)
            end
        end
    end

    @testset "Radial dipole moments" begin
        # We cannot compare with Table 13 of Bethe and Salpeter 1977,
        # since it provides so few figures and some of the values
        # appear to be wrong. Instead, we compare with the specific
        # formulas for the Lyman and Balmer series.

        for (j,(o,ℓ′)) in enumerate([(o"1s",1), (o"2s",1), (o"2p",0), (o"2p",2)])
            for n′ = 1+ℓ′:8
                @test radial_dipole_moment(o, Orbital(n′, ℓ′))^2 ≈
                    lyman_balmer²(o, n′, ℓ′)
            end
        end
    end

    @testset "Orbitals" begin
        # Error conditions:
        @test_throws ArgumentError non_relativistic_orbital(0, 0)
        @test_throws ArgumentError non_relativistic_orbital(-1, 0)
        @test_throws ArgumentError non_relativistic_orbital(1, 1)
        @test_throws ArgumentError non_relativistic_orbital(2, -1)
        @test_throws MethodError non_relativistic_orbital(2.5, -1)
        # Check values:
        P10 = non_relativistic_orbital(1, 0)
        @test P10(0) ≈ 0.0
    end
end


text:
# Julia IO subclasses are converted into Python objects implementing
# the IOBase + RawIOBase interface.
#
# (A useful model here is the Python FileIO class, which is implemented
#  in Modules/_io/fileio.c in the Python source code.)

##########################################################################
# IOBase methods:

# IO objects should raise IOError for unsupported operations or failed IO
function ioraise(e)
    if isa(e, MethodError) || isa(e, ErrorException)
        ccall((@pysym :PyErr_SetString), Void, (PyPtr, Cstring),
              (pyexc::Dict)[PyIOError],
              string("Julia exception: ", e))
    else
        pyraise(e)
    end
end

macro with_ioraise(expr)
    :(try
        $(esc(expr))
      catch e
        ioraise(e)
      end)
end

function jl_io_readline(io::IO, nb)
    d = readline(io)
    if length(d) > nb
        resize!(d, nb)
    end
    return d
end

function jl_io_readlines(io::IO, nb)
    ret = Any[]  # should it be some other type?
    nread = 0
    while nread < nb && !eof(io)
        d = readline(io)
        nread += length(d)
        push!(ret, PyObject(d))
    end
    return ret
end

isseekable(io) = method_exists(seek, Tuple{typeof(io), Int64})
isseekable(io::IOBuffer) = io.seekable

function jl_io_seek(io::IO, offset, whence)
    if whence == 0
        seek(io, offset)
    elseif whence == 1
        skip(io, offset)
    elseif whence == 2
        seekend(io)
        skip(io, offset)
    else
        throw(ArgumentError("unrecognized whence=$n argument to seek"))
    end
    return position(io)
end

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

pyio_jl(self::PyObject) = unsafe_pyjlwrap_to_objref(self["io"].o)::IO

const PyIO = PyNULL()

pyio_initialized = false
function pyio_initialize()
    global pyio_initialized
    if !pyio_initialized::Bool
        copy!(PyIO, @pydef_object type PyIO
            __init__(self, io::IO; istextio=false) = begin
                self[:io] = pyjlwrap_new(io) # avoid recursion
                self[:istextio] = istextio
            end
            close(self) = @with_ioraise(close(pyio_jl(self)))
            closed.get(self) = @with_ioraise(!isopen(pyio_jl(self)))
            encoding.get(self) = "UTF-8"
            fileno(self) = @with_ioraise(fd(pyio_jl(self)))
            flush(self) = begin
                @with_ioraise begin
                    io = pyio_jl(self)
                    if method_exists(flush, Tuple{typeof(io)})
                        flush(io)
                    end
                end
            end
            isatty(self) = isa(pyio_jl(self), Base.TTY)
            readable(self) = isreadable(pyio_jl(self))
            writable(self) = iswritable(pyio_jl(self))
            readline(self, size=typemax(Int)) =
                @with_ioraise(jl_io_readline(pyio_jl(self), size))
            readlines(self, size=typemax(Int)) =
                @with_ioraise(jl_io_readlines(pyio_jl(self), size))
            seek(self, offset, whence=1) =
                @with_ioraise(jl_io_seek(pyio_jl(self), offset, whence))
            seekable(self) = isseekable(pyio_jl(self))
            tell(self) = @with_ioraise(position(pyio_jl(self)))
            writelines(self, seq) =
                @with_ioraise(for s in seq write(pyio_jl(self), s) end)
            read(self, nb=typemax(Int)) =
                @with_ioraise(self[:istextio] ?
                              String(read(pyio_jl(self), nb)) :
                              read(pyio_jl(self), nb))
            readall(self) =
                @with_ioraise(self[:istextio] ? readstring(pyio_jl(self)) :
                                                read(pyio_jl(self)))
            readinto(self, b) = @with_ioraise(readbytes!(pyio_jl(self), b))
            write(self, b) = @with_ioraise(write(pyio_jl(self), b))
        end)
        pyio_initialized = true
    end
    return
end

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

function PyObject(io::IO)
    pyio_initialize()
    # pyjlwrap_new is necessary to avoid PyIO(io) calling PyObject(::IO)
    PyIO(pyjlwrap_new(io))
end

"""
    PyTextIO(io::IO)
    PyObject(io::IO)

Julia IO streams are converted into Python objects implementing the RawIOBase interface,
so they can be used for binary I/O in Python
"""
function PyTextIO(io::IO)
    pyio_initialize()
    PyIO(pyjlwrap_new(io); istextio=true)
end


text:
# This file was a part of Julia. License is MIT: http://julialang.org/license

import Base: similar, copy, copy!, eltype, push!, pop!, delete!, shift!,
             empty!, isempty, union, union!, intersect, intersect!,
             setdiff, setdiff!, symdiff, symdiff!, in, start, next, done,
             last, length, show, hash, issubset, ==, <=, <, unsafe_getindex,
             unsafe_setindex!, findnextnot, first
if !isdefined(Base, :complement)
    export complement, complement!
else
    import Base: complement, complement!
end

type IntSet
    bits::BitVector
    inverse::Bool
    IntSet() = new(fill!(BitVector(256), false), false)
end
IntSet(itr) = union!(IntSet(), itr)

similar(s::IntSet) = IntSet()
copy(s1::IntSet) = copy!(IntSet(), s1)
function copy!(to::IntSet, from::IntSet)
    resize!(to.bits, length(from.bits))
    copy!(to.bits, from.bits)
    to.inverse = from.inverse
    to
end
eltype(s::IntSet) = Int
sizehint!(s::IntSet, n::Integer) = (_resize0!(s.bits, n+1); s)

# only required on 0.3:
function first(itr::IntSet)
    state = start(itr)
    done(itr, state) && throw(ArgumentError("collection must be non-empty"))
    next(itr, state)[1]
end

# An internal function for setting the inclusion bit for a given integer n >= 0
@inline function _setint!(s::IntSet, n::Integer, b::Bool)
    idx = n+1
    if idx > length(s.bits)
        !b && return s # setting a bit to zero outside the set's bits is a no-op
        newlen = idx + idx>>1 # This operation may overflow; we want saturation
        _resize0!(s.bits, ifelse(newlen<0, typemax(Int), newlen))
    end
    unsafe_setindex!(s.bits, b, idx) # Use @inbounds once available
    s
end

# An internal function to resize a bitarray and ensure the newly allocated
# elements are zeroed (will become unnecessary if this behavior changes)
@inline function _resize0!(b::BitVector, newlen::Integer)
    len = length(b)
    resize!(b, newlen)
    len < newlen && unsafe_setindex!(b, false, len+1:newlen) # resize! gives dirty memory
    b
end

# An internal function that resizes a bitarray so it matches the length newlen
# Returns a bitvector of the removed elements (empty if none were removed)
function _matchlength!(b::BitArray, newlen::Integer)
    len = length(b)
    len > newlen && return splice!(b, newlen+1:len)
    len < newlen && _resize0!(b, newlen)
    return BitVector(0)
end

const _intset_bounds_err_msg = "elements of IntSet must be between 0 and typemax(Int)-1"

function push!(s::IntSet, n::Integer)
    0 <= n < typemax(Int) || throw(ArgumentError(_intset_bounds_err_msg))
    _setint!(s, n, !s.inverse)
end
push!(s::IntSet, ns::Integer...) = (for n in ns; push!(s, n); end; s)

function pop!(s::IntSet)
    s.inverse && throw(ArgumentError("cannot pop the last element of complement IntSet"))
    pop!(s, last(s))
end
function pop!(s::IntSet, n::Integer)
    0 <= n < typemax(Int) || throw(ArgumentError(_intset_bounds_err_msg))
    n in s ? (_delete!(s, n); n) : throw(KeyError(n))
end
function pop!(s::IntSet, n::Integer, default)
    0 <= n < typemax(Int) || throw(ArgumentError(_intset_bounds_err_msg))
    n in s ? (_delete!(s, n); n) : default
end
function pop!(f::Function, s::IntSet, n::Integer)
    0 <= n < typemax(Int) || throw(ArgumentError(_intset_bounds_err_msg))
    n in s ? (_delete!(s, n); n) : f()
end
_delete!(s::IntSet, n::Integer) = _setint!(s, n, s.inverse)
delete!(s::IntSet, n::Integer) = n < 0 ? s : _delete!(s, n)
shift!(s::IntSet) = pop!(s, first(s))

empty!(s::IntSet) = (fill!(s.bits, false); s.inverse = false; s)
isempty(s::IntSet) = s.inverse ? length(s.bits) == typemax(Int) && all(s.bits) : !any(s.bits)

# Mathematical set functions: union!, intersect!, setdiff!, symdiff!
# When applied to two intsets, these all have a similar form:
# - Reshape s1 to match s2, occasionally grabbing the bits that were removed
# - Use map to apply some bitwise operation across the entire bitvector
#   - These operations use functors to work on the bitvector chunks, so are
#     very efficient... but a little untraditional. E.g., (p > q) => (p & ~q)
# - If needed, append the removed bits back to s1 or invert the array

union(s::IntSet, ns) = union!(copy(s), ns)
union!(s::IntSet, ns) = (for n in ns; push!(s, n); end; s)
function union!(s1::IntSet, s2::IntSet)
    l = length(s2.bits)
    if     !s1.inverse & !s2.inverse;  e = _matchlength!(s1.bits, l); map!(|, s1.bits, s1.bits, s2.bits); append!(s1.bits, e)
    elseif  s1.inverse & !s2.inverse;  e = _matchlength!(s1.bits, l); map!(>, s1.bits, s1.bits, s2.bits); append!(s1.bits, e)
    elseif !s1.inverse &  s2.inverse;  _resize0!(s1.bits, l);         map!(<, s1.bits, s1.bits, s2.bits); s1.inverse = true
    else #= s1.inverse &  s2.inverse=# _resize0!(s1.bits, l);         map!(&, s1.bits, s1.bits, s2.bits)
    end
    s1
end

intersect(s1::IntSet) = copy(s1)
intersect(s1::IntSet, ss...) = intersect(s1, intersect(ss...))
function intersect(s1::IntSet, ns)
    s = IntSet()
    for n in ns
        n in s1 && push!(s, n)
    end
    s
end
intersect(s1::IntSet, s2::IntSet) = intersect!(copy(s1), s2)
function intersect!(s1::IntSet, s2::IntSet)
    l = length(s2.bits)
    if     !s1.inverse & !s2.inverse;  _resize0!(s1.bits, l);         map!(&, s1.bits, s1.bits, s2.bits)
    elseif  s1.inverse & !s2.inverse;  _resize0!(s1.bits, l);         map!(<, s1.bits, s1.bits, s2.bits); s1.inverse = false
    elseif !s1.inverse &  s2.inverse;  e = _matchlength!(s1.bits, l); map!(>, s1.bits, s1.bits, s2.bits); append!(s1.bits, e)
    else #= s1.inverse &  s2.inverse=# e = _matchlength!(s1.bits, l); map!(|, s1.bits, s1.bits, s2.bits); append!(s1.bits, e)
    end
    s1
end

setdiff(s::IntSet, ns) = setdiff!(copy(s), ns)
setdiff!(s::IntSet, ns) = (for n in ns; _delete!(s, n); end; s)
function setdiff!(s1::IntSet, s2::IntSet)
    l = length(s2.bits)
    if     !s1.inverse & !s2.inverse;  e = _matchlength!(s1.bits, l); map!(>, s1.bits, s1.bits, s2.bits); append!(s1.bits, e)
    elseif  s1.inverse & !s2.inverse;  e = _matchlength!(s1.bits, l); map!(|, s1.bits, s1.bits, s2.bits); append!(s1.bits, e)
    elseif !s1.inverse &  s2.inverse;  _resize0!(s1.bits, l);         map!(&, s1.bits, s1.bits, s2.bits)
    else #= s1.inverse &  s2.inverse=# _resize0!(s1.bits, l);         map!(<, s1.bits, s1.bits, s2.bits); s1.inverse = false
    end
    s1
end

symdiff(s::IntSet, ns) = symdiff!(copy(s), ns)
symdiff!(s::IntSet, ns) = (for n in ns; symdiff!(s, n); end; s)
function symdiff!(s::IntSet, n::Integer)
    0 <= n < typemax(Int) || throw(ArgumentError(_intset_bounds_err_msg))
    val = (n in s) ⊻ !s.inverse
    _setint!(s, n, val)
    s
end
function symdiff!(s1::IntSet, s2::IntSet)
    e = _matchlength!(s1.bits, length(s2.bits))
    map!(⊻, s1.bits, s1.bits, s2.bits)
    s2.inverse && (s1.inverse = !s1.inverse)
    append!(s1.bits, e)
    s1
end

function in(n::Integer, s::IntSet)
    idx = n+1
    if 1 <= idx <= length(s.bits)
        unsafe_getindex(s.bits, idx) != s.inverse
    else
        ifelse((idx <= 0) | (idx > typemax(Int)), false, s.inverse)
    end
end

# Use the next-set index as the state to prevent looking it up again in done
start(s::IntSet) = next(s, 0)[2]
function next(s::IntSet, i, invert=false)
    if s.inverse ⊻ invert
        # i+1 could rollover causing a BoundsError in findnext/findnextnot
        nextidx = i == typemax(Int) ? 0 : findnextnot(s.bits, i+1)
        # Extend indices beyond the length of the bits since it is inverted
        nextidx = nextidx == 0 ? max(i, length(s.bits))+1 : nextidx
    else
        nextidx = i == typemax(Int) ? 0 : findnext(s.bits, i+1)
    end
    (i-1, nextidx)
end
done(s::IntSet, i) = i <= 0

# Nextnot iterates through elements *not* in the set
nextnot(s::IntSet, i) = next(s, i, true)

function last(s::IntSet)
    l = length(s.bits)
    if s.inverse
        idx = l < typemax(Int) ? typemax(Int) : findprevnot(s.bits, l)
    else
        idx = findprev(s.bits, l)
    end
    idx == 0 ? throw(ArgumentError("collection must be non-empty")) : idx - 1
end

length(s::IntSet) = (n = sum(s.bits); ifelse(s.inverse, typemax(Int) - n, n))

complement(s::IntSet) = complement!(copy(s))
complement!(s::IntSet) = (s.inverse = !s.inverse; s)

function show(io::IO, s::IntSet)
    print(io, "IntSet([")
    first = true
    for n in s
        if s.inverse && n > 2 && done(s, nextnot(s, n-3)[2])
             print(io, ", ..., ", typemax(Int)-1)
             break
         end
        !first && print(io, ", ")
        print(io, n)
        first = false
    end
    print(io, "])")
end

function ==(s1::IntSet, s2::IntSet)
    l1 = length(s1.bits)
    l2 = length(s2.bits)
    l1 < l2 && return ==(s2, s1) # Swap so s1 is always equal-length or longer

    # Try to do this without allocating memory or checking bit-by-bit
    if s1.inverse == s2.inverse
        # If the lengths are the same, simply punt to bitarray comparison
        l1 == l2 && return s1.bits == s2.bits
        # Otherwise check the last bit. If equal, we only need to check up to l2
        return findprev(s1.bits, l1) == findprev(s2.bits, l2) &&
               unsafe_getindex(s1.bits, 1:l2) == s2.bits
    else
        # one complement, one not. Could feasibly be true on 32 bit machines
        # Only if all non-overlapping bits are set and overlaps are inverted
        return l1 == typemax(Int) &&
               map!(!, unsafe_getindex(s1.bits, 1:l2)) == s2.bits &&
               (l1 == l2 || all(unsafe_getindex(s1.bits, l2+1:l1)))
    end
end

const hashis_seed = UInt === UInt64 ? 0x88989f1fc7dea67d : 0xc7dea67d
function hash(s::IntSet, h::UInt)
    # Only hash the bits array up to the last-set bit to prevent extra empty
    # bits from changing the hash result
    l = findprev(s.bits, length(s.bits))
    hash(unsafe_getindex(s.bits, 1:l), h) ⊻ hash(s.inverse) ⊻ hashis_seed
end

issubset(a::IntSet, b::IntSet) = isequal(a, intersect(a,b))
<(a::IntSet, b::IntSet) = (a<=b) && !isequal(a,b)
<=(a::IntSet, b::IntSet) = issubset(a, b)


text:
using JLD

include("perfplots.jl")

timeshift = 10.0

# produced by process_csv.jl in perfprofile mode
jldsource = joinpath(dirname(@__FILE__),"..",ARGS[1],"misocp_perf.jld")

d = load(jldsource)
solvers = d["solvers"]
@assert solvers == ["BONMIN","PAJ_CBC_ECOS","PAJ_GLPK_ECOS","SCIP_MISOCP","CPLEX_MISOCP","PAJ_CPLEX_MOSEK","PAJ_CPLEX_MOSEK_msd"]
instance_names = d["instance_names"]

comm_solvers = solvers[[5,7]]
comm_time_table = d["time_table"][:,[5,7]] + timeshift

open_solvers = solvers[[1,2]]
open_time_table = d["time_table"][:,[1,2]] + timeshift

performance_profile(comm_time_table, comm_solvers, logscale=false, ymax=1.0, xmax=8, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"..",ARGS[1],"misocp_comm_time.tex"))

performance_profile(open_time_table, open_solvers, logscale=false, ymax=1.0, xmax=8, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"..",ARGS[1],"misocp_open_time.tex"))

max_stars_repo_name: chriscoey/PajaritoSupplement

max_issues_repo_name: chriscoey/PajaritoSupplement

max_forks_repo_name: chriscoey/PajaritoSupplement


text:
using Socp: Problem, State, Cone, POC, SOC, vprod, iprod, make_e, SparseSolver, DenseSolver
using Socp: compute_scaling, cgt, setup_iter, solve_kkt
using Socp: deg, max_step, solve_socp, Scaling, SqrScaling, scale!, iscale!, SolverState, modify_factors!, compute_full_scaling
using Test
using LinearAlgebra
using SparseArrays
using SuiteSparse
using SuiteSparse.CHOLMOD

@testset "Vector operations" begin
    tv1 = Float64[1,1,1,1,2,3]
    tv2 = Float64[1,1,1,1,5,6]
    tcone = (POC(0,3), SOC(3,3))
    tid = make_e(tcone)
    @test vprod(tcone, tv1, tv1) == [1,1,1,14,4,6]
    @test vprod(tcone, tv1, tv2) == [1,1,1,29,7,9]
    @test norm(vprod(tcone, tv1, iprod(tcone, tv1, tv2)) .- tv2) < 0.0001
    @test vprod(tcone, tid, tv1) == tv1
    @test vprod(tcone, tid, tv2) == tv2
    @test cgt(POC(0,3), [1,2,3], [0,0,0])
    @test !cgt(POC(0,3), [-1,2,3], [0,0,0])
    @test cgt(SOC(0,3), [3,2,2], [0,0,0])
    @test !cgt(SOC(0,3), [2,2,2], [0,0,0])

    @test max_step(POC(0,3), [1,2,3]) == -1
    @test max_step(SOC(0,3), [1,2,3]) == sqrt(2^2 + 3^2)-1.0
    @test max_step((POC(0,3), SOC(0,3)), [1,2,3]) == sqrt(2^2 + 3^2)-1.0
end

@testset "Nesterov-Todd Scalings" begin
    tv1 = Float64[1,1,1,9,2,3]
    tv2 = Float64[1,1,1,22,5,6]
    tcone = (POC(0,3), SOC(3,3))
    s = Scaling(zeros(6,6),zeros(6,6),zeros(6,6),zeros(6),tcone)
    compute_scaling(tcone, s, tv1, tv2)
    sca,isca,pt = s.W,s.iW,s.l
    tvp1 = isca' * tv1
    tvp2 = sca * tv2
    @test sum(isca*sca .- Matrix{Float64}(I, size(sca)...)) < 0.001
    @test norm(isca' * tv1 .- sca * tv2) < 0.001
    @test norm(isca' * tv1 .- pt) < 0.001
    op,op2 = similar(tv2),similar(tv2)
    scale!(tcone, s, tv2, op)
    @test norm(sca*tv2 .- op) < 0.001
    iscale!(tcone, s, op, op2)
    println(tv2 .- op2)
    @test norm(tv2 .- op2) < 0.001
end

@testset "Squared NT Scalings" begin
    tv1 = Float64[1,1,1,9,2,3]
    tv2 = Float64[1,1,1,22,5,6]
    tcone = (POC(0,3), SOC(3,3))
    s = Scaling(zeros(6,6),zeros(6,6),zeros(6,6),zeros(6),tcone)
    compute_scaling(tcone, s, tv1, tv2)
    s2 = SqrScaling(Diagonal(zeros(6)), Diagonal(zeros(6)), zeros(6), tcone, 3, sparse(I, 6, 6))
    compute_scaling(tcone, s2, tv1, tv2)
    @test norm(s.l .- s2.l) < 0.0001
    @test norm(s.wbs .- s2.wbs) < 0.0001
    @test norm(s.mu .- s2.mu) < 0.0001
    #@test norm(norm(s.iWiW .- s2.iWiW)) < 0.01 # this doesn't hold because we've factored out the low rank terms

	G = [0 0 1.0; 0 0 -1; 0 -1 0; -1 0 0]
	cones = (POC(0,1), SOC(1,3))
    u=[3.414213562373095, 2.414213562373095, 1.0, 1.0]
    v=[1.414213562373095, 2.414213562373095, -1.0, -1.0]
    s = Scaling(zeros(4,4),zeros(4,4),zeros(4,4),zeros(4),cones)
    compute_scaling(cones, s, u, v)
    s2 = SqrScaling(Diagonal(zeros(4)), Diagonal(zeros(4)), zeros(4), cones, 3, G)
    compute_scaling(cones, s2, u, v)
    @test norm(s.l .- s2.l) < 0.0001
    @test norm(s.wbs .- s2.wbs) < 0.0001
    @test norm(s.mu .- s2.mu) < 0.0001
    f1 = cholesky(Sparse(sparse(G' * s.iWiW * G)))
    f2 = cholesky(Sparse(sparse(G' * s2.iWiW * G)))
    modify_factors!(cones, f2, s2, sparse(G))
    @test norm(norm(f2\Matrix(I, 3, 3) .- f1\Matrix(I, 3, 3))) < 0.01

    fs = compute_full_scaling(cones, s2)
    @test norm(norm(fs .- s.iWiW)) < 0.01

    u=[2.5571536140045033, 2.594521365784194, 1.6419234525608073, 1.6419234525608069]
    v=[0.42421356237309515, 1.424213562373095, -1.0, -0.9999999999999997]
    s = Scaling(zeros(4,4),zeros(4,4),zeros(4,4),zeros(4),cones)
    compute_scaling(cones, s, u, v)
    s2 = SqrScaling(Diagonal(zeros(4)), Diagonal(zeros(4)), zeros(4), cones, 3, G)
    compute_scaling(cones, s2, u, v)
    fs = compute_full_scaling(cones, s2)
    @test norm(norm(fs .- s.iWiW)) < 0.01
    @test norm(norm(s2.iW*s2.iW .- s2.iWiW)) < 0.01

    #@test norm(norm(s.iWiW .- s2.iWiW)) < 0.01
end

@testset "KKT reference solution" begin

	A = Array{Float64,2}(UndefInitializer(), 0, 3)
	b = Array{Float64,1}(UndefInitializer(), 0)

	c = [-1.0,-1.0,1.0]
	G = [0 0 1.0; 0 0 -1; 0 -1 0; -1 0 0]
	h = [5.0,0.0,0.0,0.0]
	cones = (POC(0,1), SOC(1,3))
	x=[3.5093936289670493, 3.5093936289670546, 4.984484369850202]
	y=Float64[]
	z=[0.4416936891219317, 1.4416936891219376, -1.0000000000000004, -0.9999999999999993]
	s=[0.025571536140045037, 4.994540275840449, 3.5093936289670546, 3.5093936289670484]
	dx=[6.661338147750939e-16, -4.440892098500626e-16, 5.995204332975845e-15]
	dy=Float64[]
	dz=[-0.010055905990247638, -0.01005590599024675, -0.0, 8.881784197001252e-16]
	ds=[-0.011294786134211294, -0.18180993781041443, -0.06493037168608469, -0.06493037168608427]
	cxr=[0.02479569536244916, 0.02479569536250157, 0.013918468357446631]
	cyr=Float64[]
	czr=[-0.02758755986830461, -0.02758755986827908, 1.3092496892856528e-14, 3.9336377471191126e-14]
	csr=[-0.02397437434769427, 0.003862562367179831, 0.024795695362486953, 0.02479569536243503]
	cx,cy,cz,cs = zeros(size(cxr)), zeros(size(cyr)), zeros(size(czr)), zeros(size(csr))
	mehrotra=true
	prob = Problem(c, A, b, G, h, cones)
    s2 = SqrScaling(Diagonal(zeros(4)), Diagonal(zeros(4)), zeros(4), cones, 3, G)
    compute_scaling(cones, s2, s, z)
    spsolver = SparseSolver(prob)
    setup_iter(spsolver, prob, State(prob, x,y,z,s), s2)
    solve_kkt(spsolver, prob, State(prob, x,y,z,s), s2, dx, dy, dz, ds, cx, cy, cz, cs)
    @test norm(cx .- cxr) < 0.001
    @test norm(cy .- cyr) < 0.001
    @test norm(cz .- czr) < 0.001
    @test norm(cs .- csr) < 0.001
end

@testset "SOC programming 1" begin
	c = [-1.0,-1.0,1.0]

	A = Array{Float64,2}(UndefInitializer(), 0, 3)
	b = Array{Float64,1}(UndefInitializer(), 0)

	G = [0 0 1.0; 0 0 -1; 0 -1 0; -1 0 0]
	h = [5.0,0.0,0.0,0.0]
	cones = (POC(0,1), SOC(1,3))
	prob = Problem(c, A, b, G, h, cones)
	ss = SolverState(prob, SparseSolver(prob))
	soln = solve_socp(prob, ss)
	@test norm(soln.x .- [3.53553, 3.53553, 5.0]) < 0.001
	ssd = SolverState(prob, DenseSolver(prob))
	soln = solve_socp(prob, ss)
	@test norm(soln.x .- [3.53553, 3.53553, 5.0]) < 0.001
end

@testset "SOC programming 2" begin
	c = [-2.0, 1.0, 5.0]

	A = Array{Float64,2}(UndefInitializer(), 0, 3)
	b = Array{Float64,1}(UndefInitializer(), 0)

	G = [12.0 6.0 -5.0;
		 13.0 -3.0 -5.0;
		 12.0 -12.0 6.0;
		 3.0 -6.0 10.0;
		 3.0 -6.0 -2.0;
		 -1.0 -9.0 -2.0;
		 1.0 19.0 -3.0]
	h = [-12.0, -3.0, -2.0, 27.0, 0.0, 3.0, -42.0]
	cones = (SOC(0,3), SOC(3,4))
	prob = Problem(c, A, b, G, h, cones)
	ss = SolverState(prob, SparseSolver(prob))
	soln = solve_socp(prob, ss)
	@test norm(soln.x .- [-5.01467, -5.7669, -8.52176]) < 0.001
end

@testset "SOC programming 3" begin
	c = [-2.0, 1.0, 5.0]

	A = [1.0 0.0 0.0] # Array{Float64,2}(UndefInitializer(), 0, 3)
	b = [-3.0] # Array{Float64,1}(UndefInitializer(), 0)

	G = [12.0 6.0 -5.0;
		 13.0 -3.0 -5.0;
		 12.0 -12.0 6.0;
		 3.0 -6.0 10.0;
		 3.0 -6.0 -2.0;
		 -1.0 -9.0 -2.0;
		 1.0 19.0 -3.0]
	h = [-12.0, -3.0, -2.0, 27.0, 0.0, 3.0, -42.0]
	cones = (SOC(0,3), SOC(3,4))
	prob = Problem(c, A, b, G, h, cones)
	ss = SolverState(prob, SparseSolver(prob))
	soln = solve_socp(prob, ss)
	@test norm(soln.x .- [-3.0, -4.82569, -6.64011]) < 0.001
	ssd = SolverState(prob, DenseSolver(prob))
	soln = solve_socp(prob, ss)
	@test norm(soln.x .- [-3.0, -4.82569, -6.64011]) < 0.001
end

# problem: mass 1 travelling in +x direction with velocity 1 starting from position 0
# objective: min applied force for velocity of 0 at position 0
# state: velocity, position
# velocity_0 = 1, position_0 = 0
# velocity_i = force_i-1 + velocity_i-1
# position_i = velocity_i-1 + position_i-1
# velocity_n = 0, position_0 = 0
# |force_i| <= tot
# n vel, n pos, n-1 force, total
# x = [vel..., pos..., force...,tot] length 3n

@testset "Linear optimal control" begin
	n = 50;
	c = [zeros(3n-1); 1.0] # min tot ;
	vel = 1:n; pos = n+1:2*n; force = 2*n+1:3*n-1;
	A = zeros(2*n+2, 3*n);
	b = zeros(2*n+2);
	A[vel[1],vel[1]] = 1.0 # vel[1] = 1.0;
	b[vel[1]] = 1.0;
	A[pos[1],pos[1]] = 1.0 # pos[1] = 0.0;
	b[pos[1]] = 0.0;
	for stp=2:n;
		A[vel[stp], vel[stp]] = -1.0 # vel[stp] - vel[stp-1] - force[stp-1] = 0.0;
		A[vel[stp], vel[stp-1]] = 1.0;
		A[vel[stp], force[stp-1]] = 1.0;
		b[vel[stp]] = 0.0;

		A[pos[stp], pos[stp]] = -1.0;
		A[pos[stp], pos[stp-1]] = 1.0;
		A[pos[stp], vel[stp-1]] = 1.0;
		b[pos[stp]] = 0.0;
	end;
	A[2*n+1, vel[n]] = 1.0 # vel[n] = 0.0;
	b[2*n+1] = 0.0;
	A[2*n+2, pos[n]] = 1.0;
	b[2*n+2] = 0.0;

	# G = [0 0 0 -1.0; 0 0 -I 0];
	G = zeros(n, 3*n);
	G[1,3*n] = -1.0;
	for i = 1:n-1;
		G[i+1,2*n+i] = -1.0;
	end	;
	h = zeros(n);
	cones = (SOC(0,n),);
	prob = Problem(c, A, b, G, h, cones);
	ss = SolverState(prob, SparseSolver(prob));
	soln = solve_socp(prob, ss);
	println(soln.x[force])
	println("Performance evaluation")
	@time for i=1:100 soln = solve_socp(prob, ss) end
end


text:
using Documenter
using DataDeps

using Documenter

makedocs(
    format = :html,
    sitename = "DataDeps.jl",
    modules = [DataDeps]
)

deploydocs(
    repo = "github.com/oxinabox/DataDeps.jl.git",
    julia  = "0.7",
    latest = "master",
    target = "build",
    deps = nothing,  # we use the `format = :html`, without `mkdocs`
    make = nothing,  # we use the `format = :html`, without `mkdocs`
)


text:
#!/usr/bin/env julia

using FeedbackParticleFilters, Test, SafeTestsets

println("Starting tests")

@safetestset "Basics" begin
    println("")
    println("")
    println("")
    printstyled("Basics", bold=true)
    println("")
    include("Basics/test_Abstractions.jl")
    include("Basics/test_AbstractModel.jl")
    include("Basics/test_AbstractFilteringProblem.jl")
    include("Basics/test_AbstractFilterRepresentation.jl")
    include("Basics/test_AbstractFilteringAlgorithm.jl")
    include("Basics/test_HiddenStateModel.jl")
    include("Basics/test_ObservationModel.jl")
    include("Basics/test_ParametricRepresentation.jl")
    include("Basics/test_ParticleRepresentation.jl")
    include("Basics/test_UnweightedParticleRepresentation.jl")
    include("Basics/test_WeightedParticleRepresentation.jl")
    include("Basics/test_GainEstimation.jl")
    include("Basics/test_FilteringProblem.jl")
end

@safetestset "State models" begin
    println("")
    println("")
    println("")
    printstyled("State models", bold=true)
    println("")
    include("StateModels/test_DiffusionStateModel.jl")
    include("StateModels/test_LinearDiffusionStateModel.jl")
    include("StateModels/test_BrownianMotionTorus.jl")
end

@safetestset "Observation models" begin
    println("")
    println("")
    println("")
    printstyled("Observation models", bold=true)
    println("")
    include("ObservationModels/test_DiffusionObservationModel.jl")
    include("ObservationModels/test_LinearDiffusionObservationModel.jl")
end

@safetestset "Filter representations" begin
    println("")
    println("")
    println("")
    printstyled("Filter representations", bold=true)
    println("")
    include("FilterRepresentations/test_MultivariateGaussian.jl")
    include("FilterRepresentations/test_UnweightedParticleEnsemble.jl")
    include("FilterRepresentations/test_WeightedParticleEnsemble.jl")
end

@safetestset "Gain equations" begin
    println("")
    println("")
    println("")
    printstyled("Gain equations", bold=true)
    println("")
    include("GainEquations/test_PoissonEquation.jl")
end

@safetestset "Gain estimation methods" begin
    println("")
    println("")
    println("")
    printstyled("Gain estimation methods", bold=true)
    println("")
    include("GainEstimationMethods/test_ConstantGainApproximation.jl")
    include("GainEstimationMethods/test_ConstantGainEKSPF.jl")
    include("GainEstimationMethods/test_SemigroupMethod.jl")
    include("GainEstimationMethods/test_DifferentialRKHSMethod.jl")
    include("GainEstimationMethods/test_DifferentialRKHSMethodS1.jl")
end

@safetestset "Filtering algorithms" begin
    println("")
    println("")
    println("")
    printstyled("Filtering algorithms", bold=true)
    println("")
    include("FilteringAlgorithms/test_BPF.jl")
    include("FilteringAlgorithms/test_KBF.jl")
    include("FilteringAlgorithms/test_FPF.jl")
    include("FilteringAlgorithms/test_ppFPF.jl")
end

@safetestset "Simulation" begin
    println("")
    println("")
    println("")
    printstyled("Simulation", bold=true)
    println("")
    include("Simulation/test_Simulation.jl")
    include("Simulation/test_SimulationState.jl")
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
#!/usr/bin/env julia

#=
This node listens to and records messages of following topics:
Raw GPS (/hedge_pos)
Raw IMU (/imu/data)
Commands (/ecu)
Raw velocity estimation (/vel_est)
Position info (/pos_info)
=#

using RobotOS
@rosimport barc.msg: ECU, pos_info, Vel_est
#@rosimport data_service.msg: TimeData
@rosimport geometry_msgs.msg: Vector3
@rosimport sensor_msgs.msg: Imu
@rosimport marvelmind_nav.msg: hedge_pos
rostypegen()
using barc.msg
#using data_service.msg
using geometry_msgs.msg
using sensor_msgs.msg
using std_msgs.msg
using marvelmind_nav.msg
using JLD

include("barc_lib/log_functions.jl")

function main()

    buffersize      = 60000
    gps_meas        = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
    imu_meas        = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,9))
    cmd_log         = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
    cmd_pwm_log     = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
    vel_est_log     = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,5))
    pos_info_log    = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,20))

    # initiate node, set up publisher / subscriber topics
    init_node("barc_record")
    s1  = Subscriber("ecu", ECU, ECU_callback, (cmd_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.ECU}
    s2  = Subscriber("imu/data", Imu, IMU_callback, (imu_meas,), queue_size=1)::RobotOS.Subscriber{sensor_msgs.msg.Imu}
    s3  = Subscriber("ecu_pwm", ECU, ECU_PWM_callback, (cmd_pwm_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.ECU}
    s4  = Subscriber("hedge_pos", hedge_pos, GPS_callback, (gps_meas,), queue_size=30)::RobotOS.Subscriber{marvelmind_nav.msg.hedge_pos}
    s5  = Subscriber("pos_info", pos_info, pos_info_callback, (pos_info_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.pos_info}
    s6  = Subscriber("vel_est", Vel_est, vel_est_callback, (vel_est_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.Vel_est}

    run_id = get_param("run_id")

    println("Recorder running.")
    spin()                              # wait for callbacks until shutdown

    # Clean up buffers
    clean_up(gps_meas)
    clean_up(imu_meas)
    clean_up(cmd_log)
    clean_up(cmd_pwm_log)
    clean_up(pos_info_log)
    clean_up(vel_est_log)

    # Save simulation data to file
    #log_path = "$(homedir())/simulations/record-$(Dates.format(now(),"yyyy-mm-dd-HH-MM-SS")).jld"
    log_path = "$(homedir())/simulations/output-record-$(run_id[1:4]).jld"
    if isfile(log_path)
        log_path = "$(homedir())/simulations/output-record-$(run_id[1:4])-2.jld"
        warn("Warning: File already exists.")
    end
    save(log_path,"gps_meas",gps_meas,"imu_meas",imu_meas,"cmd_log",cmd_log,"cmd_pwm_log",cmd_pwm_log,"pos_info",pos_info_log,"vel_est",vel_est_log)
    println("Exiting node... Saving recorded data to $log_path.")
end

if ! isinteractive()
    main()
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
mutable struct C2
    Wc1::AbstractArray{T,4} where T
    Wc2::AbstractArray{T,4} where T
    W1::AbstractMatrix
    b::AbstractVector
end

C2() = C2(rand(1,1,1,1), rand(1,1,1,1), rand(1,1), rand(1))

function predict(m, X::AbstractArray)
    c1 = conv2d(X, m.Wc1)
    p1 = pool(c1)
    c2 = conv2d(p1, m.Wc2)
    p2 = pool(c2)
    cc = reshape(p2, 4*4*6, size(p1, 4))
    res = sigmoid.(m.W1 * cc .+ m.b)
    return res
end

function loss(m, X::AbstractArray, y::AbstractMatrix)
    ŷ = predict(m, X)
    # cost = sum((y .- ŷ) .^ 2)
    cost = pytorch_cross_entropy(ŷ, y)
    return cost
end

function partial_fit!(m, X::AbstractArray, y::AbstractMatrix; mem=Dict())
    # opt = ModelOptimizer(m, Momentum(η=0.001, γ=0.9))
    opt = ModelOptimizer(m, Adam(α=0.001))
    cost, dm, dX, dy = xgrad(loss; mem=mem, m=m, X=X, y=y)
    update_params!(opt, m, dm)
    return cost
end

function fit!(m, X::AbstractArray, Y::AbstractMatrix; cuda=false, n_epochs=100)
    mem = Dict()
    if cuda
        m = to_cuda(m)
    end
    @time for epoch=1:n_epochs
        print("Epoch: $epoch: ")
        cost = 0
        for (x, y) in batchview((X, Y))
            x, y = map(copy, (x, y))
            if cuda
                x, y = map(to_cuda, (x, y))
            end
            cost = partial_fit!(m, x, y)
        end
        println(cost)
    end
    return m
end

function predict_label(m, X)
    findmax(predict(m, X))[2]
end