text:
__precompile__()
module BlossomV

if isfile(joinpath(Pkg.dir("BlossomV"),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("BlossomV not properly installed. Please run Pkg.build(\"BlossomV\")")
end

export
    PerfectMatchingCtx,
    Matching, solve, add_edge, get_match, get_all_matches

const CostT = Union{Int32, Float64}

type PerfectMatchingCtx{T<:CostT}
    ptr::Ptr{Void}
    function (::Type{PerfectMatchingCtx{T}}){T<:CostT}(ptr::Ptr{Void})
        self = new{T}(ptr)
        finalizer(self, destroy)
        verbose(self, false)
        self
    end
end

dense_num_edges(node_num) =  node_num*(node_num-1)÷2

Matching(node_num::Integer) = Matching(Int32, node_num)
Matching(node_num::Integer, edge_num_max::Integer) = Matching(Int32, node_num, edge_num_max)

Matching{T}(::Type{T}, node_num::Integer) = Matching(T, node_num, dense_num_edges(node_num))

function Matching{T}(::Type{T}, node_num::Integer, edge_num_max::Integer)
    if T <: Integer
        return Matching(Int32, Int32(node_num), Int32(edge_num_max))
    elseif T <: AbstractFloat
        return Matching(Float64, Int32(node_num), Int32(edge_num_max))
    end
end

function add_edge{T}(matching::PerfectMatchingCtx{T}, first_node::Integer, second_node::Integer, cost::Number)
	add_edge(matching, Int32(first_node), Int32(second_node), T(cost))
end

get_match(matching::PerfectMatchingCtx, node::Integer) = get_match(matching, Int32(node))

function get_all_matches(m::PerfectMatchingCtx, n_nodes::Integer)
    ret = Matrix{Int32}(2, n_nodes÷2)
    assigned = falses(n_nodes)
    kk = 1
    for ii in 0:n_nodes-1
        assigned[ii+1] && continue
        jj = get_match(m, ii)
        @assert(!assigned[jj+1])
        assigned[ii+1] = true
        assigned[jj+1] = true

        ret[1,kk] = ii
        ret[2,kk] = jj
        kk+=1
    end
    @assert(kk-1 == n_nodes÷2)
    return ret
end

for (T, _jl_blossom5) in  ((Int32, :_jl_blossom5int32),  (Float64, :_jl_blossom5float64))
@eval begin

function destroy(matching::PerfectMatchingCtx{$T})
    if matching.ptr == C_NULL
        return
    end
    ccall((:matching_destruct, $_jl_blossom5), Void, (Ptr{Void},), matching.ptr)
    matching.ptr = C_NULL
    nothing
end

function Matching(::Type{$T}, node_num::Int32, edge_num_max::Int32)
    assert(node_num % 2 == 0)
    ptr = ccall((:matching_construct, $_jl_blossom5), Ptr{Void}, (Int32, Int32), node_num, edge_num_max)
    return PerfectMatchingCtx{$T}(ptr)
end

function add_edge(matching::PerfectMatchingCtx{$T}, first_node::Int32, second_node::Int32, cost::$T)
    first_node != second_node || error("Can not have an edge between $(first_node) and itself.")
    first_node >= 0  || error("first_node less than zero (value: $(first_node)). Indexes are zero-based.")
    second_node >= 0  || error("second_node less than zero (value: $(second_node)). Indexes are zero-based.")
    cost >= 0  || error("Cost must be positive. Edge between $(first_node) and $(second_node) has cost $cost.")
    ccall((:matching_add_edge, $_jl_blossom5), Int32, (Ptr{Void}, Int32, Int32, $T), matching.ptr, first_node, second_node, cost)
end

function solve(matching::PerfectMatchingCtx{$T})
    ccall((:matching_solve, $_jl_blossom5), Void, (Ptr{Void},), matching.ptr)
    nothing
end

function get_match(matching::PerfectMatchingCtx{$T}, node::Int32)
    ccall((:matching_get_match, $_jl_blossom5), Int32, (Ptr{Void}, Int32), matching.ptr, node)
end

function verbose(matching::PerfectMatchingCtx{$T}, verbose::Bool)
    ccall((:matching_verbose, $_jl_blossom5), Void, (Ptr{Void}, Bool), matching.ptr, verbose)
end

end
end

end # module


text:
# for parameters
function parameters()
    g = 10.; rho = 1000.; mu = 1.3076e-3
    return [g rho mu]
end

ext: jl


text:
# reduction

#################################################
#
#    Folding
#
#################################################

macro compose_foldfuns(AN, foldlf, foldrf)
    # argument preparation

    _foldlf = symbol("_$(foldlf)")
    _foldrf = symbol("_$(foldrf)")

    h = codegen_helper(AN)
    aparams = h.contiguous_aparams
    args = h.args
    ti = h.term(:i)
    t1 = h.term(1)
    tn = h.term(:n)

    # code skeletons

    quote
        # foldl & foldr

        global $_foldlf
        function $(_foldlf)(ifirst::Int, ilast::Int, op::Functor{2}, s::Number, $(aparams...))
            i = ifirst
            while i <= ilast
                @inbounds vi = $(ti)
                s = evaluate(op, s, vi)
                i += 1
            end
            return s
        end

        global $foldlf
        $(foldlf)(op::Functor{2}, s::Number, $(aparams...)) = $(_foldlf)(1, $(h.inputlen), op, s, $(args...))
        function $(foldlf)(op::Functor{2}, $(aparams...))
            n = $(h.inputlen)
            n > 0 || error("Empty argument not allowed.")
            s = $(t1)
            $(_foldlf)(2, n, op, s, $(args...))
        end

        global $_foldrf
        function $(_foldrf)(ifirst::Int, ilast::Int, op::Functor{2}, s::Number, $(aparams...))
            i = ilast
            while i >= ifirst
                @inbounds vi = $(ti)
                s = evaluate(op, vi, s)
                i -= 1
            end
            return s
        end

        global $foldrf
        $(foldrf)(op::Functor{2}, s::Number, $(aparams...)) = $(_foldrf)(1, $(h.inputlen), op, s, $(args...))
        function $(foldrf)(op::Functor{2}, $(aparams...))
            n = $(h.inputlen)
            n > 0 || error("Empty argument not allowed.")
            s = $(tn)
            $(_foldrf)(1, n-1, op, s, $(args...))
        end

    end
end

@compose_foldfuns 0 foldl foldr
@compose_foldfuns 1 foldl foldr
@compose_foldfuns 2 foldl foldr
@compose_foldfuns 3 foldl foldr
@compose_foldfuns (-2) foldl_fdiff foldr_fdiff

#################################################
#
#    Sum, Maximum, Minimum, and Mean
#
#################################################

# init value

reduceinit{T<:Number}(::Add, ::Type{T}) = zero(T)
reduceinit{T<:Unsigned}(::Add, ::Type{T}) = uint(0)
reduceinit{T<:Signed}(::Add, ::Type{T}) = 0
reduceinit(::Add, ::Type{Int64}) = @compat(Int64(0))
reduceinit(::Add, ::Type{Uint64}) = @compat(Uint64(0))
reduceinit(::Add, ::Type{Int128}) = @compat(Int128(0))
reduceinit(::Add, ::Type{Uint128}) = @compat(Uint128(0))
reduceinit(::Add, ::Type{Bool}) = 0

reduceinit{T<:Real}(::_Max, ::Type{T}) = typemin(T)
reduceinit{T<:Real}(::_Min, ::Type{T}) = typemax(T)
reduceinit{T<:Real}(::NonnegMax, ::Type{T}) = zero(T)

# functions

function reduce{T<:Number}(op::Functor{2}, a::ContiguousNumericArray{T})
    n = length(a)
    n > 0 ? saccum(op, n, parent(a), offset(a) + 1) : reduceinit(op, T)
end

function mapreduce(fun::Functor{1}, op::Functor{2}, a::ContiguousNumericArray)
    n = length(a)
    n > 0 ? saccum(op, n, fun, parent(a), offset(a) + 1) : reduceinit(op, result_type(fun, eltype(a)))
end

function mapreduce(fun::Functor{2}, op::Functor{2}, a::ContiguousArrOrNum, b::ContiguousArrOrNum)
    n = maplength(a, b);
    n > 0 ? saccum(op, n, fun, parent(a), offset(a) + 1, parent(b), offset(b) + 1) :
            reduceinit(op, result_type(fun, eltype(a), eltype(b)))
end

function mapreduce(fun::Functor{3}, op::Functor{2}, a::ContiguousArrOrNum, b::ContiguousArrOrNum, c::ContiguousArrOrNum)
    n = maplength(a, b, c);
    n > 0 ? saccum(op, n, fun, parent(a), offset(a) + 1, parent(b), offset(b) + 1, parent(c), offset(c) + 1) :
            reduceinit(op, result_type(fun, eltype(a), eltype(b), eltype(c)))
end

function mapreduce_fdiff(fun::Functor{1}, op::Functor{2}, a::ContiguousArrOrNum, b::ContiguousArrOrNum)
    n = maplength(a, b);
    n > 0 ? saccum_fdiff(op, n, fun, parent(a), offset(a) + 1, parent(b), offset(b) + 1) :
            reduceinit(op, result_type(fun, promote_type(eltype(a), eltype(b))))
end

macro compose_mapreduce_funs(rf, rfdiff, OT)
    quote
        global $(rf)
        $(rf)(fun::Functor{1}, a::ContiguousNumericArray) = mapreduce(fun, $(OT)(), a)
        $(rf)(fun::Functor{2}, a::ContiguousArrOrNum, b::ContiguousArrOrNum) =
            mapreduce(fun, $(OT)(), a, b)
        $(rf)(fun::Functor{3}, a::ContiguousArrOrNum, b::ContiguousArrOrNum, c::ContiguousArrOrNum) =
            mapreduce(fun, $(OT)(), a, b, c)

        global $(rfdiff)
        $(rfdiff)(fun::Functor{1}, a::ContiguousArrOrNum, b::ContiguousArrOrNum) =
            mapreduce_fdiff(fun, $(OT)(), a, b)
    end
end

@compose_mapreduce_funs sum sumfdiff Add
@compose_mapreduce_funs maximum maxfdiff _Max
@compose_mapreduce_funs minimum minfdiff _Min
@compose_mapreduce_funs nonneg_maximum nonneg_maxfdiff NonnegMax

mean(fun::Functor{1}, a::ContiguousNumericArray) = sum(fun, a) / length(a)
mean(fun::Functor{2}, a::ContiguousArrOrNum, b::ContiguousArrOrNum) =
    sum(fun, a, b) / maplength(a, b)
mean(fun::Functor{3}, a::ContiguousArrOrNum, b::ContiguousArrOrNum, c::ContiguousArrOrNum) =
    sum(fun, a, b, c) / maplength(a, b, c)
meanfdiff(fun::Functor{1}, a::ContiguousArrOrNum, b::ContiguousArrOrNum) =
    sumfdiff(fun, a, b) / maplength(a, b)

#################################################
#
#    Derived functions
#
#################################################

sumabs{T<:Number}(a::ContiguousArray{T}) = sum(AbsFun(), a)
maxabs{T<:Number}(a::ContiguousArray{T}) = nonneg_maximum(AbsFun(), a)
minabs{T<:Number}(a::ContiguousArray{T}) = minimum(AbsFun(), a)
meanabs{T<:Number}(a::ContiguousArray{T}) = mean(AbsFun(), a)

sumsq{T<:Real}(a::ContiguousArray{T}) = sum(Abs2Fun(), a)
meansq{T<:Real}(a::ContiguousArray{T}) = mean(Abs2Fun(), a)
dot{T<:Real}(a::ContiguousArray{T}, b::ContiguousArray{T}) = sum(Multiply(), a, b)

sumabsdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = sumfdiff(AbsFun(), a, b)
maxabsdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = nonneg_maxfdiff(AbsFun(), a, b)
minabsdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = minfdiff(AbsFun(), a, b)
meanabsdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = meanfdiff(AbsFun(), a, b)

sumsqdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = sumfdiff(Abs2Fun(), a, b)
maxsqdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = nonneg_maxfdiff(Abs2Fun(), a, b)
minsqdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = minfdiff(Abs2Fun(), a, b)
meansqdiff(a::ContiguousArrOrNum, b::ContiguousArrOrNum) = meanfdiff(Abs2Fun(), a, b)

sumxlogx{T<:Real}(a::ContiguousArray{T}) = sum(XlogxFun(), a)
sumxlogy{T<:Real}(a::ContiguousArray{T}, b::ContiguousArray{T}) = sum(XlogyFun(), a, b)
entropy{T<:Real}(a::ContiguousArray{T}) = -sumxlogx(a)


text:
using ImageInpainting
using ColorTypes
using TestImages
using Plots, VisualRegressionTests
using Test, Pkg

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

# helper functions
function plot_before_after(before, mask, after)
  gr(size=(800,400),yflip=true,colorbar=false,ticks=false)
  B = copy(before); B[mask] .= NaN
  p1 = heatmap(B,title="before")
  p2 = heatmap(after,title="after")
  plot(p1, p2, aspect_ratio=:equal)
end

# list of tests
testfiles = [
  "criminisi.jl"
]

@testset "ImageInpainting.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
# ----------
#  camelize
# ----------

@test camelize("innerHTML") == "innerHTML"
@test camelize("action_name") == "actionName"
@test camelize("css-class-name") == "cssClassName"
@test camelize("my favorite items") == "myFavoriteItems"
@test camelize("My Favorite Items") == "myFavoriteItems"
@test camelize("private-docs/owner-invoice") == "privateDocs/ownerInvoice"

# ------------
#  capitalize
# ------------

@test capitalize("innerHTML") == "InnerHTML"
@test capitalize("action_name") == "Action_name"
@test capitalize("css-class-name") == "Css-class-name"
@test capitalize("my favorite items") == "My favorite items"
@test capitalize("privateDocs/ownerInvoice") == "PrivateDocs/OwnerInvoice"

# ----------
#  classify
# ----------

@test classify("innerHTML") == "InnerHTML"
@test classify("action_name") == "ActionName"
@test classify("css-class-name") == "CssClassName"
@test classify("my favorite items") == "MyFavoriteItems"
@test classify("private-docs/owner-invoice") == "PrivateDocs/OwnerInvoice"

# -----------
#  dasherize
# -----------
@test dasherize("innerHTML") == "inner-html"
@test dasherize("action_name") == "action-name"
@test dasherize("css-class-name") == "css-class-name"
@test dasherize("my favorite items") == "my-favorite-items"
@test dasherize("privateDocs/ownerInvoice") == "private-docs/owner-invoice"

# ------------
#  decamelize
# ------------
@test decamelize("innerHTML") == "inner_html"
@test decamelize("action_name") == "action_name"
@test decamelize("css-class-name") == "css-class-name"
@test decamelize("my favorite items") == "my favorite items"

# ------------
#  underscore
# ------------

@test underscore("innerHTML") == "inner_html"
@test underscore("action_name") == "action_name"
@test underscore("css-class-name") == "css_class_name"
@test underscore("my favorite items") == "my_favorite_items"
@test underscore("privateDocs/ownerInvoice") == "private_docs/owner_invoice"


text:
"""
    subsampling_test(folder_name::String, subsampling_method::Function, args...)

Run a series of tests to check whether the subsampling functions in subsampling.jl work.
"""
function subsampling_test(folder_name::String, subsampling_method::Function, args...)

    # Load data
    Y = read_test_input("./input/subsampling/data");

    # Copy Y and subsample
    Y_copy = deepcopy(Y);
    args_copy = deepcopy(args);

    # Run `subsampling_method` with default random seed
    output = subsampling_method(Y, args...);

    # Load benchmark output
    benchmark = read_test_input("./input/subsampling/$(folder_name)/output_chunk1");
    benchmark_size = length(readdir("./input/subsampling/$(folder_name)"));

    for i=2:size(output,3)
        benchmark = cat(dims=3, benchmark, read_test_input("./input/subsampling/$(folder_name)/output_chunk$(i)"));
    end

    # Run tests
    @test Y_copy == Y;
    @test args_copy == args;
    @test size(output,3) == benchmark_size;
    @test sum(output .=== benchmark) == prod(size(output));
end

@testset "block jackknife" begin
    subsampling_test("block_jackknife", block_jackknife, 0.2)
end

@testset "artificial jackknife" begin
    subsampling_test("artificial_jackknife", artificial_jackknife, 0.2, 100)
end

@testset "moving block bootstrap" begin
    subsampling_test("moving_block_bootstrap", moving_block_bootstrap, 0.2, 100)
end

@testset "stationary block bootstrap" begin
    subsampling_test("stationary_block_bootstrap", stationary_block_bootstrap, 0.2, 100)
end


text:
using FlightSims
using Test
using Plots
gr()
using Transducers
using FSimZoo
using UnPack
using Convex, Mosek, MosekTools

function main()
    env = TwoDimensionalNonlinearDTSystem()
    x10, x20 = 0.5, 0.5
    x0 = State(env)(x10, x20)
    tf = 10
    simulator = Simulator(x0, apply_inputs(Dynamics!(env); u=(x, p, t) -> FSimZoo.OptimalControl(env)(x));
                          tf=tf,
                          Problem=:Discrete,
                         )
    df = solve(simulator)
    ts = df.time
    states = df.sol |> Map(datum -> datum.state) |> collect
    inputs = df.sol |> Map(datum -> datum.input) |> collect
    xs = states
    us = inputs
    fig_x = plot(ts, hcat(xs...)';
                 st=:scatter,
                 label=["x1" "x2"],
                 color=[:blue :red],
                 ylim=(-1, 1),
                )
    fig_u = plot(ts, hcat(us...)';
                 st=:scatter,
                 label="u",
                 ylim=(-1, 1),
                )
    fig_traj = plot(fig_x, fig_u; layout=(2, 1))
    # value functions
    x1s = -1:0.1:1
    x2s = x1s
    us = x1s
    Qs_x_fixed = us |> Map(u -> FSimZoo.OptimalQValue(env)(x0, u)) |> collect
    Q_u_fixed_func(x1, x2) = FSimZoo.OptimalQValue(env)(State(env)(x1, x2), zeros(1))
    Qs_u_and_x1_fixed = x2s |> Map(x2 -> Q_u_fixed_func(0.0, x2)) |> collect
    fig_Q = plot(x1s, x2s, Q_u_fixed_func;
                 st=:surface,
                 xlabel="x1",
                 ylabel="x2",
                 # zlim=(-1e6, 1e6),
                )
    # fig_V_true = plot(
    #                   x1s, x2s,
    #                   (x1, x2) -> FSimZoo.OptimalValue(env)(State(env)(x1, x2));
    #                   st=:surface,
    #                   zlim=(0, 5),
    #                  )
    function min_Q_numerical(x)
        u = Convex.Variable(1)
        dx = copy(x)
        @unpack x1, x2 = x
        @unpack c, d = env
        x1_next = c*(x1+x2)
        x2_next = c*x2 + u
        # FSimZoo.OptimalValue(env)(State(env)(x1_next, x2_next))
        V_next = d*FSimZoo.V1(env)(x1_next) + square(x2_next)
        a = FSimZoo.CubicSolution(env)(x2)
        r = (
             (1/4)*square(square(u))
             + (3/4)*a^4 + a^2
             - d*FSimZoo.V1(env)(x1_next) + d*FSimZoo.V1(env)(x1)
             + (1-c^2)*x2^2
            )  # for V(x) = d*sqrt(abs(x1)) + x2^2
        problem = minimize(r + V_next)
        solve!(problem, () -> Mosek.Optimizer(); silent_solver=true)
        problem.optval
    end
    function min_Q_exact(x)
        u = FSimZoo.OptimalControl(env)(x)
        optval = FSimZoo.OptimalQValue(env)(x, u)
    end
    fig_min_Q_exact_minus_V = plot(
                             x1s, x2s,
                             (x1, x2) -> min_Q_exact(State(env)(x1, x2)) - FSimZoo.OptimalValue(env)(State(env)(x1, x2));
                             st=:surface,
                             zlim=(-5, 5),
                            )
    fig_min_Q_numerical_minus_V = plot(
                             x1s, x2s,
                             (x1, x2) -> min_Q_numerical(State(env)(x1, x2)) - FSimZoo.OptimalValue(env)(State(env)(x1, x2));
                             st=:surface,
                             zlim=(-5, 5),
                            )
    fig_V = plot(
                 fig_min_Q_numerical_minus_V,
                 fig_min_Q_exact_minus_V;
                 layout=(1, 2),
                )
    display(fig_traj)
    # display(fig_Q)
    # display(fig_V)
end

@testset "two_dimensional_nonlinear_dt_system" begin
    main()
end


text:
__precompile__()

module DynProg

# include("rodcut.jl")
# include("bitstream.jl")
include("huffman.jl")
include("lcs.jl")

end


text:
immutable HermiteInterpolation{T} <: Interpolator{T}

    x₀::T
    x₁::T
    a₀::T
    a₁::T
    b₀::T
    b₁::T

    num::Vector{T}

    function HermiteInterpolation(x₀, x₁, d)
        a₀ = 2/(x₀-x₁)
        a₁ = 2/(x₁-x₀)

        b₀ = one(T)/(x₀-x₁)^2
        b₁ = one(T)/(x₁-x₀)^2

        new(x₀, x₁, a₀, a₁, b₀, b₁, zeros(T, d))
    end
end

function HermiteInterpolation{T}(x₀::T, x₁::T, d::Int)
    HermiteInterpolation{T}(x₀, x₁, d)
end

function evaluate{T}(int::HermiteInterpolation{T}, y₀::Vector{T}, y₁::Vector{T}, f₀::Vector{T}, f₁::Vector{T}, x::T, y::Vector{T})
    local d₀::T
    local d₁::T
    local c₀::T
    local c₁::T
    local den::T

    # Interpolate y values at required locations
    if x == int.x₀
        simd_copy!(y₀, y)
    elseif x == int.x₁
        simd_copy!(y₁, y)
    else
        d₀ = one(T)/(x-int.x₀)
        c₁ = d₀*int.b₀
        c₀ = c₁*(d₀-int.a₀)

        simd_copy_scale!(c₀, y₀, int.num)
        simd_axpy!(c₁, f₀, int.num)
        den = c₀

        d₁ = one(T)/(x-int.x₁)
        c₁ = d₁*int.b₁
        c₀ = c₁*(d₁-int.a₁)

        simd_axpy!(c₀, y₁, int.num)
        simd_axpy!(c₁, f₁, int.num)
        den += c₀

        simd_copy_scale!(one(T)/den, int.num, y)
    end
end


text:
mutable struct planar_lidar_t <: LCMType
    utime::Int64
    nranges::Int32
    ranges::Vector{Float32}

    nintensities::Int32
    intensities::Vector{Float32}

    rad0::Float32
    radstep::Float32
end

@lcmtypesetup(planar_lidar_t,
    ranges => (nranges,),
    intensities => (nintensities,),
)


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

mutable struct IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects <: SwaggerModel
    nonResourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule} } # spec name: nonResourceRules
    resourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule} } # spec name: resourceRules
    subjects::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1Subject} } # spec name: subjects

    function IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects(;nonResourceRules=nothing, resourceRules=nothing, subjects=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects, Symbol("nonResourceRules"), nonResourceRules)
        setfield!(o, Symbol("nonResourceRules"), nonResourceRules)
        validate_property(IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects, Symbol("resourceRules"), resourceRules)
        setfield!(o, Symbol("resourceRules"), resourceRules)
        validate_property(IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects, Symbol("subjects"), subjects)
        setfield!(o, Symbol("subjects"), subjects)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects

const _property_map_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects = Dict{Symbol,Symbol}(Symbol("nonResourceRules")=>Symbol("nonResourceRules"), Symbol("resourceRules")=>Symbol("resourceRules"), Symbol("subjects")=>Symbol("subjects"))
const _property_types_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects = Dict{Symbol,String}(Symbol("nonResourceRules")=>"Vector{IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule}", Symbol("resourceRules")=>"Vector{IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule}", Symbol("subjects")=>"Vector{IoK8sApiFlowcontrolV1alpha1Subject}")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects)
    (getproperty(o, Symbol("subjects")) === nothing) && (return false)
    true
end

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


text:
module TestGrid

using Test
using MLJBase
using MLJTuning
# include("../test/models.jl")
# using .Models
using ..Models
import Random.seed!
seed!(1234)

x1 = rand(100);
x2 = rand(100);
x3 = rand(100)
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.2*rand(100);

mutable struct DummyModel <: Deterministic
    lambda::Float64
    metric::Float64
    kernel::Char
end

dummy_model = DummyModel(4, 9.5, 'k')

mutable struct SuperModel <: Deterministic
    K::Int64
    model1::DummyModel
    model2::DummyModel
end

dummy_model = DummyModel(1.2, 9.5, 'k')
super_model = SuperModel(4, dummy_model, deepcopy(dummy_model))

s = range(super_model, :(model1.kernel), values=['c', 'd'])
r1 = range(super_model, :(model1.lambda), lower=20, upper=31)
r2 = range(super_model, :K, lower=1, upper=11, scale=:log10)

@testset "setup, default_n" begin
    user_range = [r1, (r2, 3), s]

    # with method:
    tuning = Grid(resolution=2, shuffle=false)
    @test MLJTuning.default_n(tuning, user_range) == 12
    models1 =
        params.((MLJTuning.setup(tuning, super_model, user_range, 3)).models)
    tuning = Grid(resolution=2, rng=123)
    @test MLJTuning.default_n(tuning, user_range) == 12
    models1r =
        params.((MLJTuning.setup(tuning, super_model, user_range, 3)).models)

    # by hand:
    m1 = [(K = 1, model1 = (lambda = 20.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 1, model1 = (lambda = 31.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 3, model1 = (lambda = 20.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 3, model1 = (lambda = 31.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 20.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 31.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 1, model1 = (lambda = 20.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 1, model1 = (lambda = 31.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 3, model1 = (lambda = 20.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 3, model1 = (lambda = 31.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 20.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 31.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))]

    @test m1 == models1
    @test m1 != models1r
    @test Set(m1) == Set(models1r)

    # with method:
    tuning = Grid(goal=9, shuffle=false)
    @test MLJTuning.default_n(tuning, user_range) == 8
    models2 =
        params.((MLJTuning.setup(tuning, super_model, user_range, 3)).models)

    # by hand:
    m2 = [(K = 1, model1 = (lambda = 20.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 1, model1 = (lambda = 31.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 20.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 31.0, metric = 9.5, kernel = 'c'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 1, model1 = (lambda = 20.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 1, model1 = (lambda = 31.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 20.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))
          (K = 11, model1 = (lambda = 31.0, metric = 9.5, kernel = 'd'),
           model2 = (lambda = 1.2, metric = 9.5, kernel = 'k'))]

    @test models2 == m2

end

@testset "2-parameter tune, with nesting" begin

    sel = FeatureSelector()
    stand = UnivariateStandardizer()
    ridge = FooBarRegressor()
    composite = SimpleDeterministicCompositeModel(transformer=sel, model=ridge)

    features_ = range(composite, :(transformer.features),
                      values=[[:x1], [:x1, :x2], [:x2, :x3], [:x1, :x2, :x3]])
    lambda_ = range(composite, :(model.lambda),
                    lower=1e-6, upper=1e-1, scale=:log10)

    r = [features_, lambda_]

    holdout = Holdout(fraction_train=0.8)
    grid = Grid(resolution=10)

    tuned_model = TunedModel(model=composite, tuning=grid,
                             resampling=holdout, measure=rms,
                             range=r)

    MLJBase.info_dict(tuned_model)

    tuned = machine(tuned_model, X, y)

    fit!(tuned, verbosity=0)
    r = report(tuned)
    @test r.best_report isa NamedTuple{(:machines, :reports)}
    fit!(tuned, verbosity=0)
    rep = report(tuned)
    fp = fitted_params(tuned)
    @test fp.best_fitted_params isa NamedTuple{(:machines, :fitted_params)}
    b = fp.best_model
    @test b isa SimpleDeterministicCompositeModel

    measurements = map(x->x.measurement[1], last.(tuned.report.history))
    # should be all different:
    @test length(unique(measurements)) == length(measurements)

    @test length(b.transformer.features) == 3
    @test abs(b.model.lambda - 0.027825) < 1e-6

    # get the training error of the tuned_model:
    e = rms(y, predict(tuned, X))

    # check this error has same order of magnitude as best measurement
    # during tuning:
    e_training = tuned.report.best_result.measurement[1]
    ratio = e/e_training
    @test ratio < 10 && ratio > 0.1

    # test weights:
    tuned_model.weights = rand(length(y))
    fit!(tuned, verbosity=0)
    @test e_training != tuned.report.best_result.measurement[1]

    # test plotting part of report:
    @test r.plotting.parameter_names ==
        ["transformer.features", "model.lambda"]
    @test r.plotting.parameter_scales == [:none, :log10]
    @test r.plotting.measurements == measurements
    @test size(r.plotting.parameter_values) == (40, 2)

end

@testset "basic tuning with training weights" begin

    seed!(1234)
    N = 100
    X = (x = rand(3N), );
    y = categorical(rand("abc", 3N));
    model = KNNClassifier()
    r = range(model, :K, lower=2, upper=N)
    tuned_model = TunedModel(model=model,
                             tuning=Grid(), measure=BrierScore(),
                             resampling=Holdout(fraction_train=2/3),
                             range=r)

    # no weights:
    tuned = machine(tuned_model, X, y)
    fit!(tuned, verbosity=0)
    best1 = fitted_params(tuned).best_model
    posterior1 = average([predict(tuned, X)...])

    # uniform weights:
    tuned = machine(tuned_model, X, y, fill(1, 3N))
    fit!(tuned, verbosity=0)
    best2 = fitted_params(tuned).best_model
    posterior2 = average([predict(tuned, X)...])

    @test best1 == best2
    @test all([pdf(posterior1, c) ≈ pdf(posterior2, c) for c in levels(y)])

    # skewed weights:
    w = map(y) do η
        if η == 'a'
            return 2
        elseif η == 'b'
            return 4
        else
            return 1
        end
    end
    tuned = machine(tuned_model, X, y, w)
    fit!(tuned, verbosity=0)
    best3 = fitted_params(tuned).best_model
    posterior3 = average([predict(tuned, X)...])

    # different tuning outcome:
    @test best1.K != best3.K

    # "posterior" is skewed appropriately in weighted case:
    @test abs(pdf(posterior3, 'b')/(2*pdf(posterior3, 'a'))  - 1) < 0.15
    @test abs(pdf(posterior3, 'b')/(4*pdf(posterior3, 'c'))  - 1) < 0.15

end

# ## LEARNING CURVE

# @testset "learning curves" begin
#     atom = FooBarRegressor()
#     ensemble = EnsembleModel(atom=atom, n=50, rng=1)
#     mach = machine(ensemble, X, y)
#     r_lambda = range(ensemble, :(atom.lambda),
#                      lower=0.0001, upper=0.1, scale=:log10)
#     curve = MLJ.learning_curve!(mach; range=r_lambda)
#     atom.lambda=0.3
#     r_n = range(ensemble, :n, lower=10, upper=100)
#     curve2 = MLJ.learning_curve!(mach; range=r_n)
#     curve3 = learning_curve(ensemble, X, y; range=r_n)
#     @test curve2.measurements ≈ curve3.measurements
# end

end # module
true


text:
using Bokeh

plotfile("functions.html")

plot(sin) # this uses the default range of -10 to 10
plot(sin, 0:20) # supply the rante
plot(x -> x^3 - 2x^2 + 3x, -5, 5) # also supplies the range
# the three plots above won't be shown since we haven't set hold, see below

plot([sin, cos, tan, sinh, cosh, tanh], -1:1)
showplot()


text:
# Traits
struct HasE end
struct HasCN end
struct HasCNE end

has_reserves(o) =
  if typeof(catabolism_pars(o)) <: AbstractCatabolismCN
      HasCN()
  elseif typeof(catabolism_pars(o)) <: AbstractCatabolismCNE
      HasCNE()
  elseif typeof(catabolism_pars(o)) <: AbstractCatabolism
      HasE()
  else
      nothing
  end


text:
"""
A wrapper for the PyAMG solver.

~~~julia
 amgSolver{Tv,Ti}(ddmat::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, pcgIts=Int[], verbose=false)
~~~
"""
function AMGSolver(ddmat::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, pcgIts=Int[], verbose=false) where {Tv,Ti}

  amg = PyAMG.RugeStubenSolver(ddmat);
  M = PyAMG.aspreconditioner(amg);
  function F(b)
    return M \ b;
  end

    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

  f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, pcgIts=pcgIts_, verbose=verbose_) = pcg(ddmat, b, F, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, verbose=verbose)

  return f

end

"""
A wrapper for the PyAMG solver. In line with our other solvers, takes in an adjacency matrix.

~~~julia
 amgSolver{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, pcgIts=Int[], verbose=false)
~~~
"""
function AMGLapSolver(a::SparseMatrixCSC{Tv,Ti}; tol::Float64=1e-6, maxits=Inf, maxtime=Inf, pcgIts=Int[], verbose=false) where {Tv,Ti}

  la = lap(a)

  amg = PyAMG.RugeStubenSolver(la);
  M = PyAMG.aspreconditioner(amg);
  function F(b)
    return M \ b;
  end

    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

  f(b;tol=tol_,maxits=maxits_, maxtime=maxtime_, pcgIts=pcgIts_, verbose=verbose_) = pcg(la, b-mean(b), F, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, verbose=verbose)

  return f

end


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

mutable struct NetworkInterfaceDnsSettings <: SwaggerModel
    dnsServers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: dnsServers
    appliedDnsServers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: appliedDnsServers
    internalDnsNameLabel::Any # spec type: Union{ Nothing, String } # spec name: internalDnsNameLabel
    internalFqdn::Any # spec type: Union{ Nothing, String } # spec name: internalFqdn
    internalDomainNameSuffix::Any # spec type: Union{ Nothing, String } # spec name: internalDomainNameSuffix

    function NetworkInterfaceDnsSettings(;dnsServers=nothing, appliedDnsServers=nothing, internalDnsNameLabel=nothing, internalFqdn=nothing, internalDomainNameSuffix=nothing)
        o = new()
        validate_property(NetworkInterfaceDnsSettings, Symbol("dnsServers"), dnsServers)
        setfield!(o, Symbol("dnsServers"), dnsServers)
        validate_property(NetworkInterfaceDnsSettings, Symbol("appliedDnsServers"), appliedDnsServers)
        setfield!(o, Symbol("appliedDnsServers"), appliedDnsServers)
        validate_property(NetworkInterfaceDnsSettings, Symbol("internalDnsNameLabel"), internalDnsNameLabel)
        setfield!(o, Symbol("internalDnsNameLabel"), internalDnsNameLabel)
        validate_property(NetworkInterfaceDnsSettings, Symbol("internalFqdn"), internalFqdn)
        setfield!(o, Symbol("internalFqdn"), internalFqdn)
        validate_property(NetworkInterfaceDnsSettings, Symbol("internalDomainNameSuffix"), internalDomainNameSuffix)
        setfield!(o, Symbol("internalDomainNameSuffix"), internalDomainNameSuffix)
        o
    end
end # type NetworkInterfaceDnsSettings

const _property_map_NetworkInterfaceDnsSettings = Dict{Symbol,Symbol}(Symbol("dnsServers")=>Symbol("dnsServers"), Symbol("appliedDnsServers")=>Symbol("appliedDnsServers"), Symbol("internalDnsNameLabel")=>Symbol("internalDnsNameLabel"), Symbol("internalFqdn")=>Symbol("internalFqdn"), Symbol("internalDomainNameSuffix")=>Symbol("internalDomainNameSuffix"))
const _property_types_NetworkInterfaceDnsSettings = Dict{Symbol,String}(Symbol("dnsServers")=>"Vector{String}", Symbol("appliedDnsServers")=>"Vector{String}", Symbol("internalDnsNameLabel")=>"String", Symbol("internalFqdn")=>"String", Symbol("internalDomainNameSuffix")=>"String")
Base.propertynames(::Type{ NetworkInterfaceDnsSettings }) = collect(keys(_property_map_NetworkInterfaceDnsSettings))
Swagger.property_type(::Type{ NetworkInterfaceDnsSettings }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_NetworkInterfaceDnsSettings[name]))}
Swagger.field_name(::Type{ NetworkInterfaceDnsSettings }, property_name::Symbol) =  _property_map_NetworkInterfaceDnsSettings[property_name]

function check_required(o::NetworkInterfaceDnsSettings)
    true
end

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


text:
# This file comtains some useful functions to import logfiles for use in either julia or matlab.
# In julia, the data object from orcalog2mat can be used together with the function getData
# In matlab, the function csv2mat can be used to convert the textfile to a .mat-file which is very fast to read
using MAT

"""`findData(pattern::AbstractString, data)`"""
function findData(pattern::AbstractString, data)
    indexes = falses(length(keys(data)))
    for (i,key) in enumerate(keys(data))
        if match(Regex(pattern), key) != nothing
            indexes[i] = true
        end
    end
    return indexes
end

"""
`data = orcalog2mat(pathopen, pathsave)`
`pathopen` is the full path to the file
`pathsave` is the full path to the saved file, should typically end with .mat
the returned `data` object is a `Dict{ByteString,Any}`
"""
function orcalog2mat(pathopen, pathsave; separator = ',')
    csv2mat(pathopen, pathsave, writeNames = false, separator = separator)
    data = MAT.matread(pathsave)
    return data
end

"""
`data = readmat(pathopen)`
`pathopen` is the full path to the file
the returned `data` object is a `Dict{ByteString,Any}`
"""
function readmat(pathopen)
    data = MAT.matread(pathopen)
    return data
end

"""
`array = getData(pattern, data, ds=1; removeNaN = false)`
pattern is a string representation of a regex used to match the data
e.g. \"posRawAbs\"
"""
function getData(pattern, data, ds=1; removeNaN = false)
    indexes = findData(pattern, data)
    if sum(indexes) < 1
        warn("No data found using the string $pattern")
        return
    end
    names = collect(keys(data))[indexes]
    function lt(x,y) # String comparison function, this soring places 2 in front of 10
        if length(x) < length(y)
            return true
        elseif length(x) > length(y)
            return false
        else
            return x < y
        end
    end
    sort!(names,lt=lt)
    retmat = Matrix{Float64}(ceil(Int,length(data[names[1]])/ds),length(names))
    for (j,key) in enumerate(names)
        for (i,k)= enumerate(1:ds:length(data[key]))
            retmat[i,j] = data[key][k]
        end
    end

    if removeNaN
        retmat[any(isnan(retmat),2),:] = []
    end
    return retmat
end


text:
using MPIClusterManagers
using Distributed

nprocs = length(readlines(ENV["PBS_NODEFILE"]))
manager = MPIManager(np=nprocs, master_tcp_interface=gethostname())

addprocs(manager)
# addprocs(manager, exename=`$(Sys.BINDIR)/julia -J$(pwd())/PWS.jl/GMcmcSysimage.so`)

using Logging
@info "Successfully launched MPI processes"

include(ARGS[2])


text:
module LogicalKoans

function not(x)
    !x
end

function equals(x, y)
    x == y
end

function not_equals(x, y)
    x != y
end

function less_than(x, y)
    x < y
end

function greater_than(x, y)
    x > y
end

function less_than_or_equal_to(x, y)
    x <=y
end

function greater_than_or_equal_to(x, y)
    x >= y
end

end


text:
using ApproxFun

 # Semicircle law
x=Fun(identity,-2..2)
w=sqrt(4-x^2)/(2π)

lanczos(w,5)

# Marchenco–Pastur law

r = .5
lmax = (1+sqrt(r))^2
lmin = (1-sqrt(r))^2

x= Fun(identity,lmin..lmax)

w=sqrt((lmax-x)*(x-lmin))/(π*x)

lanczos(w,5)

# Wachter law

a = 5; b= 10;
c = sqrt(a/(a+b)*(1-1/(a+b)))
d = sqrt(1/(a+b)*(1-a/(a+b)))

lmax = (c+d)^2
lmin = (c-d)^2

x = Fun(identity,lmin..lmax)

w = (a+b)*sqrt((x-lmin).*(lmax-x))/(2π*x*(1-x))

lanczos(w,5)


text:
module ARFIMA

using Random, LinearAlgebra, StaticArrays
export arfima, SVector, @SVector, arma

"""
    arfima([rng,] N, σ, d, φ=nothing, θ=nothing) -> Xₜ
Create a stochastic timeseries of length `N` that follows the ARFIMA
process, or any of its subclasses, like e.g. ARMA, AR, ARIMA, etc., see below.
`σ` is the standard deviation of the white noise used to generate the
process. The first optional argument is an `AbstractRNG`, a random
number generator to establish reproducibility.

The generating equation for `Xₜ` is:
```math
\\left( 1 - \\sum_{i=1}^p \\phi_i B^i \\right)
\\left( 1-B \\right)^d X_t
=
\\left( 1 + \\sum_{i=1}^q \\theta_i B^i \\right) \\varepsilon_t
```
with ``B`` the backshift operator and ``\\varepsilon_t`` white noise.

This equation encapsulates all possible variants of ARFIMA and Julia's
multiple dispatch system decides which will be the simulated variant,
based on the types of `d, φ, θ`.

## Variants
The ARFIMA parameters are (p, d, q) with `p = length(φ)` and `q = length(θ)`,
with `p, q` describing the autoregressive or moving average "orders" while
`d` is the differencing "order".
Both `φ, θ` can be of two types: `Nothing` or `SVector`. If they are `Nothing`
the corresponding components of autoregressive (φ) and moving average (θ)
are not done. Otherwise, the static vectors simply contain their values.

If `d` is `Nothing`, then the differencing (integrated)
part is not done and the process is in fact AR/MA/ARMA.
If `d` is of type `Int`, then the simulated process is in fact ARIMA,
while if `d` is `AbstractFloat` then the process is AR**F**IMA.
In the last case it must hold that `d ∈ (-0.5, 0.5)`.
If all `d, φ, θ` are `nothing`, white noise is returned.

The function `arma(N, σ, φ, θ = nothing)` is provided for convienience.

## Examples
```julia
N, σ = 10_000, 0.5
arfima(N, σ, 0.4)                                   # ARFIMA(0,d,0)
arfima(N, σ, 0.4, SVector(0.8))                     # ARFIMA(1,d,0)
arfima(N, σ, 1, SVector(0.8))                       # ARIMA(1,d,0)
arfima(N, σ, 1, SVector(0.8), SVector(1.2))         # ARIMA(1,d,1)
arfima(N, σ, 0.4, SVector(0.8), SVector(1.2))       # ARFIMA(1,d,1)
arfima(N, σ, nothing, SVector(0.8))                 # AR(1)
arfima(N, σ, nothing,  SVector(0.8), SVector(1.2))  # ARMA(1,1)
```
"""
arfima(N::Int, args...) = arfima(Random.GLOBAL_RNG, N, args...)
arfima(rng::AbstractRNG, N, σ, d) = arfima(rng, N, σ, d, nothing, nothing)
arfima(rng::AbstractRNG, N, σ, d, φ) = arfima(rng, N, σ, d, φ, nothing)

arfima(rng, N, σ, ::Nothing, ::Nothing, θ) = generate_noise(rng, N, σ, θ) # MA
arma(N::Int, args...) = arma(Random.GLOBAL_RNG, N, args...)
arma(rng, N, σ, φ, θ = nothing) = arfima(rng, N, σ, nothing, φ, θ)

function arfima(rng, N, σ, d, φ::SVector{P}, θ) where {P} # AR(F)IMA
    L = estimate_past(φ)
    Z = arfima(rng, L+N, σ, d, nothing, θ) # Generate FIMA or IMA
    X = autoregressive(N, Z, φ) # autoregress on the result
end

function arfima(rng, N, σ, d::AbstractFloat, φ::Nothing, θ) # FIMA
    @assert -0.5 ≤ d ≤ 0.5 "For (AR)FIMA, it must be d ∈ (-0.5, 0.5)"
    M = 2N # infinite summation becomes summation over 2N
    noise = generate_noise(rng, M+N-1, σ, θ)
    # coefficients of truncated sum in fractional derivative
    ψ = zeros(M); ψ[1] = 1.0
    for k in 0:M-2
        @inbounds @fastmath ψ[k+2] = ψ[k+1] * (d+k)/(k+1)
    end
    # Here we create the actual arfima process Xₜ
    X = zeros(N)
    for k in 1:N
        @inbounds X[k] = dot(ψ, view(noise, k:k+M-1))
    end
    return X
end

@inline function arfima(rng, N, σ, d::Int, φ::Nothing, θ) # IMA
    @assert d>0
    M = N + 2d
    noise = generate_noise(rng, M, σ, θ)
    differencing = SVector{d}([(-1)^(k+1) * binomial(d, k) for k in 1:d]...)
    X = zeros(N+d)
    _hotloop1!(X, noise, differencing, d, N)
    return deleteat!(X, 1:d)
end
# hotloop1: so that dispach doesn't happen dynamically on SVector{d}:
function _hotloop1!(X, noise, differencing, d, N)
    for i in d+1:N+d
        @inbounds X[i] = bdp(differencing, X, i) + noise[d+i]
    end
end

function arfima(rng, N, σ, d::Nothing, φ::SVector{P}, θ) where {P} # AR(MA)
     noise = generate_noise(rng, N + P, σ, θ)
     L = estimate_past(φ)
     Z = generate_noise(rng, N + L, σ, θ) # white noise
     X = autoregressive(N, Z, φ)
end

function generate_noise(rng, N, σ, θ::Nothing)
    noise = randn(rng, N)
    if σ == 1.0
        return noise
    else
        return σ .* noise
    end
end

function generate_noise(rng, N, σ, θ::SVector{Q}) where {Q} # MA
    @warn "There is a critical in the Moving Average code. Please help us fix it. "*
    "https://github.com/JuliaDynamics/ARFIMA.jl/issues/3"
    ε = generate_noise(rng, N+Q, σ, nothing)
    noise = zeros(N)
    θ = -θ # this change is necessary due to the defining equation
    # simply now do the average process
    @inbounds for i in 1:N
        noise[i] = bdp(θ, ε, i+Q) + ε[i]
    end
    return noise
end

"""
    bdp(φ::SVector, X::AbstractVector, t)
Perform the backshift dot product between `φ` and `X`, with starting index `t`:
``\\sum_i \\phi_i X_{t-i}``.
"""
@generated function bdp(φ::SVector{P}, X::AbstractVector, t) where {P}
    exprs = [:(φ[$i]*X[t-$i]) for i in 1:P]
    quote
        Base.@_inline_meta
        @inbounds @fastmath +($(exprs...))
    end
end

"Estimate how long into the past to go for accurate AR process."
estimate_past(φ::SVector{P}) where {P} =
max(P+1, ceil(Int, log(0.001)/log(maximum(abs, φ))))

"""
    autoregressive(N, Z, φ::SVector{P}) -> X
Generate an autoregressive process of length `N` based on input noise term `Z`.
This is used in both ARFIMA and ARMA.
The noise term must have at least `N+P` elements.
"""
function autoregressive(N, Z, φ::SVector{P}) where {P}
    L = length(Z) - N; @assert L > P
    tmp = zeros(P)
    # Generate correct inital condition: the first P values of X
    for i = 1:L-P;
        y = bdp(φ, tmp, P+1) + Z[i];
        tmp[1:end-1] .= tmp[2:end] # shift values and add the new value
        tmp[end] = y
    end
    # X0 is now the "correct" X0, after L steps in advance
    X = zeros(N); X[1:P] .= tmp
    for i in (P+1):N
        @inbounds X[i] = bdp(φ, X, i) + Z[i+L]
    end
    return X
end

end # module


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_translated)" --startup-file=no %s 2>&1 | FileCheck %s

function branches(c)
    if c
        return c
    else
        return !c
    end
end
emit(branches, Bool)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.branches), Bool}, svec(), branches(c) in Main.Main at /{{.*}}/test/Codegen/translate/branches..jl:3, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.branches), Bool}"(%arg0: !jlir<"typeof(Main.branches)">, %arg1: !jlir.Bool) -> !jlir.Bool attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     "jlir.gotoifnot"(%arg1)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     "jlir.return"(%arg1) : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %1 = "jlir.call"(%0, %arg1) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.return"(%1) : (!jlir.Bool) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: stephenneuendorffer/brutus

max_issues_repo_name: stephenneuendorffer/brutus

max_forks_repo_name: stephenneuendorffer/brutus


text:
#structure to hold the solver parameters
@with_kw struct SolverParams{T <: Real, N <: Integer}
    n_iter_viscosity::N = 2;  @assert n_iter_viscosity ==2
    tol_picard::T = 1e-5
    maxiter_picard::N = 30
    tol_coarse::T = 1e-5
    maxiter_coarse::N = 1000
    levels::N = 3
    wavelet_threshold::T = 10.0
    nsmooth::N = 5
    smoother_omega::T = 1.0
    stencil_margin::N = 3
end


text:
module TriplotRecipes

using PlotUtils,RecipesBase,TriplotBase

export tricontour,tricontour!,tripcolor,tripcolor!,dgtripcolor,dgtripcolor!,trimesh,trimesh!

function append_with_nan!(a,b)
    append!(a,b)
    push!(a,NaN)
end

@recipe function f(contours::Vector{TriplotBase.Contour{T}}) where {T}
    color = get(plotattributes, :seriescolor, :auto)
    set_line_z = (color == :auto || plot_color(color) isa ColorGradient)
    for c=contours
        x = T[]
        y = T[]
        z = T[]
        for polyline=c.polylines
            append_with_nan!(x,first.(polyline))
            append_with_nan!(y,last.(polyline))
            append!(z,fill(c.level,length(polyline)))
        end
        @series begin
            label := nothing
            if set_line_z
                line_z := z
            end
            x,y
        end
    end
end

tricontour(x,y,z,t,levels;kw...) = RecipesBase.plot(TriplotBase.tricontour(x,y,z,t,levels);kw...)
tricontour!(x,y,z,t,levels;kw...) = RecipesBase.plot!(TriplotBase.tricontour(x,y,z,t,levels);kw...)

struct TriPseudocolor{X,Y,Z,T} x::X; y::Y; z::Z; t::T; end

@recipe function f(p::TriPseudocolor;px=512,py=512,ncolors=256)
    cmap = range(extrema(p.z)...,length=ncolors)
    x = range(extrema(p.x)...,length=px)
    y = range(extrema(p.y)...,length=py)
    z = TriplotBase.tripcolor(p.x,p.y,p.z,p.t,cmap;bg=NaN,px=px,py=py)
    seriestype := :heatmap
    x,y,z'
end

tripcolor(x,y,z,t;kw...) = RecipesBase.plot(TriPseudocolor(x,y,z,t);kw...)
tripcolor!(x,y,z,t;kw...) = RecipesBase.plot!(TriPseudocolor(x,y,z,t);kw...)

struct DGTriPseudocolor{X,Y,Z,T} x::X; y::Y; z::Z; t::T; end

@recipe function f(p::DGTriPseudocolor;px=512,py=512,ncolors=256)
    cmap = range(extrema(p.z)...,length=ncolors)
    x = range(extrema(p.x)...,length=px)
    y = range(extrema(p.y)...,length=py)
    z = TriplotBase.dgtripcolor(p.x,p.y,p.z,p.t,cmap;bg=NaN,px=px,py=py)
    seriestype := :heatmap
    x,y,z'
end

dgtripcolor(x,y,z,t;kw...) = RecipesBase.plot(DGTriPseudocolor(x,y,z,t);kw...)
dgtripcolor!(x,y,z,t;kw...) = RecipesBase.plot!(DGTriPseudocolor(x,y,z,t);kw...)

struct TriMesh{X,Y,T} x::X; y::Y; t::T; end

@recipe function f(m::TriMesh)
    x = Vector{eltype(m.x)}()
    y = Vector{eltype(m.y)}()
    for t=eachcol(m.t)
        append_with_nan!(x,[m.x[t];m.x[t[1]]])
        append_with_nan!(y,[m.y[t];m.y[t[1]]])
    end
    seriestype := :shape
    seriescolor --> RGB(0.7,1.0,0.8)
    label --> nothing
    x,y
end

trimesh(x,y,t;kw...) = RecipesBase.plot(TriMesh(x,y,t);kw...)
trimesh!(x,y,t;kw...) = RecipesBase.plot!(TriMesh(x,y,t);kw...)

end


text:
const extmap = Dict{String, Val}(
    ".dta" => Val(:dta),
    ".sav" => Val(:sav),
    ".por" => Val(:por),
    ".sas7bdat" => Val(:sas7bdat),
    ".xpt" => Val(:xport)
)

"""
    ColumnIndex

A type union for values accepted by [`readstat`](@ref) for selecting a column.
A column can be selected either with the column name as `Symbol` or `String`;
or with an integer index based on the position in a table.
See also [`ColumnSelector`](@ref).
"""
const ColumnIndex = Union{Symbol, String, Integer}

"""
    ColumnSelector

A type union for values accepted by [`readstat`](@ref)
for selecting either a single column or multiple columns.
The accepted values must be either a [`ColumnIndex`](@ref)
or a vector of [`ColumnIndex`](@ref).
"""
const ColumnSelector = Union{ColumnIndex, AbstractVector{<:Union{ColumnIndex}}}

function _parse_usecols(file, usecols::Union{Symbol, String})
    c = findfirst(x->x==Symbol(usecols), file.headers)
    c === nothing && throw(ArgumentError("column name $usecols is not found"))
    return (c,)
end

function _parse_usecols(file, usecols::AbstractVector{<:Union{Symbol, String}})
    lookup = Dict(n=>i for (i, n) in enumerate(file.headers))
    icols = Vector{Int}(undef, length(usecols))
    for (i, c) in enumerate(usecols)
        k = get(lookup, Symbol(c), 0)
        k == 0 && throw(ArgumentError("column name $c is not found"))
        icols[i] = k
    end
    return icols
end

function _parse_usecols(file, usecols::Integer)
    N = length(file.data)
    0 < usecols <= N || throw(ArgumentError("invalid column index $usecols"))
    return (usecols,)
end

function _parse_usecols(file, usecols::AbstractVector{<:Integer})
    N = length(file.data)
    for c in usecols
        0 < c <= N || throw(ArgumentError("invalid column index $c"))
    end
    return usecols
end

_selected(i::Int, n::Symbol, sel::Bool) = sel
_selected(i::Int, n::Symbol, sel::Integer) = i == sel
_selected(i::Int, n::Symbol, sel::Symbol) = n == sel
_selected(i::Int, n::Symbol, sel::AbstractVector{<:Integer}) = i in sel
_selected(i::Int, n::Symbol, sel::AbstractVector{Symbol}) = n in sel

function _to_array!(d::DataValueVector, a::Vector, imissing::Int, missingvalue)
    if imissing > 1
        @inbounds for i in 1:imissing-1
            a[i] = d.values[i]
        end
    end
    @inbounds for i in imissing:length(d)
        a[i] = d.isna[i] ? missingvalue : d.values[i]
    end
end

function _to_array(d::DataValueVector, missingvalue)
    N = length(d)
    imissing = findfirst(d.isna)
    if imissing === nothing
        return d.values, false
    else
        V = promote_type(typeof(d).parameters[1], typeof(missingvalue))
        a = Vector{V}(undef, N)
        _to_array!(d, a, imissing, missingvalue)
        return a, true
    end
end

"""
    readstat(filepath::AbstractString; kwargs...)

Return a [`ReadStatTable`](@ref) that collects data from a supported data file
located at `filepath`.

# Accepted File Extensions
- Stata: `.dta`.
- SAS: `.sas7bdat` and `.xpt`.
- SPSS: `.sav` and `por`.

# Keywords
- `usecols::Union{ColumnSelector, Nothing} = nothing`: only keep data from the specified columns (variables); keep all columns if `usecols=nothing`.
- `convert_datetime::Union{Bool, ColumnSelector} = true`: convert data from the specified columns to `Date` or `DateTime` if they are recorded in supported time formats; if specified as `true` (`false`), always (never) convert the data whenever possible.
- `apply_value_labels::Union{Bool, ColumnSelector} = true`: convert data from the specified columns to [`LabeledArray`](@ref) with their value labels; if specified as `true` (`false`), always (never) convert the data whenever possible.
- `missingvalue = missing`: value used to fill any missing value (should be `missing` unless in special circumstances).
"""
function readstat(filepath::AbstractString;
        usecols::Union{ColumnSelector, Nothing} = nothing,
        convert_datetime::Union{Bool, ColumnSelector} = true,
        apply_value_labels::Union{Bool, ColumnSelector} = true,
        missingvalue = missing)

    ext = lowercase(splitext(filepath)[2])
    filetype = get(extmap, ext, nothing)
    filetype === nothing && throw(ArgumentError("file extension $ext is not supported"))
    file = read_data_file(filepath, filetype)
    if usecols === nothing
        icols = 1:length(file.data)
        names = file.headers
    else
        icols = _parse_usecols(file, usecols)
        names = [file.headers[i] for i in icols]
    end
    convert_datetime isa AbstractVector{String} &&
        (convert_datetime = Symbol.(convert_datetime))
    # ReadStat.jl does not handle value labels for SAS at this moment
    (filetype == Val(:sas7bdat) || filetype == Val(:xport)) &&
        (apply_value_labels = false)
    apply_value_labels isa AbstractVector{String} &&
        (apply_value_labels = Symbol.(apply_value_labels))

    cols = Vector{AbstractVector}(undef, length(icols))
    varlabels = Dict{Symbol, String}()
    formats = Dict{Symbol, String}()
    val_label_keys = Dict{Symbol, String}()
    @inbounds for (i, c) in enumerate(icols)
        col, hasmissing = _to_array(file.data[c], missingvalue)
        n = names[i]
        if _selected(c, n, convert_datetime)
            format = file.formats[c]
            filetype == Val(:dta) && (format = first(format, 3))
            dtpara = get(dt_formats[filetype], format, nothing)
            if dtpara !== nothing
                if hasmissing
                    col = parse_datetime(col, dtpara..., missingvalue)
                else
                    col = parse_datetime(col, dtpara...)
                end
            end
        end
        if _selected(c, n, apply_value_labels)
            lblname = file.val_label_keys[c]
            if lblname != ""
                lbls = convert(Dict{eltype(col), String}, file.val_label_dict[lblname])
                col = LabeledArray(col, lbls)
            end
        end
        cols[i] = col
        varlabels[n] = file.labels[c]
        formats[n] = file.formats[c]
        val_label_keys[n] = file.val_label_keys[c]
    end
    meta = ReadStatMeta(varlabels, formats, val_label_keys,
        file.val_label_dict, file.filelabel, file.timestamp, ext)
    return ReadStatTable(cols, names, meta)
end


text:
"""
    struct ASTContext <: Any
Hold a pointer to a `clang::ASTContext` object.
"""
struct ASTContext
    ptr::CXASTContext
end

Base.unsafe_convert(::Type{CXASTContext}, x::ASTContext) = x.ptr
Base.cconvert(::Type{CXASTContext}, x::ASTContext) = x


text:
"""
A sub-module for a specific idea for how the basics of Breadth works.
"""

include("Types.jl")
include("Singletons.jl")
include("Proving.jl")
include("truths/__init__.jl")
include("Provers.jl")

P = symbol("P", Statements)
Q = symbol("Q", Statements)
R = symbol("R", Statements)

x = symbol("x", Statements)
A = symbol("A", Sets)
B = symbol("B", Sets)
C = symbol("C", Sets)


text:
using Distributions, Random
using Test, LinearAlgebra, PDMats

#  dimensions

n = 6
p = 3

#  parameters

v    = 10
M    = randn(n, p)
Σ    = rand(InverseWishart(n + 2, Matrix(1.0I, n, n)))
Ω    = rand(InverseWishart(p + 2, Matrix(1.0I, p, p)))
ΣF32 = Matrix{Float32}(Σ)
ΩBF  = Matrix{BigFloat}(Ω)
PDΣ  = PDMat(Σ)
PDΩ  = PDMat(Ω)
μ    = randn(n)
θ    = randn(p)
σ²   = 2.0rand()

#  distribution instances

D = MatrixTDist(v, M, Σ, Ω)                                   #  n x p
L = MatrixTDist(v, reshape(μ, n, 1), Σ, reshape([σ²], 1, 1))  #  n x 1
H = MatrixTDist(v, reshape(θ, 1, p), reshape([σ²], 1, 1), Ω)  #  1 x p
K = MatrixTDist(1000, M, 1000Σ, Ω)

l  = MvTDist(L)
h  = MvTDist(H)
MN = MatrixNormal(M, Σ, Ω)

#  draws

A = rand(D)

X = rand(L)
Y = rand(H)

x = vec(X)
y = vec(Y)

@testset "Check all MatrixTDist constructors" begin
    @test MatrixTDist(v, M, Σ,        PDΩ.chol) isa MatrixTDist
    @test MatrixTDist(v, M, Σ,        PDΩ)      isa MatrixTDist
    @test MatrixTDist(v, M, PDΣ.chol, Ω)        isa MatrixTDist
    @test MatrixTDist(v, M, PDΣ.chol, PDΩ.chol) isa MatrixTDist
    @test MatrixTDist(v, M, PDΣ.chol, PDΩ)      isa MatrixTDist
    @test MatrixTDist(v, M, PDΣ,      Ω)        isa MatrixTDist
    @test MatrixTDist(v, M, PDΣ,      PDΩ.chol) isa MatrixTDist
    @test MatrixTDist(v, M, PDΣ,      PDΩ)      isa MatrixTDist
end

@testset "MatrixTDist promotion during construction" begin
    R = MatrixTDist(v, M, ΣF32, ΩBF)
    @test partype(R) == BigFloat
end

@testset "MatrixTDist construction errors" begin
    @test_throws ArgumentError MatrixTDist(-1, M, Σ, Ω)
    @test_throws ArgumentError MatrixTDist(Inf, M, Σ, Ω)
    @test_throws ArgumentError MatrixTDist(v, M, Ω, Σ)
    @test_throws ArgumentError MatrixTDist(v, M, Ω, Ω)
    @test_throws ArgumentError MatrixTDist(v, M, Σ, Σ)
end

@testset "MatrixTDist params" begin
    vv, MM, PDΣΣ, PDΩΩ = params(D)
    @test v == vv
    @test M == MM
    @test Σ == PDΣΣ.mat
    @test Ω == PDΩΩ.mat
end

@testset "MatrixTDist size" begin
    @test size(D) == (n, p)
    @test size(L) == (n, 1)
    @test size(H) == (1, p)
end

@testset "MatrixTDist rank" begin
    @test rank(D) == min(n, p)
    @test rank(L) == 1
    @test rank(H) == 1
    @test rank(D) == rank(rand(D))
    @test rank(L) == rank(rand(L))
    @test rank(H) == rank(rand(H))
end

@testset "MatrixTDist insupport" begin
    @test insupport(D, rand(D))
    @test insupport(L, rand(L))
    @test insupport(H, rand(H))

    @test !insupport(D, Matrix(rand(D)'))
    @test !insupport(D, M + M * im)
    @test !insupport(L, rand(H))
    @test !insupport(H, rand(L))
end

@testset "MatrixTDist mean" begin
    @test mean(D) == M
    @test mean(L) == reshape(μ, n, 1)
    @test mean(H) == reshape(θ, 1, p)

    @test_throws ArgumentError mean( MatrixTDist(1, M, Σ, Ω) )
end

@testset "MatrixTDist mode" begin
    @test mode(D) == M
    @test mode(L) == reshape(μ, n, 1)
    @test mode(H) == reshape(θ, 1, p)
end

@testset "MatrixTDist cov and var" begin
    @test vec(var(D)) ≈ diag(cov(D))
    @test reshape(cov(D), n, p, n, p) ≈ cov(D, Val(false))
    @test cov(L) ≈ cov(l)
    @test cov(H) ≈ cov(h)

    @test_throws ArgumentError cov( MatrixTDist(1, M, Σ, Ω) )
    @test_throws ArgumentError var( MatrixTDist(1, M, Σ, Ω) )
end

@testset "MatrixTDist logpdf" begin
    #  Check against MvTDist
    @test logpdf(L, X) ≈ logpdf(l, x)
    @test logpdf(H, Y) ≈ logpdf(h, y)

    #  MT(v, M, vΣ, Ω) -> MN(M, Σ, Ω) as v -> ∞
    @test isapprox(logpdf(K, A), logpdf(MN, A), atol = 0.1)
end

@testset "MatrixTDist sample moments" begin
    @test isapprox(mean(rand(D, 100000)), mean(D) , atol = 0.1)
    @test isapprox(cov(hcat(vec.(rand(D, 100000))...)'), cov(D), atol = 0.1)
end

@testset "MatrixTDist conversion" for elty in (Float32, Float64, BigFloat)
    Del1 = convert(MatrixTDist{elty}, D)
    Del2 = convert(MatrixTDist{elty}, v, M, PDΣ, PDΩ, D.logc0)

    @test partype(Del1) == elty
    @test partype(Del2) == elty
end


text:
abstract type SpectralEstimate{D,T} <: AbstractVector{T} end

@doc raw"""
    Periodogram(timeseries, Δ)

Compute the periodogram for the provided timeseries with sampling rate Δ.

# Arguments
- `timeseries`: A `Vector` if univariate and an `n` by `d` `Matrix` if multivariate, where `n` is the number of observations and `d` is the dimension of the timeseries.
- `Δ`: A positive real number.

The periodogram is defined as
```math
\boldsymbol I(ω)&=\boldsymbol J(ω) \boldsymbol J(ω)^H \quad \text{where} \quad \boldsymbol J(ω) = \sqrt{\frac{Δ}{2π n}}\sum_{t=0}^{n-1} \boldsymbol{P}_{tΔ}e^{-itΔ ω}
```

Note the periodogram is in terms of angular frequency here, and uses the normalisation ``Δ/2π``.
The choice of normalisation is essentially arbitrary; however, this matches our definition for the spectral density function.
"""
struct Periodogram{D,T,V} <: SpectralEstimate{D,T}
    Ω::V
    ordinates::Vector{T}
    function Periodogram(timeseries::Matrix{T}, Δ::Real) where {T}
        Δ > 0 || throw(ArgumentError("Δ should be a positive."))
        J = fftshift(fft(timeseries, 1),1)
        n = size(J,1)
        D = size(J,2)
        ordinates = [(J[ii, :]*J[ii, :]'.*(Δ / (2π * n))) for ii ∈ 1:size(J,1)]
        Ω = fftshift(fftfreq(n, 2π/Δ))
        new{D, eltype(ordinates), typeof(Ω)}(Ω, ordinates)
    end
    function Periodogram(timeseries::Vector{T}, Δ::Real) where {T}
        Δ > 0 || throw(ArgumentError("Δ should be a positive."))
        J = fftshift(fft(timeseries, 1),1)
        n = size(J,1)
        ordinates = abs.(J.^2).*(Δ / (2π * n))
        Ω = fftshift(fftfreq(n, 2π/Δ))
        new{1, eltype(ordinates), typeof(Ω)}(Ω, ordinates)
    end
    function Periodogram(Ω, ordinates::Vector{T}) where {T}
        new{size(ordinates[1],1), T, typeof(Ω)}(Ω, ordinates)
    end
end

@doc raw"""
    BartlettPeriodogram(timeseries, Δ, segmentlength)

Compute the Bartlett periodogram for the provided timeseries with sampling rate Δ.

# Arguments
- `timeseries`: A `Vector` if univariate and an `n` by `d` `Matrix` if multivariate, where `n` is the number of observations and `d` is the dimension of the timeseries.
- `Δ`: A positive real number.
- `segmentlength`: the length of series used in each segment.

Computes an estimate of the spectral density function using Bartlett's method. Using the same normalisation as `Periodogram`.

# External links

* [Bartlett's method on Wikipedia](https://en.wikipedia.org/wiki/Bartlett%27s_method)

"""
struct BartlettPeriodogram{D,T,V} <: SpectralEstimate{D,T}
    Ω::V
    ordinates::Vector{T}
    function BartlettPeriodogram(timeseries::Matrix{T}, Δ::Real, segmentlength::Int = Int(100÷Δ)) where {T}
        nsegments = size(timeseries, 1) ÷ segmentlength
        P = Periodogram(timeseries[1:segmentlength, :], Δ)

        for ii = 1:nsegments-1
            P += Periodogram(timeseries[segmentlength*ii.+(1:segmentlength), :], Δ)
        end

        BartlettPeriodogram(P/nsegments)
    end
    function BartlettPeriodogram(timeseries::Vector{T}, Δ::Real, segmentlength::Int = Int(100÷Δ)) where {T}
        nsegments = size(timeseries, 1) ÷ segmentlength
        P = Periodogram(timeseries[1:segmentlength], Δ)

        for ii = 1:nsegments-1
            P += Periodogram(timeseries[segmentlength*ii.+(1:segmentlength)], Δ)
        end

        BartlettPeriodogram(P/nsegments)
    end
    function BartlettPeriodogram(p::Periodogram{D,T}) where {D,T}
        new{D, T, typeof(p.Ω)}(p.Ω, p.ordinates)
    end
    function BartlettPeriodogram(Ω, ordinates::Vector{T}) where {T}
        new{size(ordinates[1],1), T, typeof(Ω)}(Ω, ordinates)
    end
end

struct CoherancyEstimate{D,T,V} <: SpectralEstimate{D,T}
    Ω::V
    ordinates::Vector{T}

    function CoherancyEstimate(timeseries::Matrix{T}, Δ::Real, segmentlength::Int = Int(100÷Δ)) where {T}
        S = BartlettPeriodogram(timeseries, Δ, segmentlength)
        CoherancyEstimate(S)
    end

    function CoherancyEstimate(S::SpectralEstimate{D,T}) where {D,T}
        newordinates = [[s[i,j] / sqrt(s[i,i]*s[j,j]) for i ∈ 1:D, j ∈ 1:D] for s ∈ S.ordinates]
        new{D,eltype(newordinates),typeof(S.Ω)}(S.Ω, newordinates)
    end

    function CoherancyEstimate(Ω, ordinates::Vector{T}) where {T}
        new{size(ordinates[1],1), T, typeof(Ω)}(Ω, ordinates)
    end
end

ndims(::SpectralEstimate{D,T}) where {D,T} = D
size(ŝ::SpectralEstimate) = (length(getfreq(ŝ)),)
getindex(ŝ::T, inds) where {T<:SpectralEstimate} = getconstructor(T)(getfreq(ŝ)[inds], getordinate(ŝ)[inds])
getindex(ŝ::SpectralEstimate, ind::Int) = (getfreq(ŝ)[ind],getordinate(ŝ)[ind])
log10(ŝ::T) where {T<:SpectralEstimate} = getconstructor(T)(getfreq(ŝ), log10.(getordinate(ŝ)))
getfreq(p::Periodogram) = p.Ω
getfreq(b::BartlettPeriodogram) = b.Ω
getfreq(c::CoherancyEstimate) = c.Ω
getordinate(p::Periodogram) = p.ordinates
getordinate(b::BartlettPeriodogram) = b.ordinates
getordinate(c::CoherancyEstimate) = c.ordinates
getconstructor(::Type{<:Periodogram}) = Periodogram
getconstructor(::Type{<:BartlettPeriodogram}) = BartlettPeriodogram
getconstructor(::Type{<:CoherancyEstimate}) = CoherancyEstimate

@recipe function f(ŝ::SpectralEstimate)
    HermitianPlot(getfreq(ŝ), getordinate(ŝ))
end

@recipe function f(ŝ::SpectralEstimate{1,T}) where {T} # univariate plotting
    getfreq(ŝ), getordinate(ŝ)
end

function Base.:+(ŝ₁::T, ŝ₂::T) where {T<:SpectralEstimate}
    getfreq(ŝ₁) == getfreq(ŝ₂) || error("Frequencies must be the same to add spectral estimates.")
    getconstructor(T)(getfreq(ŝ₁), getordinate(ŝ₁) .+ getordinate(ŝ₂))
end

function Base.:/(ŝ::T, a::Real) where {T<:SpectralEstimate}
    getconstructor(T)(getfreq(ŝ), getordinate(ŝ)./a)
end


text:
export GP

###################### Implementation of AbstractGaussianProcess interface #################

"""
    GP{Tμ<:MeanFunction, Tk<:Kernel}

A Gaussian Process (GP) object. Either constructed using an Affine Transformation of
existing GPs or by providing a mean function `μ`, a kernel `k`, and a `GPC` `gpc`.
"""
struct GP{Tμ<:MeanFunction, Tk<:Kernel} <: AbstractGaussianProcess
    args::Any
    μ::Tμ
    k::Tk
    n::Int
    gpc::GPC
    function GP{Tμ, Tk}(args, μ::Tμ, k::Tk, gpc::GPC) where {Tμ, Tk<:Kernel}
        gp = new{Tμ, Tk}(args, μ, k, gpc.n, gpc)
        gpc.n += 1
        return gp
    end
    GP{Tμ, Tk}(μ::Tμ, k::Tk, gpc::GPC) where {Tμ, Tk} = GP{Tμ, Tk}(nothing, μ, k, gpc)
end
GP(μ::Tμ, k::Tk, gpc::GPC) where {Tμ<:MeanFunction, Tk<:Kernel} = GP{Tμ, Tk}(μ, k, gpc)

# GP initialised with a constant mean. Zero is specially handled.
function GP(m::Real, k::Kernel, gpc::GPC)
    if isfinite(length(k))
        if iszero(m)
            return GP(zero(EmpiricalMean(Zeros(length(k)))), k, gpc)  # hack
        else
            return GP(FiniteMean(ConstantMean(m), eachindex(k)), k, gpc)
        end
    else
        if iszero(m)
            return GP(zero(ConstantMean(m)), k, gpc)  # hack
        else
            return GP(ConstantMean(m), k, gpc)
        end
    end
end
GP(m, k::Kernel, gpc::GPC) = GP(CustomMean(m), k, gpc)
GP(k::Kernel, gpc::GPC) = GP(ZeroMean{Float64}(), k, gpc)
function GP(args...)
    μ, k, gpc = μ_p′(args...), k_p′(args...), get_check_gpc(args...)
    return GP{typeof(μ), typeof(k)}(args, μ, k, gpc)
end

function get_check_gpc(args...)
    gpc = args[findfirst(map(arg->arg isa GP, args))].gpc
    @assert all([!(arg isa GP) || arg.gpc == gpc for arg in args])
    return gpc
end

show(io::IO, gp::GP) = print(io, gp)

function Base.print(io::IO, gp::GP)
    println(io, "GP with mean:")
    println(io, Shunted(Shunt(4, ' '), gp.μ))
    println(io, "and kernel:")
    print(io, Shunted(Shunt(4, ' '), gp.k))
end

mean(f::GP) = f.μ

"""
    kernel(f::Union{Real, Function})
    kernel(f::AbstractGP)
    kernel(f::Union{Real, Function}, g::AbstractGP)
    kernel(f::AbstractGP, g::Union{Real, Function})
    kernel(fa::AbstractGP, fb::AbstractGP)

Get the cross-kernel between `GP`s `fa` and `fb`, and . If either argument is deterministic
then the zero-kernel is returned. Also, `kernel(f) === kernel(f, f)`.
"""
kernel(f::GP) = f.k
function kernel(fa::GP, fb::GP)
    @assert fa.gpc === fb.gpc
    if fa === fb
        return kernel(fa)
    elseif fa.args == nothing && fa.n > fb.n || fb.args == nothing && fb.n > fa.n
        if isfinite(length(fa)) && isfinite(length(fb))
            return FiniteZeroCrossKernel(eachindex(fa), eachindex(fb))
        elseif isfinite(length(fa)) && !isfinite(length(fb))
            return LhsFiniteZeroCrossKernel(eachindex(fa))
        elseif !isfinite(length(fa)) && isfinite(length(fb))
            return RhsFiniteZeroCrossKernel(eachindex(fb))
        else # Both processes are infinite dimensional.
            return ZeroKernel{Float64}()
        end
    elseif fa.n > fb.n
        return k_p′p(fa.args..., fb)
    else
        return k_pp′(fa, fb.args...)
    end
end
kernel(::Union{Real, Function}) = ZeroKernel{Float64}()
kernel(::Union{Real, Function}, ::GP) = ZeroKernel{Float64}()
kernel(::GP, ::Union{Real, Function}) = ZeroKernel{Float64}()

##################################### Syntactic Sugar#######################################

import Base: promote, convert

get_zero(p::Int) = get_zero(p, p)
function get_zero(p::Real, q::Real)
    if isfinite(p) && isfinite(q)
        return p == q ? FiniteZeroKernel(1:p) : FiniteZeroCrossKernel(p, q)
    elseif isfinite(p)
        return LhsFiniteZeroCrossKernel(1:p)
    elseif isfinite(q)
        return RhsFiniteZeroCrossKernel(1:q)
    else
        return ZeroKernel{Float64}()
    end
end

promote(f::GP, x::Union{Real, Function, MeanFunction}) = (f, convert(GP, x, f.gpc))
promote(x::Union{Real, Function, MeanFunction}, f::GP) = reverse(promote(f, x))
convert(::Type{<:GP}, x::Real, gpc::GPC) = GP(ConstantMean(x), ZeroKernel{typeof(x)}(), gpc)
convert(::Type{<:GP}, f::Function, gpc::GPC) = GP(CustomMean(f), ZeroKernel{Float64}(), gpc)
convert(::Type{<:GP}, μ::MeanFunction, gpc::GPC) = GP(μ, get_zero(length(μ)), gpc)


text:
using Test
using LinearAlgebra
using QBase
using Base.Iterators
using Suppressor
using CVChannel

@testset "./src/see-saw_optimization.jl" begin

id_kraus_ops = [[1 0;0 1]]
bit_flip_kraus_ops = [[1 0;0 1]*sqrt(0.7), [0 1;1 0]*sqrt(0.3)]

x_states = [[0.5 0.5;0.5 0.5], [0.5 -0.5;-0.5 0.5]]
z_states = [[1 0;0 0],[0 0;0 1]]
trine_states = trine_qubit_states()

@testset "fixedStateCV" begin
    @testset "identity chanel" begin
        z_cv, z_povm = fixedStateCV(z_states, id_kraus_ops)

        @test z_cv ≈ 2 atol=1e-6
        @test is_povm(z_povm, atol=1e-6)
        @test isapprox(z_povm, z_states, atol=1e-6)

        trine_cv, trine_povm = fixedStateCV(trine_states, id_kraus_ops)

        @test trine_cv ≈ 2 atol=1e-6
        @test is_povm(trine_povm, atol=1e-6)
        @test isapprox(trine_povm, 2/3*trine_states, atol=1e-6)
    end

    @testset "bit-flip channel" begin
        x_cv, x_povm = fixedStateCV(x_states, bit_flip_kraus_ops)

        @test x_cv ≈ 2 atol=1e-6
        @test is_povm(x_povm, atol=1e-6)
        @test isapprox(x_povm, x_states, atol=1e-6)

        z_cv, z_povm = fixedStateCV(z_states, bit_flip_kraus_ops)

        @test z_cv ≈ 1.4 atol=1e-6
        @test is_povm(z_povm)
        @test isapprox(z_povm, z_states, atol=1e-6)
    end
end

@testset "fixedMeasurementCV" begin
    @testset "identity channel" begin
        z_cv, opt_z_states = fixedMeasurementCV(z_states, id_kraus_ops)

        @test z_cv ≈ 2 atol=1e-6
        @test all(ρ -> is_density_matrix(ρ, atol=1e-6), opt_z_states)
        @test isapprox(opt_z_states, z_states, atol=1e-6)

        trine_povm = 2/3 * trine_states
        trine_cv, opt_trine_states = fixedMeasurementCV(trine_povm, id_kraus_ops)

        @test trine_cv ≈ 2 atol=1e-6
        @test all(ρ -> is_density_matrix(ρ, atol=1e-6), opt_trine_states)
        @test isapprox(opt_trine_states, trine_states, atol=1e-6)
    end

    @testset "bit-flip channel" begin
        x_cv, opt_x_states = fixedMeasurementCV(x_states, bit_flip_kraus_ops)

        @test x_cv ≈ 2 atol=1e-6
        @test all(ρ -> is_density_matrix(ρ), opt_x_states)
        @test isapprox(opt_x_states, x_states, atol=1e-6)

        z_cv, opt_z_states = fixedMeasurementCV(z_states, bit_flip_kraus_ops)

        @test z_cv ≈ 1.4 atol=1e-6
        @test all(ρ -> is_density_matrix(ρ, atol=1e-6), opt_z_states)
        @test isapprox(opt_z_states, z_states, atol=1e-6)
    end
end

@testset "seesawCV" begin
    d = 3
    anti_sym_choi = wernerState(d, 0) * d
    anti_sym_kraus_ops = map(i -> reshape(anti_sym_choi[:,i], (d,d)), 1:d^2)
    par_anti_sym_kraus_ops = collect(flatten(
        map(k1 -> map(k2 -> kron(k1,k2), anti_sym_kraus_ops), anti_sym_kraus_ops)
    ))

    # verify kraus operators are trace-preserving
    @test sum(k -> k' * k, anti_sym_kraus_ops) ≈ I
    @test sum(k -> k' * k, par_anti_sym_kraus_ops) ≈ I
    @test choi(anti_sym_kraus_ops) ≈ anti_sym_choi

    @testset "singular qutrit anti-symmetric Werner-Holevo channel" begin
        init_states = [
            [1/2 0 0;0 1/3 0;0 0 1/6],[1/6 0 0;0 1/2 0;0 0 1/3],[1/3 0 0;0 1/6 0;0 0 1/2]
        ]

        num_steps = 3
        max_cv_tuple, cvs, opt_ensembles, opt_povms = seesawCV(
            init_states, anti_sym_kraus_ops, num_steps
        )

        @test max_cv_tuple[1] ≈ 1.5 atol=1e-6
        @test all(ρ -> is_density_matrix(ρ), max_cv_tuple[2])
        @test is_povm(max_cv_tuple[3])
        @test length(max_cv_tuple) == 3

        @test length(cvs) == 2 * num_steps
        @test max(cvs...) ≈ 1.5 atol=1e-6
        @test min(cvs...) ≈ 1.25 atol=1e-5

        @test length(opt_ensembles) == num_steps
        @test all(states -> all(ρ -> is_density_matrix(ρ, atol=1e-6), states), opt_ensembles)

        @test length(opt_povms) == num_steps
        @test all(Π -> is_povm(Π, atol=1e-6), opt_povms)
    end

    @testset "parallel qutrit anti-symmetric Werner-Holevo channel" begin
        z_basis_states = computational_basis_states(9)
        max_mix_state = Matrix{Float64}(I, (9,9))/9
        init_states = map(ρ -> ρ/16 + 15/16 * max_mix_state, z_basis_states)

        num_steps = 2
        max_cv_tuple, cvs, opt_ensembles, opt_povms = seesawCV(
            init_states, par_anti_sym_kraus_ops, num_steps
        )

        @test max_cv_tuple[1] ≈ 2.25 atol=1e-5
        @test all(ρ -> is_density_matrix(ρ, atol=1e-5), max_cv_tuple[2])
        @test is_povm(max_cv_tuple[3], atol=1e-5)
        @test length(max_cv_tuple) == 3

        @test length(cvs) == 2 * num_steps
        @test max(cvs...) ≈ 2.25 atol=1e-5
        @test min(cvs...) ≈ 1.078 atol=1e-3

        @test length(opt_ensembles) == num_steps
        @test all(states -> all(ρ -> is_density_matrix(ρ, atol=1e-5), states), opt_ensembles)

        @test length(opt_povms) == num_steps
        @test all(Π -> is_povm(Π, atol=1e-5), opt_povms)
    end

    @testset "verbose printout" begin
        init_states = [
            [1/2 0 0;0 1/3 0;0 0 1/6],[1/6 0 0;0 1/2 0;0 0 1/3],[1/3 0 0;0 1/6 0;0 0 1/2]
        ]

        num_steps = 1
        empty_printout = @capture_out seesawCV(
            init_states, anti_sym_kraus_ops, num_steps
        )
        @test empty_printout == ""

        printout = @capture_out seesawCV(
            init_states, anti_sym_kraus_ops, num_steps, verbose=true
        )

        regex_match = r"^i = 1\nfixed_state_cv = \d\.\d+\nfixed_povm_cv = \d\.\d+\nmax_cv = \d\.\d+\n$"
        @test occursin(regex_match, printout)
    end
end

end

max_forks_repo_name: ChitambarLab/cv-channel


text:
function load_Ham_medium(; inf_medium = false)
    ## ---------------------------------------------------------------
    # Ham's F-12
    # but mediam componets from Human1 Zenodo https://zenodo.org/record/3583004#.Xu8KdmpKiCh
    # task (Growth in Ham's medium)
    # Final concentrations in mM

    # Nutrient Mixture F-12 Ham N6760
    # Source concentration where are g/L, from https://www.sigmaaldrich.com/life-science/cell-culture/learning-center/media-formulations/f-12-ham.html

    Ham_medium = Dict()

    ## ---------------------------------------------------------------
    # Inorganic Salts

    # Calcium Chloride (CaCl2)   0.0333 g/L
    if inf_medium
        Ham_medium["Ca2+[s]"] = INF
        Ham_medium["chloride[s]"] = INF
    else
        Ham_medium["Ca2+[s]"] = 0.0333 * 1e3 / 110.98 # mM
        Ham_medium["chloride[s]"] = 2 * 0.0333 * 1e3 / 110.98 # mM
    end

    # Cupric Sulfate • 5H2O (CuSO4 + 5H2O)   0.0000025 g/L
    if inf_medium
        Ham_medium["Cu2+[s]"] = INF
        Ham_medium["sulfate[s]"] = INF
    else
        Ham_medium["Cu2+[s]"] = 0.0000025 * 1e3 / 249.685 # mM
        Ham_medium["sulfate[s]"] = 0.0000025 * 1e3 / 249.685 # mM
    end

    # Ferrous Sulfate • 7H2O (FeSO4.7H2O)   0.000834 g/L
    if inf_medium
        Ham_medium["Fe2+[s]"] = INF
    else
        Ham_medium["Fe2+[s]"] = 0.000834 * 1e3 / 151.91 # mM
        Ham_medium["sulfate[s]"] += 0.000834 * 1e3 / 151.91 # mM
    end

    # Magnesium Chloride (MgCl2)   0.0576 g/L
    if inf_medium
        Ham_medium["Mg2+[s]"] = INF
    else
        Ham_medium["Mg2+[s]"] = 0.0576 * 1e3 / 95.211 # mM
        Ham_medium["chloride[s]"] += 0.0576 * 1e3 / 95.211 # mM
    end

    # Potassium Chloride (KCl)  0.224 g/L
    if inf_medium
        Ham_medium["K+[s]"] = INF
    else
        Ham_medium["K+[s]"] = 0.224 * 1e3 / 74.5513 # mM
        Ham_medium["chloride[s]"] += 0.224 * 1e3 / 74.5513 # mM
    end

    # Sodium Bicarbonate (NaHCO3)  1.176 g/L
    if inf_medium
        Ham_medium["Na+[s]"] = INF
        Ham_medium["HCO3-[s]"] = INF
    else
        Ham_medium["Na+[s]"] = 1.176 * 1e3 / 84.007 # mM
        Ham_medium["HCO3-[s]"] = 1.176 * 1e3 / 84.007 # mM
    end

    # Sodium Chloride (NaCl)  7.599 g/L
    if !inf_medium
        Ham_medium["Na+[s]"] += 7.599 * 1e3 / 58.44 # mM
        Ham_medium["chloride[s]"] += 7.599 * 1e3 / 58.44 # mM
    end

    # Sodium Phosphate Dibasic (anhydrous) (Na2HPO4)  0.14204 g/L
    if inf_medium
        Ham_medium["Pi[s]"] = INF
    else
        Ham_medium["Pi[s]"] = 0.14204 * 1e3 / 141.96 # mM
        Ham_medium["Na+[s]"] += 2 * 0.14204 * 1e3 / 141.96 # mM
    end

    # Zinc Sulfate • 7H2O (ZnSO4. 7H2O)  0.000863 g/L
    if inf_medium
        Ham_medium["zinc[s]"] = INF
    else
        Ham_medium["zinc[s]"] = 0.000863 * 1e3 / 287.6 # mM
        Ham_medium["sulfate[s]"] += 0.000863 * 1e3 / 287.6 # mM
    end

    # Sodium Bicarbonate  1.176 g/L
    if !inf_medium
        Ham_medium["Na+[s]"] += 1.176 * 1e3 / 58.44 # mM
        Ham_medium["HCO3-[s]"] += 1.176 * 1e3 / 84.007 # mM
    end

    ## ---------------------------------------------------------------
    # Vitamins

    # D-Biotin   0.0000073 g/L
    if inf_medium
        Ham_medium["biotin[s]"] = INF
    else
        Ham_medium["biotin[s]"] = 0.0000073 * 1e3 / 244.31 # mM
    end

    # Choline Chloride (C5H14NO.Cl) 0.01396 g/L
    if inf_medium
        Ham_medium["choline[s]"] = INF
    else
        Ham_medium["choline[s]"] = 0.01396 * 1e3 / 139.62 # mM
        Ham_medium["chloride[s]"] += 0.01396 * 1e3 / 139.62 # mM
    end

    # Folic Acid   0.00132 g/L
    if inf_medium
        Ham_medium["folate[s]"] = INF
    else
        Ham_medium["folate[s]"] = 0.00132 * 1e3 / 441.4 # mM
    end

    # myo-Inositol   0.018 g/L
    if inf_medium
        Ham_medium["inositol[s]"] = INF
    else
        Ham_medium["inositol[s]"] = 0.018 * 1e3 / 180.16 # mM
    end

    # Niacinamide   0.000037 g/L
    if inf_medium
        Ham_medium["nicotinamide[s]"] = INF
    else
        Ham_medium["nicotinamide[s]"] = 0.000037 * 1e3 / 122.12 # mM
    end

    # D-Pantothenic Acid (hemicalcium)   0.00048 g/L
    if inf_medium
        Ham_medium["pantothenate[s]"] = INF
    else
        Ham_medium["pantothenate[s]"] = 0.00048 * 1e3 / 219.23 # mM
    end

    # Pyridoxine • HCl   0.000062 g/L
    if inf_medium
        Ham_medium["pyridoxine[s]"] = INF
    else
        Ham_medium["pyridoxine[s]"] = 0.000062 * 1e3 / 169.18 # mM
        Ham_medium["chloride[s]"] += 0.000062 * 1e3 / 169.18 # mM
    end

    # Riboflavin   0.000038 g/L
    if inf_medium
        Ham_medium["riboflavin[s]"] = INF
    else
        Ham_medium["riboflavin[s]"] = 0.000038 * 1e3 / 376.36 # mM
    end

    # Thiamine • HCl   0.00034 g/L
    if inf_medium
        Ham_medium["thiamin[s]"] = INF
    else
        Ham_medium["thiamin[s]"] = 0.00034 * 1e3 / 265.355 # mM
        Ham_medium["chloride[s]"] += 0.00034 * 1e3 / 265.355 # mM
    end

    # Vitamin B12   0.00136 g/L
    if inf_medium
        Ham_medium["cobamide-coenzyme[s]"] = INF
    else
        Ham_medium["cobamide-coenzyme[s]"] = 0.00136 * 1e3 / 1355.365 # mM
    end

    ## ---------------------------------------------------------------
    # Other
    # D-Glucose   1.802  g/L
    Ham_medium["glucose[s]"] = 1.802 * 1e3 / 180.156 # mM

    # Hypoxanthine   0.00408 g/L
    if inf_medium
        Ham_medium["hypoxanthine[s]"] = INF
    else
        Ham_medium["hypoxanthine[s]"] = 0.00408 * 1e3 / 136.1115 # mM
    end

    # Linoleic Acid   0.000084 g/L
    if inf_medium
        Ham_medium["linolenate[s]"] = INF
    else
        Ham_medium["linolenate[s]"] = 0.000084 * 1e3 / 280.4472 # mM
    end

    # Phenol Red • Na   0.0013 g/L # PH Indicator

    # Putrescine • HCl   0.000161 g/L
    if inf_medium
        Ham_medium["putrescine[s]"] = INF
    else
        Ham_medium["putrescine[s]"] = 0.000161 * 1e3 / 88.15 # mM
        Ham_medium["chloride[s]"] += 0.000161 * 1e3 / 88.15 # mM
    end

    # Pyruvic Acid • Na   0.11 g/L # Carbon source, not included

    # Thioctic Acid   0.00021 g/L
    if inf_medium
        Ham_medium["lipoic acid[s]"] = INF
    else
        Ham_medium["lipoic acid[s]"] = 0.00021 * 1e3 / 206.33 # mM
    end

    # Thymidine   0.00073   g/L
    if inf_medium
        Ham_medium["thymidine[s]"] = INF
    else
        Ham_medium["thymidine[s]"] = 0.00073 * 1e3 / 242.2286 # mM
    end

    ## ---------------------------------------------------------------
    # Add
    # L-Glutamine   0.146 g/L
    Ham_medium["glutamine[s]"] = 0.146 * 1e3 / 146.14 # mM

    ## ---------------------------------------------------------------
    # Amino Acids

    # L-Alanine   0.009 g/L
    Ham_medium["alanine[s]"] = 0.009 * 1e3 / 89.09 # mM

    # L-Arginine • HCl   0.211 g/L
    Ham_medium["arginine[s]"] = 0.211 * 1e3 / 174.2 # mM
    Ham_medium["chloride[s]"] += 0.211 * 1e3 / 174.2 # mM

    # L-Asparagine • H2O   0.01501 g/L
    Ham_medium["asparagine[s]"] = 0.01501 * 1e3 / 132.12 # mM

    # L-Aspartic Acid   0.0133 g/L
    Ham_medium["aspartate[s]"] = 0.0133 * 1e3 / 133.11 # mM

    # L-Cysteine • HCl • H2O   0.035 g/L
    Ham_medium["cysteine[s]"] = 0.035 * 1e3 / 121.16 # mM
    Ham_medium["chloride[s]"] += 0.035 * 1e3 / 121.16  # mM

    # L-Glutamic Acid   0.0147 g/L
    Ham_medium["glutamate[s]"] = 0.0147 * 1e3 / 147.13 # mM

    # L-Glutamine   0.146 g/L
    Ham_medium["glutamine[s]"] += 0.146 * 1e3 / 146.14 # mM

    # Glycine   0.00751 g/L
    Ham_medium["glycine[s]"] = 0.00751 * 1e3 / 75.07 # mM

    # L-Histidine • 3HCl • H2O   0.02096 g/L
    Ham_medium["histidine[s]"] = 0.02096 * 1e3 / 155.1546 # mM
    Ham_medium["chloride[s]"] += 3 * 0.02096 * 1e3 / 155.1546 # mM

    # L-Isoleucine   0.00394 g/L
    Ham_medium["isoleucine[s]"] = 0.00394 * 1e3 / 131.17 # mM

    # L-Leucine   0.0131 g/L
    Ham_medium["leucine[s]"] = 0.0131 * 1e3 / 131.17 # mM

    # L-Lysine • HCl   0.0365 g/L
    Ham_medium["lysine[s]"] = 0.0365 * 1e3 / 146.19 # mM
    Ham_medium["chloride[s]"] += 0.0365 * 1e3 / 146.19 # mM

    # L-Methionine   0.00448 g/L
    Ham_medium["methionine[s]"] = 0.00448 * 1e3 / 149.21 # mM

    # L-Phenylalanine   0.00496 g/L
    Ham_medium["phenylalanine[s]"] = 0.00496 * 1e3 / 165.19 # mM

    # L-Proline   0.0345 g/L
    Ham_medium["proline[s]"] = 0.0345 * 1e3 / 115.13 # mM

    # L-Serine   0.0105 g/L
    Ham_medium["serine[s]"] = 0.0105 * 1e3 / 105.09 # mM

    # L-Threonine   0.0119 g/L
    Ham_medium["threonine[s]"] = 0.0119 * 1e3 / 119.1192 # mM

    # L-Tryptophan   0.00204 g/L
    Ham_medium["tryptophan[s]"] = 0.00204 * 1e3 / 204.23 # mM

    # L-Tyrosine • 2Na • 2H2O   0.00778 g/L
    Ham_medium["tyrosine[s]"] =  0.00778 * 1e3 / 181.19 # mM
    Ham_medium["Na+[s]"] += 2 * 0.00778 * 1e3 / 181.19 # mM

    # L-Valine   0.0117 g/L
    Ham_medium["valine[s]"] = 0.0117 * 1e3 / 117.151 # mM
    # -

    # ## Extras

    Ham_medium["O2[s]"] = INF
    Ham_medium["H2O[s]"] = INF
    # Not found in concentration source
    # TODO search real conc for this mets
    # The concentration of this components was selected to ensure
    # suficient growth rate compared to the experimental values
    Ham_medium["alpha-tocopherol[s]"] = INF
    Ham_medium["gamma-tocopherol[s]"] = INF
    Ham_medium["aquacob(III)alamin[s]"] = INF
    Ham_medium["retinoate[s]"] = INF

    # Testing
    All_components = Any["valine[s]", "putrescine[s]", "H2O[s]", "Fe2+[s]", "arginine[s]", "tryptophan[s]",
                        "asparagine[s]", "glutamine[s]", "serine[s]", "linolenate[s]", "cysteine[s]",
                        "HCO3-[s]", "Ca2+[s]", "gamma-tocopherol[s]", "chloride[s]", "histidine[s]",
                        "hypoxanthine[s]", "pyridoxine[s]", "Na+[s]", "thiamin[s]", "phenylalanine[s]",
                        "riboflavin[s]", "Cu2+[s]", "K+[s]", "lipoic acid[s]", "alanine[s]", "Mg2+[s]",
                        "zinc[s]", "sulfate[s]", "lysine[s]", "glycine[s]", "leucine[s]", "threonine[s]",
                        "tyrosine[s]", "aspartate[s]", "inositol[s]", "thymidine[s]", "cobamide-coenzyme[s]",
                        "choline[s]", "proline[s]", "retinoate[s]", "glutamate[s]", "pantothenate[s]", "folate[s]",
                        "O2[s]", "alpha-tocopherol[s]", "Pi[s]", "isoleucine[s]", "aquacob(III)alamin[s]",
                        "glucose[s]", "nicotinamide[s]", "methionine[s]", "biotin[s]"] |> sort
    @assert all(Ham_medium |> keys |> collect |> sort .== All_components)
    @assert all(c >= 0.0 for (met, c) in Ham_medium)

    return Ham_medium
end


text:
using AstroTime
import Stripeline
const Sl = Stripeline

function test_split_into_n()
    @test split_into_n(20, 3)  == [6, 7, 7]
end

function test_allocate_tod(time_range)
    t0 = minimum(time_range)
    t1 = maximum(time_range)

    # Number of samples acquired during the observation (this is usually huge!)
    nsamples = length(time_range)

    # Three fake polarimeters
    polarimeters = 1:3

    mpi_size = 10  # Number of fake MPI processes

    # Allocate all the TODs in memory, assuming we are running `mpi_size` processes
    tods = [allocate_tod(
        Float32,
        time_range,
        polarimeters,
        mpi_rank = rank,
        mpi_size = mpi_size,
    ) for rank in 0:(mpi_size - 1)]

    # Have the N TODs been really allocated?
    @test length(tods) == 10

    # The first TOD must start when the observation begins…
    @test minimum(tods[1].time_range) == t0
    # …and the last TOD must end when the observation ends
    @test maximum(tods[end].time_range) == t1

    # Check that the total number of samples in each TOD matches what we have asked for
    @test sum([length(cur_tod.time_range) for cur_tod in tods]) == length(time_range)

    # Check that there are no overlaps between consecutive TODs
    for i in 2:length(tods)
        @test maximum(tods[i - 1].time_range) < minimum(tods[i].time_range)
    end

    for i in eachindex(tods)
        cur_shape = size(tods[i].samples)
        @test cur_shape[1] == length(tods[i].time_range)
        @test cur_shape[2] == 8
        @test cur_shape[3] == length(polarimeters)
    end
end

test_split_into_n()

# Run tests using a simple floating-point range
test_allocate_tod(0.0:0.01:1000.0)

# Re-run `test_allocate_tod` using an AstroTime range
begin
    fsamp_hz = 100.0
    τ_s = 1 / fsamp_hz  # Integration time for a sample

    # Fake date when the observation starts
    t0 = TAIEpoch(2022, 1, 1, 12, 0, 0.0)
    # End of the observation
    t1 = t0 + 1hours

    # This is a range that covers *all* the samples during the observation
    time_range = t0 : (τ_s * seconds) : t1

    test_allocate_tod(time_range)
end


text:
# ---
# title: 334. Increasing Triplet Subsequence
# id: problem334
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/increasing-triplet-subsequence/description/>
# hidden: true
# ---
#
# Given an unsorted array return whether an increasing subsequence of length 3
# exists or not in the array.
#
# Formally the function should:
#
# > Return true if there exists _i, j, k_
# >  such that _arr[i]_ < _arr[j]_ < _arr[k]_ given 0 ≤ _i_ < _j_ < _k_ ≤ _n_ -1
# > else return false.
#
# **Note:** Your algorithm should run in O( _n_ ) time complexity and O( _1_ )
# space complexity.
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,4,5]
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: [5,4,3,2,1]
#     Output: false
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export override, OverrideHandler

import Base: show

"""
    OverrideHandler(::Type{T}, value::Union{Nothing, T}) where T
    OverrideHandler(::Type{T}) where T

Handler used in `substitute` operator.

See also: [`substitute`](@ref)
"""
mutable struct OverrideHandler{T}
    value :: Union{Nothing, T}
end

OverrideHandler(::Type{T}, value::Union{Nothing, T}) where T = OverrideHandler{T}(value)
OverrideHandler(::Type{T})                           where T = OverrideHandler{T}(nothing)

setvalue!(handler::OverrideHandler, value) = handler.value = value
getvalue(handler::OverrideHandler)         = handler.value

"""
    override(handler::OverrideHandler)

Creates an override operator that overrides each emission from source observable with value
provided in `handler`. If handler contains `nothing` source observable emits as usual.
For constant override see [`map_to`](@ref). Use `Rocket.setvalue!` to set new value for `handler`.

# Producing

Stream of type `<: Subscribable{Union{L, T}}` where `L` refers to the type of source stream and `T` referes to the type of handler's value

# Examples

```jldoctest
using Rocket

subject = Subject(Int)
handler = OverrideHandler(Int, -1)

source = subject |> override(handler)

subscription = subscribe!(source, logger())

next!(subject, 1)
next!(subject, 2)

Rocket.setvalue!(handler, nothing)

next!(subject, 3)
next!(subject, 4)

Rocket.setvalue!(handler, -2)

next!(subject, 5)
next!(subject, 6)

unsubscribe!(subscription)

# output
[LogActor] Data: -1
[LogActor] Data: -1
[LogActor] Data: 3
[LogActor] Data: 4
[LogActor] Data: -2
[LogActor] Data: -2
```

See also: [`OverrideHandler`](@ref)
"""
override(handler::OverrideHandler{T}) where T = OverrideOperator{T}(handler)

struct OverrideOperator{T} <: InferableOperator
    handler :: OverrideHandler{T}
end

operator_right(::OverrideOperator{T}, ::Type{L}) where { T, L } = Union{T, L}

function on_call!(::Type{L}, ::Type{Union{L, T}}, operator::OverrideOperator{T}, source) where { L, T }
    return proxy(Union{L, T}, source, OverrideProxy{T}(operator.handler))
end

struct OverrideProxy{T} <: ActorProxy
    handler :: OverrideHandler{T}
end

actor_proxy!(::Type{L}, proxy::OverrideProxy{T}, actor::A) where { L, T, A } = OverrideActor{L, T, A}(proxy.handler, actor)

struct OverrideActor{L, T, A} <: Actor{L}
    handler :: OverrideHandler{T}
    actor   :: A
end

function on_next!(actor::OverrideActor{L}, data::L) where L
    override_value = getvalue(actor.handler)
    if override_value !== nothing
        next!(actor.actor, override_value)
    else
        next!(actor.actor, data)
    end
end

on_error!(actor::OverrideActor, err) = error!(actor.actor, err)
on_complete!(actor::OverrideActor)   = complete!(actor.actor)

Base.show(io::IO, ::OverrideOperator)          = print(io, "OverrideOperator()")
Base.show(io::IO, ::OverrideProxy)             = print(io, "OverrideProxy()")
Base.show(io::IO, ::OverrideActor{L})  where L = print(io, "OverrideActor($L)")


text:
# https://adventofcode.com/2020/day/23

using ProgressMeter

example = "389125467"

include("LinkedList.jl")

function play(cups, rounds::Int)
    max = maximum(cups)

    head = tolinkedlist(cups)
    tail = last(head)

    # Without storing the address book with pointers to the
    # nodes, this code is painfully slow, significantly slower
    # than the naive implementation I used to solve the problem.
    # The trick does the job, though it's not my idea, it is
    # borrowed from the code mentioned on Julia's Zulip by
    # Nicolas Viennot (2020-12-23 17:50).

    pointers = Array{LinkedList}(undef, length(cups))
    pos = head
    while !isnothing(pos)
        pointers[pos.value] = pos
        pos = pos.next
    end

    tail.next = head

    p = Progress(rounds)
    for i in 1:rounds
        current = head.value

        picked = peek(head.next, 3)
        head.next = head.next.next.next.next

        destination = current - 1
        destination = destination < 1 ? max : destination
        while destination in picked
            destination -= 1
            destination = destination < 1 ? max : destination
        end

        pos = pointers[destination]
        tail = pos.next
        for x in picked
            pos.next = LinkedList(x, tail)
            pointers[x] = pos.next
            pos = pos.next
        end

        head = head.next
        next!(p)
    end

    n = length(cups)
    return peek(head, n + (n - 1))[n:(n + n - 1)]
end

@assert all(play(Int[3, 8, 9, 1, 2, 5, 4, 6, 7], 10) .== [5, 8, 3, 7, 4, 1, 9, 2, 6])

function part1(input; rounds::Integer=100, verbose=false)
    cups = parse.(Int, split(input, ""))

    cups = play(cups, rounds)

    pos = findfirst(cups .== 1)
    n = length(cups)
    idx = [mod1(pos + i, n) for i in 1:(n - 1)]
    return join(cups[idx])
end

@assert part1(example, rounds=10, verbose=false) == "92658374"
@assert part1(example) == "67384529"

function part2(input; rounds=10_000_000, size=1_000_000, verbose=false)
    cups = parse.(Int, split(input, ""))
    max = maximum(cups)
    for x in (max + 1):size
        push!(cups, x)
    end
    n = length(cups)

    cups = play(cups, rounds)

    pos = findfirst(cups .== 1)
    return (*)([cups[mod1(pos + i, n)] for i in 1:2]...)
end

@assert part2(example) == 149245887792

test = "784235916"
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == "53248976"
@assert result2 == 418819514477

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code