text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.
##################################################################################

"""
    sim_nr!(x, sd, maxiter, tol, verbose [, sparse_solver])

Solve the simulation problem.
  * `x` - the array of data for the simulation. All initial, final and exogenous
    conditions are already in place.
  * `sd::AbstractSolverData` - the solver data constructed for the simulation
    problem.
  * `maxiter` - maximum number of iterations.
  * `tol` - desired accuracy.
  * `verbose` - whether or not to print progress information.
  * `sparse_solver` (optional) - a function called to solve the linear system A
    x = b for x. Defaults to A\\b

"""
function sim_nr!(x::AbstractArray{Float64}, sd::StackedTimeSolverData,
    maxiter::Int64, tol::Float64, verbose::Bool,
    sparse_solver::Function = (A, b) -> A \ b)
    for it = 1:maxiter
        @timer Fx, Jx = global_RJ(x, x, sd)
        @timer nFx = norm(Fx, Inf)
        if nFx < tol
            if verbose
                @info "$it, || Fx || = $(nFx)"
            end
            break
        end
        @timer Δx = sparse_solver(Jx, Fx)
        @timer nΔx = norm(vec(Δx), Inf)
        @timer assign_update_step!(x, -1.0, Δx, sd)
        if verbose
            @info "$it, || Fx || = $(nFx), || Δx || = $(nΔx)"
        end
        if nΔx < tol
            break
        end
    end
    return nothing
end

export simulate
"""
    simulate(model, plan, data; <options>)

Run a simulation for the given model, simulation plan and exogenous data.

### Arguments
  * `model` - the [`Model`](@ref) instance to simulate.
  * `plan` - the [`Plan`](@ref) for the simulation.
  * `data` - a 2D `Array` containing the exogenous data. This includes the
    initial and final conditions.

### Options as keyword arguments
  * `fctype::`[`FCType`](@ref) - set the desired final condition type for the
    simulation. The default value is `fcgiven`. Other possible values
    include `fclevel` and `fcslope`.
  * `initial_guess::AbstractMatrix{Float64}` - a 2D `Array` containing the
    initial guess for the solution. This is used to start the Newton-Raphson
    algorithm. The default value is an empty array (`zeros(0,0)`), in which case
    we use the exogenous data for the initial condition. You can use the steady
    state solution using [`steadystatearray`](@ref).
  * `deviation::Bool` - set to `true` if the `data` is given in deviations from
    the steady state. In this case the simulation result is also returned as a
    deviation from the steady state. Default value is `false`.
  * `anticipate::Bool` - set to `false` to instruct the solver that all shocks
    are unanticilated by the agents. Default value is `true`.
  * `verbose::Bool` - control whether or not to print progress information.
    Default value is taken from `model.options`.
  * `tol::Float64` - set the desired accuracy. Default value is taken from
    `model.options`.
  * `maxiter::Int` - algorithm fails if the desired accuracy is not reached
    within this maximum number of iterations. Default value is taken from
    `model.options`.

### See also:

### Examples

"""
function simulate end

function simulate(m::Model, p::Plan, exog_data::AbstractArray{Float64,2};
    initial_guess::AbstractArray{Float64,2} = zeros(0, 0),
    deviation::Bool = false,
    anticipate::Bool = true,
    verbose::Bool = m.options.verbose,
    tol::Float64 = m.options.tol,
    maxiter::Int64 = m.options.maxiter,
    fctype = getoption(m, :fctype, fcgiven),
    expectation_horizon::Union{Nothing,Int64} = nothing,
    sparse_solver::Function = (A, b) -> A \ b
)
    NT = length(p.range)
    nauxs = length(m.auxvars)
    if size(exog_data) != (NT, length(m.varshks))
        error("Incorrect dimensions of exog_data. Expected $((NT, length(m.varshks))), got $(size(exog_data)).")
    end
    if !isempty(initial_guess) && size(initial_guess) != (NT, length(m.varshks))
        error("Incorrect dimensions of initial_guess. Expected $((NT, length(m.varshks))), got $(size(exog_data)).")
    end

    if deviation
        exog_data = copy(exog_data)
        local logvars = islog.(m.varshks) .| isneglog.(m.varshks)
        local ss_data = steadystatearray(m, p)
        exog_data[:, logvars] .*= ss_data[:, logvars]
        exog_data[:, .!logvars] .+= ss_data[:, .!logvars]
    end

    exog_data = @timer ModelBaseEcon.update_auxvars(transform(exog_data, m), m)

    if !isempty(initial_guess)
        x = @timer ModelBaseEcon.update_auxvars(transform(initial_guess, m), m)
    else
        x = copy(exog_data)
    end

    if anticipate
        @timer gdata = StackedTimeSolverData(m, p, fctype)
        @timer assign_exog_data!(x, exog_data, gdata)
        if verbose
            @info "Simulating $(p.range[1 + m.maxlag:NT - m.maxlead])" # anticipate gdata.FC
        end
        @timer sim_nr!(x, gdata, maxiter, tol, verbose, sparse_solver)
    else # unanticipated shocks
        init = 1:m.maxlag
        term = NT .+ (1-m.maxlead:0)
        sim = 1+m.maxlag:NT-m.maxlead
        nvars = length(m.variables)
        nshks = length(m.shocks)
        nauxs = length(m.auxvars)
        shkinds = nvars .+ (1:nshks)
        auxinds = nvars .+ nshks .+ (1:nauxs)
        varshkinds = 1:(nvars+nshks)
        allvarinds = 1:(nvars+nshks+nauxs)
        x[init, allvarinds] .= exog_data[init, allvarinds]
        # x[term, allvarinds] .= exog_data[term, allvarinds]
        x[sim, shkinds] .= 0.0
        t0 = first(sim)
        T = last(sim)
        if expectation_horizon === nothing
            # the original code, where we always simulate until the end with the true final conditions
            for t in sim
                exog_inds = p[t, Val(:inds)]
                psim = Plan(m, t:T)
                if (t != t0) && (psim[t0, Val(:inds)] == exog_inds) && (maximum(abs, exog_data[t, shkinds]) < tol)
                    continue
                end
                setexog!(psim, t0, exog_inds)
                @timer gdata = StackedTimeSolverData(m, psim, fctype)
                x[t, exog_inds] = exog_data[t, exog_inds]
                # @timer assign_exog_data!(x[psim.range,:], exog_data[psim.range,:], gdata)
                sim_range = UnitRange{Int}(psim.range)
                xx = view(x, sim_range, :)
                assign_final_condition!(xx, exog_data[sim_range, :], gdata)
                if verbose
                    @info "Simulating $(p.range[t:T])" # anticipate expectation_horizon gdata.FC
                end
                @timer sim_nr!(xx, gdata, maxiter, tol, verbose, sparse_solver)
                # if verbose
                #     @info "Result at $t: " x[[t], :]
                # end
            end
        else
            # the new code, where the first and last simulations use the true
            # simulation range and final condition, while the intermediate
            # simulations use expectation_horizon steps with fcnatural
            if expectation_horizon == 0
                expectation_horizon = length(sim)
            elseif expectation_horizon < 10 * m.maxlead
                @warn "Expectation horizon may be too short for this model. Consider setting it to at least $(10 * m.maxlead)."
            end
            x = [x; zeros(expectation_horizon, size(x, 2))]
            ninit = length(init)
            nterm = length(term)
            # first simulation
            let t = t0
                # first run is with the full range, the true fctype,
                # and only the first period is imposed
                exog_inds = p[t, Val(:inds)]
                psim = Plan(m, t:T)
                setexog!(psim, t0, exog_inds)
                @timer sdata = StackedTimeSolverData(m, psim, fctype)
                x[t, exog_inds] = exog_data[t, exog_inds]
                sim_range = UnitRange{Int}(psim.range)
                xx = view(x, sim_range, :)
                assign_final_condition!(xx, exog_data[sim_range, :], sdata)
                if verbose
                    @info "Simulating $(p.range[t:T])" # anticipate expectation_horizon sdata.FC
                end
                sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver)
                # if verbose
                #     @info "Result at $t: " x[[t], :]
                # end
            end
            # intermediate simulations
            last_t::Int64 = t0
            psim = Plan(m, 0:expectation_horizon-1)
            sdata = StackedTimeSolverData(m, psim, fcnatural)
            for t in sim[2:end]
                exog_inds = p[t, Val(:inds)]
                # we need to run a simulation if a variable is exogenous, or if a shock value is not zero
                # these intermediate simulations are always with fcnatural,
                #       have length equal to expectation_horizon and
                #       only the first period is imposed
                if (exog_inds == shkinds) && (maximum(abs, exog_data[t, shkinds]) <= tol)
                    continue
                end
                setexog!(psim, t0, exog_inds)
                update_plan!(sdata, m, psim)
                # note that the range always goes from 0 to expectation_horizon-1,
                # so we need to add t in order to get the correct set of rows of x
                sim_range = t .+ UnitRange{Int}(psim.range)
                xx = view(x, sim_range, :)
                # The initial conditions are already set
                # The exogenous values are already set as well, except for the first period
                # In other words, we only need to impose the first period here
                xx[t0, exog_inds] = exog_data[t, exog_inds]
                # Update the final conditions (the second argument is not used with fcnatural)
                @timer assign_final_condition!(xx, xx, sdata)
                if verbose
                    @info "Simulating $(p.range[t] .+ (0:expectation_horizon - 1))" # anticipate expectation_horizon sdata.FC
                end
                @timer sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver)
                # if verbose
                #     @info "Result at $t: " x[[t], :]
                # end
                last_t = t  # keep track of last simulation time
            end
            # last simulation
            if last_t > t0
                # do we need to re-run the last simulation?
                # if it didn't reach T, then yes
                # if the final condition is not fcnatural, then yes
                if (last_t + expectation_horizon < T) || (fctype != fcnatural)
                    psim = Plan(m, min(last_t + 1, T):T)
                    # there are no unanticipated shocks in this simulation
                    sdata = StackedTimeSolverData(m, psim, fctype)
                    # the initial conditions and the exogenous data are already in x
                    # we only need the final conditions
                    sim_range = UnitRange{Int}(psim.range)
                    xx = view(x, sim_range, :)
                    assign_final_condition!(xx, exog_data[sim_range, :], sdata)
                    if verbose
                        @info "Simulating $(p.range[last_t + 1:T])" # anticipate expectation_horizon sdata.FC
                    end
                    @timer sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver)
                    # if verbose
                    #     @info "Result at $(last_t + 1): " x[[last_t + 1], :]
                    # end
                end
            end
            # x = x[begin:end-expectation_horizon, :]
        end
    end

    x = x[axes(exog_data)...]
    x .= inverse_transform(x, m)
    if deviation
        x[:, logvars] ./= ss_data[:, logvars]
        x[:, .!logvars] .-= ss_data[:, .!logvars]
    end

    return x
end

# The versions of simulate with Dict/Workspace/SimData

@inline simulate(m::Model, p::Plan, data::Dict; kwargs...) = simulate(m, p, dict2data(data, m, p; copy = true); kwargs...)
@inline simulate(m::Model, p::Plan, data::Workspace; kwargs...) = simulate(m, p, workspace2data(data, m, p; copy = true); kwargs...)

function simulate(m::Model, p::Plan, data::SimData; kwargs...)
    exog = data2array(data, m, p)
    initial_guess = get(kwargs, :initial_guess, nothing)
    if initial_guess isa SimData
        kw = (; initial_guess = data2array(initial_guess, m, p))
    elseif initial_guess isa Workspace
        kw = (; initial_guess = workspace2data(initial_guess, m, p))
    elseif initial_guess isa AbstractDict
        kw = (; initial_guess = dict2array(initial_guess, m, p))
    else
        kw = (;)
    end
    result = copy(data)
    result[p.range, m.varshks] .= simulate(m, p, exog; kwargs..., kw...)
    return result
end


text:
## Qt has a QTimer class for calling values periodically.
## Here we illustrate its basic use.

## a timer instance
timer = QtCore[:QTimer]()

## add a callback
i = 1
qconnect(timer, :timeout) do
  global i
  i += 1
  i > 100 ? qinvoke(timer, :stop) : println(i)
end

qinvoke(timer, :start, 50)

ext: jl


text:
using SpecialFunctions: erfcx
import RationalApproximations: locate_extrema_idx, RationalFunc, aaa_v2, minimax
using PyPlot

G(t) = (1-t)/(1+t)
α = 0.5
β = 1.0
fα(t) = erfcx(G(t))

t = range(-1, 1, length=800)
max_m = 5
tol = 1e-10
F = fα.(t)
r, err = aaa_v2(F, t, max_m, tol)
m = length(r)
resid = F - r[m].(t)
idx = locate_extrema_idx(resid, m)

figure(1)
plot(t, resid, t[idx], resid[idx], "o")
grid(true)
xlabel(L"$t$", fontsize=12)
xylims=axis()
savefig("figure13L.pdf")

supp_pt = r[m].supp_pt
max_iterations=3
tx = t[idx]
opt_r, zmin, zmax = minimax(fα, tx, supp_pt, (-1.0, 1.0), max_iterations)
resid = F - opt_r[end].(t)
figure(2)
plot(t, resid, tx, fα.(tx) - opt_r[end].(tx), "o")
grid(true)
xlabel(L"$t$", fontsize=12)
axis(xylims)
savefig("figure13R.pdf")


text:
# Create a "generated" function to add together two ModelStates. This results in
# a function for each type that is as fast as a custom + function.
@generated function Base.:+(x::T, y::T) where {T <: ModelState}
    expr = ( :(x.$f + y.$f) for f in fieldnames(T) ) # Make a tuple of expressions to add each field
    :( T($(expr...)) ) # Splat them into an expression for a constructor for a T type
end

# Create a generated function to multiple a ModelState by a scalar.
@generated function Base.:*(x::S, y::T) where {S <: Real, T <: ModelState}
    expr = ( :(x * y.$f) for f in fieldnames(T) ) # Make a tuple of expressions to add each field
    :( T($(expr...)) ) # Splat them into an expression for a constructor for a T type
end

# We use this to ignore nothing states when adding them together.
Base.:+(a::Nothing, b)          = b
Base.:+(a, b::Nothing)          = a
Base.:+(a::Nothing, b::Nothing) = a
Base.:*(a, b::Nothing)          = b

# There are built-in types for UT date-times, but they don't work the way that I
# want. This isn't mutable. Just make a new one when something needs to be
# different.
struct UtcDateTime
    Y::Int
    M::Int
    D::Int
    h::Int
    m::Int
    s::Float64
    # UtcDateTime(v::Vector{Number}) UtcDateTime(v[1], v[2], v[3], v[4], v[5], v[6]);
end

# Create a function to convert a vector to a UtcDateType.
Base.convert(::Type{UtcDateTime}, d::Vector{T}) where {T <: Real} = UtcDateTime(d...);


text:
"""
ArithmeticBlocks
"""

export AddBlock

mutable struct AddBlock <: AbstractArithmeticBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    #op::Vector{Symbol}

    function AddBlock() #op::Vector{Symbol}
        #@createblock new(Vector{InPort}(), Vector{OutPort}(), op) length(op) 1
        @createblock new(Vector{InPort}(), Vector{OutPort}()) 2 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::AddBlock)
    #println(io, "AddBlock($(blk.op))")
    println(io, "AddBlock()")
end

"""
to expr
"""

function _toexpr(blk::AddBlock)
    #args = [Expr(:call, x[1], _toexpr(x[2])) for x = zip(blk.op, blk.inport)]
    #Expr(:call, :+, args...)
    Expr(:call, :+, _toexpr(blk.inport[1]), _toexpr(blk.inport[2]))
end


text:
using ForecastPlots
using Test

@testset "ForecastPlots.jl" begin

    x = rand(100)

    #splot
    px = splot(x;title="SPlot", plot=false)
    @test size(px) == (9,12)

    px = splot(x,"day";title="SPlot", plot=false)
    @test size(px) == (15,7)

    px = splot(x,"quarter",title="SPlot", plot=false)
    @test size(px) == (25,4)

    px = splot(x,"month",title="SPlot", plot=false)
    @test size(px) == (9,12)

    px = splot(x,10,title="SPlot", plot=false)
    @test size(px) == (10,10)

    # acf
    px = acf(x; type="cor", plot=false)
    @test size(px) == (20,)

    px = acf(x; type="cov", plot=false)
    @test size(px) == (20,)

    @test_throws AssertionError acf(x; type="", plot=false)

    # ccf
    x1 = rand(100);
    x2 = circshift(x1,6);
    px = ccf(x1, x2; type="cor", plot = false)

    @test size(px) == (41,)
    @test findmax(px)[2] == 27

    # pacf
    px = pacf(x; plot=false)
    @test size(px) == (20,2)

    px = pacf(x; type = "stepwise", plot=false)
    @test size(px) == (20,)

    px = pacf(x; type = "real", plot=false)
    @test size(px) == (20,)

    px = pacf(x; type = "stepwise-real", plot=false)
    @test size(px) == (20,2)

    @test_throws AssertionError pacf(x; type = "", plot=false)

end


text:
module MzXML

using Base64
using LightXML, Unitful, ProgressMeter
using MzCore
using MzCore.IntervalSets, MzCore.AxisArrays
using CodecZlib

"""
    MSscan(polarity::Char, msLevel::Int, retentionTime::typeof(1.0u"s"),
           lowMz::Float64, highMz::Float64, basePeakMz::Float64,
           totIonCurrent::Float64, mz::Vector, I::Vector, children::Vector{MSscan})

Represent mass spectrometry data at a single time point, scanning over all `m/z`.
`polarity` is `+` or `-` for positive-mode or negative-mode, respectively.
`msLevel` is 1 for MS, 2 for MS², and so on.
Retention time is represented using the Unitful package.

`mz` and `I` should be vectors of the same length, with an ion count of `I[i]` at `mz[i]`.
`children` refers to MS² or higher for particular ions in this scan.
"""
struct MSscan{T<:AbstractFloat,TI<:Real}
    polarity::Char
    msLevel::Int
    retentionTime::typeof(1.0u"s")
    lowMz::Float64
    highMz::Float64
    basePeakMz::Float64
    totIonCurrent::Float64
    mz::Vector{T}
    I::Vector{TI}
    children::Vector{MSscan{T,TI}}
end

MzCore.intensitytype(::Type{MSscan{T,TI}}) where {T<:AbstractFloat,TI<:Real} = TI
MzCore.mztype(::Type{MSscan{T,TI}}) where {T<:AbstractFloat,TI<:Real} = T

const empty32 = MSscan{Float32,Float32}[]
const empty64 = MSscan{Float64,Float64}[]

function Base.show(io::IO, scan::MSscan)
    if scan.msLevel > 1
        print(io, "  "^(scan.msLevel-2), "└─basePeak ", scan.basePeakMz, ": ")
    end
    print(io, scan.retentionTime, ": ")
    mzI = [mz=>I for (mz, I) in zip(scan.mz, scan.I)]
    println(io, mzI)
end

"""
    scanpos = MzXML.index(filename)

Return the vector of file positions corresponding to each scan.
"""
function index(filename)
    scanpositions = open(filename) do file
        # Find the indexOffset element
        seekend(file)
        p = position(file)
        l = 200
        if p < l
            error("Too short to be mzXML")
        end
        skip(file, -l)
        str = loadrest(file)
        m = match(r"<indexOffset>([0-9].*)</indexOffset>", str)
        if m == nothing
            error("Cannot find indexOffset element")
        end
        length(m.captures) == 1 || error("indexOffset should contain a single number")
        indexpos = parse(Int, m.captures[1])
        length_tail = l - m.offset
        # Read the index
        seek(file, indexpos)
        str = loadrest(file)
        xindex = parse_string(str[1:end-length_tail-1])
        xroot = root(xindex)
        scanpositions = Int[]
        for c in child_elements(xroot)
            name(c) == "offset" || error("index could not be parsed")
            push!(scanpositions, parse(Int, content(c)))
        end
        scanpositions
    end
    scanpositions
end

"""
    scans, info = MzXML.load(filename; productlevels=0, timeinterval=0.0u"s" .. Inf*u"s", timeshift=0.0u"s")

Load a `*.mzXML` file. `scans` is a vector of scan structures, each at a particular scan time.

`productlevels` controls the number of product-ion levels that will be loaded in combined MS/MS² or similar runs.
The "base" MS level (the first encountered in the file) will always be loaded--`productlevels` controls the
number of additional levels.
By default no product-ion data will be loaded (`productlevels=0`), but you can request it by increasing `productlevels` to 1 or higher.

`scans` will include only those scans occuring within `timeinterval`, with the default to include all times.
The `..` is from the `IntervalSets` package.
"""
function load(filename; productlevels=0, kwargs...)
    xdoc = parse_file(filename)
    xroot = root(xdoc)
    if name(xroot) != "mzXML"
        error("Not an mzXML file")
    end

    # Find the msRun node
    msRun = find_element(xroot, "msRun")

    props = Dict{Symbol,Any}()
    props[:startTime] = parse_time(attribute(msRun, "startTime"))
    props[:endTime] = parse_time(attribute(msRun, "endTime"))
    el = find_element(msRun, "msInstrument")
    el = find_element(el, "msModel")
    props[:msModel] = attribute(el, "value")
    props[:productlevels] = productlevels

    scans = load_scans(msRun, productlevels; kwargs...)
    LightXML.free(xdoc)
    return scans, props
end

function load_scans(elm, productlevels=0; kwargs...)
    # First we discover the data type
    local T
    local TI
    for c in child_elements(elm)
        n = name(c)
        if n != "scan"
            continue
        end
        peak = find_element(c, "peaks")
        TI, T, nochildren = precisiondict[attribute(peak, "precision")]
        break
    end
    # Now load the data
    scans = MSscan{T,TI}[]
    load_scans!(scans, elm, productlevels; kwargs...)
end

@noinline function load_scans!(scans, elm, productlevels; kwargs...)
    prog = ProgressUnknown("Number of scans read:")
    for c in child_elements(elm)
        n = name(c)
        if n != "scan"
            continue
        end
        scan = load_scan(c, productlevels; kwargs...)::Union{Nothing,eltype(scans)}
        if scan !== nothing
            push!(scans, scan)
            next!(prog)
        end
    end
    finish!(prog)
    scans
end

polaritydict = Dict("+" => '+', "-" => '-')
precisiondict = Dict("32" => (Float32, Float32, empty32), "64" => (Float64, Float64, empty64))

function load_scan(elm, productlevels; timeinterval=0.0u"s" .. Inf*u"s", timeshift=0.0u"s")
    polarity = polaritydict[attribute(elm, "polarity")]
    retentionTime = parse_time(attribute(elm, "retentionTime")) + timeshift
    retentionTime ∈ timeinterval || return nothing
    lMza, hMza = attribute(elm, "lowMz"), attribute(elm, "highMz")
    if lMza !== nothing && hMza !== nothing
        lowMz, highMz = parse(Float64, lMza), parse(Float64, hMza)
    else
        lowMz = highMz = NaN
    end
    msLevela = attribute(elm, "msLevel")
    msLevel = msLevela === nothing ? 1 : parse(Int, msLevela)
    basePeakMz = parse(Float64, attribute(elm, "basePeakMz"))
    totIonCurrent = parse(Float64, attribute(elm, "totIonCurrent"))
    npeaks = parse(Int, attribute(elm, "peaksCount"))
    peak = find_element(elm, "peaks")
    data = base64decode(content(peak))
    compression_type = attribute(peak, "compressionType")
    if compression_type == "zlib"
        data = transcode(ZlibDecompressor, data)
    elseif compression_type != "none"
        error("compressionType $compression_type not supported")
    end
    TI, T, nochildren = precisiondict[attribute(peak, "precision")]
    A = reinterpret(TI, data)
    bo = attribute(peak, "byteOrder")
    if bo == "network"
        ntoh!(A)
    else
        error("Don't know what to do with byteOrder $bo")
    end
    po = attribute(peak, "pairOrder")
    if po === nothing
        po = attribute(peak, "contentType")
    end
    po == "m/z-int" || error("Don't know what to do with pairOrder/contentType $po")
    I = A[2:2:end]
    mz = reinterpret(T, A[1:2:end])
    if !issorted(mz)   # it's unclear whether the standard requires these to be sorted, so better to check
        perm = sortperm(mz)
        I, mz = I[perm], mz[perm]
    end
    children = productlevels > 0 ? load_scans!(MSscan{T,TI}[], elm, productlevels-1) : nochildren
    MSscan{T,TI}(polarity, msLevel, retentionTime, lowMz, highMz, basePeakMz, totIonCurrent, mz, I, children)
end

function loadrest(file::IOStream)
    nb = filesize(file) - position(file)
    b = read(file, nb)
    return String(b)
end

function parse_time(tstr)
    if startswith(tstr, "PT") && endswith(tstr, "S")
        return parse(Float64, tstr[3:end-1])*u"s"
    end
    error("Time string $tstr not recognized")
end

function ntoh!(A)
    for i = 1:length(A)
        A[i] = ntoh(A[i])
    end
end

Base.pairs(scan::MSscan) = zip(scan.mz, scan.I)

function val2index(ax::Axis, val)
    f, s = first(ax.val), step(ax.val)
    r = (val - f) / s
    idx = floor(Int, r)
    return idx + 1, r - idx
end

"""
    copyto!(dest::AxisArray, src::Vector{<:MSscan})

Convert the "sparse-matrix" representation `src` of mass spectrometry
data into a dense 2d matrix `dest`. `dest` must have one axis named
`:mz` and the other `:time`; the `:time` axis must have physical units of
time.

Any previous contents of `dest` are erased. The values of `src` are
accumulated into `dest` using bilinear interpolation.
"""
function Base.copyto!(dest::AxisArray, src::Vector{<:MSscan})
    axmz, axt = AxisArrays.axes(dest, Axis{:mz}), AxisArrays.axes(dest, Axis{:time})
    fill!(dest, 0)
    for scan in src
        j, ft = val2index(axt, scan.retentionTime)
        firstindex(axt.val) - 1 <= j <= lastindex(axt.val) || continue
        for (mz, I) in pairs(scan)
            i, fmz = val2index(axmz, mz)
            # Accumulate into `dest`, using bilinear interpolation to
            # avoid aliasing
            if firstindex(axt.val) <= j <= lastindex(axt.val)
                if firstindex(axmz.val) <= i <= lastindex(axmz.val)
                    dest[axmz(i),axt(j)] += (1-fmz)*(1-ft)*I
                end
                if firstindex(axmz.val) <= i+1 <= lastindex(axmz.val)
                    dest[axmz(i+1),axt(j)] += fmz*(1-ft)*I
                end
            end
            if firstindex(axt.val) <= j+1 <= lastindex(axt.val)
                if firstindex(axmz.val) <= i <= lastindex(axmz.val)
                    dest[axmz(i),axt(j+1)] += (1-fmz)*ft*I
                end
                if firstindex(axmz.val) <= i+1 <= lastindex(axmz.val)
                    dest[axmz(i+1),axt(j+1)] += fmz*ft*I
                end
            end
        end
    end
    return dest
end

# An internal type for nested `extrema` calls (yielding the min of all mins, max of all maxs)
struct ExPair{T}
    min::T
    max::T
end
Base.min(a::ExPair, b::ExPair) = a.min <= b.min ? a : b
Base.max(a::ExPair, b::ExPair) = a.max >= b.max ? a : b

function MzCore.limits(scans::Vector{<:MSscan})
    function extrema_mz(scan)
        isfinite(scan.lowMz) && isfinite(scan.highMz) && return ExPair(scan.lowMz, scan.highMz)
        return ExPair(first(scan.mz), last(scan.mz))  # safe because guaranteed to be sorted
    end
    tmin, tmax   = extrema(scan->scan.retentionTime, scans)
    mzmin, mzmax = extrema(extrema_mz, scans)
    return Axis{:mz}(mzmin.min .. mzmax.max), Axis{:time}(tmin..tmax)
end

end


text:
include("fibonacci_app.jl")
app = fibonacci_app(false)
run_server(app)


text:
module Counter

struct Counter
    items::Dict{Any, Integer}
end

"""
    Bag() -> Bag
"""
function Counter()::Bag
    Bag(Dict())
end

"""
    length(bag) -> Integer

Return the number of items in the given bag.
"""
function length(bag::Bag)::Integer
    size = 0
    for items in bag.values
        size += item[2]
    end
    size
end

"""
    isempty(bag) -> Bool

Return if the given is empty.
"""
function isempty(bag::Bag)::Integer
    length(bag) == 0
end

"""
    add!(bag, item)

Add the given item to the given bag.
"""
function add!(bag::Bag, item)
    !haskey(bag.items, item) ? bag.items[item] = 1 : bag.items[item] += 1
    return
end

"""
    remove!(bag, item) -> item

Remove and return an instance of the given item from the given bag
"""
function remove!(bag::Bag, item)
    if haskey(bag.items, item)
        bag.items[item] -= 1
        if bag.items[item] == 0
            delete!(bag.items, item)
        end
    end
end

function iterate(bag::Bag)
    (iterate(keys(bag.items)))
end

"""
    ∈(bag, item) -> Bool

Return if the given bag contains the given item.
"""
function ∈(bag::Bag, item)::Bool
    item in bag.items
end

"""
    countitem(bag, item) -> Integers

Return the count of the given item in the given bag.
"""
function countitem(bag::Bag, item)::Integer
    count = 0
    for bagitem in bag.items
        item == bagitem[1] ? count = bagitem[2] : nothing
    end
    count
end


text:
using Leiden
using Test
using Random: Random
using SparseArrays: spzeros
using StatsBase: summarystats

include("reference.jl")

function generate_two_community()
    weights = spzeros(6, 6)
    edges = [
        (1, 2),
        (2, 3),
        (3, 1),

        (4, 5),
        (5, 6),
        (6, 4),

        (1, 4),
    ]
    for (u, v) in edges
        weights[u,v] = weights[v,u] = 1
    end
    return weights
end

function generate_karate_club()
    # derived from http://konect.cc/networks/ucidata-zachary/
    weights = spzeros(34, 34)
    edges = [
        ( 1,  2), ( 1,  3), ( 2,  3), ( 1,  4), ( 2,  4),
        ( 3,  4), ( 1,  5), ( 1,  6), ( 1,  7), ( 5,  7),
        ( 6,  7), ( 1,  8), ( 2,  8), ( 3,  8), ( 4,  8),
        ( 1,  9), ( 3,  9), ( 3, 10), ( 1, 11), ( 5, 11),
        ( 6, 11), ( 1, 12), ( 1, 13), ( 4, 13), ( 1, 14),
        ( 2, 14), ( 3, 14), ( 4, 14), ( 6, 17), ( 7, 17),
        ( 1, 18), ( 2, 18), ( 1, 20), ( 2, 20), ( 1, 22),
        ( 2, 22), (24, 26), (25, 26), ( 3, 28), (24, 28),
        (25, 28), ( 3, 29), (24, 30), (27, 30), ( 2, 31),
        ( 9, 31), ( 1, 32), (25, 32), (26, 32), (29, 32),
        ( 3, 33), ( 9, 33), (15, 33), (16, 33), (19, 33),
        (21, 33), (23, 33), (24, 33), (30, 33), (31, 33),
        (32, 33), ( 9, 34), (10, 34), (14, 34), (15, 34),
        (16, 34), (19, 34), (20, 34), (21, 34), (23, 34),
        (24, 34), (27, 34), (28, 34), (29, 34), (30, 34),
        (31, 34), (32, 34), (33, 34),
    ]
    for (u, v) in edges
        weights[u,v] = weights[v,u] = 1
    end
    return weights
end

# Relative error.
relerror(x, y) = abs(x - y) / max(abs(x), abs(y))

@testset "Leiden.jl" begin
    Random.seed!(1234)
    γ = 0.25
    result = Leiden.leiden(generate_two_community(), resolution = γ)
    @test result.quality == 6 * (1 - γ)
    @test result.partition == [[1,2,3], [4,5,6]]

    Random.seed!(1234)
    γ = 0.25
    result = Leiden.louvain(generate_two_community(), resolution = γ)
    @test result.quality == 6 * (1 - γ)
    @test result.partition == [[1,2,3], [4,5,6]]

    Random.seed!(1234)
    γ = 0.05
    karate = generate_karate_club()
    stats = summarystats([Leiden.leiden(karate, resolution = γ).quality for _ in 1:1000])
    stats_ref = summarystats([Reference.leiden(karate, resolution = γ).quality for _ in 1:1000])
    #@show relerror(stats.mean, stats_ref.mean) relerror(stats.max, stats_ref.max)
    #@show relerror(stats.min, stats_ref.min) relerror(stats.median, stats_ref.median)
    @test relerror(stats.mean,   stats_ref.mean)   < 1e-3
    @test relerror(stats.max,    stats_ref.max)    < 1e-3
    @test relerror(stats.min,    stats_ref.min)    < 1e-3
    @test relerror(stats.median, stats_ref.median) < 1e-3

    Random.seed!(1234)
    γ = 0.05
    karate = generate_karate_club()
    stats = summarystats([Leiden.louvain(karate, resolution = γ).quality for _ in 1:1000])
    stats_ref = summarystats([Reference.louvain(karate, resolution = γ).quality for _ in 1:1000])
    #@show relerror(stats.mean, stats_ref.mean) relerror(stats.max, stats_ref.max)
    #@show relerror(stats.min, stats_ref.min) relerror(stats.median, stats_ref.median)
    @test relerror(stats.mean,   stats_ref.mean)   < 1e-3
    @test relerror(stats.max,    stats_ref.max)    < 1e-3
    #@test relerror(stats.min,    stats_ref.min)    < 1e-3
    @test relerror(stats.median, stats_ref.median) < 1e-3
end


text:
@testset "$(rpad("Gauss Tableaus",80))" begin

    function _getTableauGLRK1(T=Float64)
        a = reshape([1//2], (1,1))
        b = [1//1]
        c = [1//2]
        o = 2

        Tableau{T}(:Gauss1, o, a, b, c)
    end

    function _getTableauGLRK2(T=Float64)
        a = [
                [1/4       1/4-√3/6]
                [1/4+√3/6  1/4     ]
            ]
        b = [1/2,      1/2     ]
        c = [1/2-√3/6, 1/2+√3/6]
        o = 4

        Tableau{T}(:Gauss2, o, a, b, c)
    end

    function _getTableauGLRK3(T=Float64)
        a = [
                [5/36         2/9-√15/15  5/36-√15/30]
                [5/36+√15/24  2/9         5/36-√15/24]
                [5/36+√15/30  2/9+√15/15  5/36       ]
            ]
        b = [5/18,        4/9,        5/18       ]
        c = [1/2-√15/10,  1/2,        1/2+√15/10 ]
        o = 6

        Tableau{T}(:Gauss3, o, a, b, c)
    end

    function _getTableauGLRK4(T=Float64)
        a = reshape([
                [(-√30/144 + 1/8)
                 (-13*√42/336 - √10/16 - √3/12 + √(-60*√30 + 450)/144 + √35/28 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
                 (-√35/28 - √10/16 - √3/12 + √(-60*√30 + 450)/144 + 13*√42/336 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
                 (-√35/14 - √(60*√30 + 450)/144 - √42/168 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
                [(-√(60*√30 + 450)/144 - √3/12 + √10/16 + √35/28 + 13*√42/336 + √(2*√30 + 15)/8)/√(2*√30 + 15)
                 (√30/144 + 1/8)
                 (-√35/14 + √42/168 + √(-60*√30 + 450)/144 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
                 (-13*√42/336 - √35/28 - √(60*√30 + 450)/144 - √3/12 + √10/16 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
                [(-√10/16 - √(60*√30 + 450)/144 + √3/12 + √35/28 + 13*√42/336 + √(2*√30 + 15)/8)/√(2*√30 + 15)
                 (-√42/168 + √(-60*√30 + 450)/144 + √(-2*√30 + 15)/8 + √35/14)/√(-2*√30 + 15)
                 (√30/144 + 1/8)
                 (-13*√42/336 - √35/28 - √10/16 - √(60*√30 + 450)/144 + √3/12 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
                [(-√(60*√30 + 450)/144 + √42/168 + √35/14 + √(2*√30 + 15)/8)/√(2*√30 + 15)
                 (-13*√42/336 + √(-60*√30 + 450)/144 + √3/12 + √10/16 + √35/28 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
                 (-√35/28 + √(-60*√30 + 450)/144 + √3/12 + √10/16 + 13*√42/336 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
                 (-√30/144 + 1/8)]
                ], (4,4))'

        b = [
                49/(216+12*√30),
                49/(216-12*√30),
                49/(216-12*√30),
                49/(216+12*√30)
            ]

        c = [
                1/2-√(3/7+2*√30/35)/2,
                1/2-√(3/7-2*√30/35)/2,
                1/2+√(3/7-2*√30/35)/2,
                1/2+√(3/7+2*√30/35)/2
            ]

        o = 8

        Tableau{T}(:Gauss4, o, a, b, c)
    end

    function _getTableauGLRK5(T=Float64)
        a = reshape([
                [-13*√70/3600 + 161/1800
                 -59*√(20*√70 + 350)/19440 - 113*√(14*√70 + 245)/34020 - 11*√(-20*√70 + 350)/6480 + 4*√(-14*√70 + 245)/2835 + 13*√70/3600 + 161/1800
                 -92*√(14*√70 + 245)/8505 + 4*√(20*√70 + 350)/1215 + 32/225
                 -59*√(20*√70 + 350)/19440 - 113*√(14*√70 + 245)/34020 - 4*√(-14*√70 + 245)/2835 + 11*√(-20*√70 + 350)/6480 + 13*√70/3600 + 161/1800
                 -2*√(14*√70 + 245)/315 - 13*√70/3600 + √(20*√70 + 350)/360 + 161/1800];
                [-113*√(-14*√70 + 245)/34020 - 13*√70/3600 + 4*√(14*√70 + 245)/2835 + 11*√(20*√70 + 350)/6480 + 59*√(-20*√70 + 350)/19440 + 161/1800
                 13*√70/3600 + 161/1800
                 -92*√(-14*√70 + 245)/8505 - 4*√(-20*√70 + 350)/1215 + 32/225
                 -2*√(-14*√70 + 245)/315 - √(-20*√70 + 350)/360 + 13*√70/3600 + 161/1800
                 -11*√(20*√70 + 350)/6480 - 113*√(-14*√70 + 245)/34020 - 13*√70/3600 - 4*√(14*√70 + 245)/2835 + 59*√(-20*√70 + 350)/19440 + 161/1800];
                [-23*√(20*√70 + 350)/5760 - 13*√70/3600 + 161/1800 + 11*√(14*√70 + 245)/1440
                 13*√70/3600 + 23*√(-20*√70 + 350)/5760 + 11*√(-14*√70 + 245)/1440 + 161/1800
                 32/225
                 -11*√(-14*√70 + 245)/1440 - 23*√(-20*√70 + 350)/5760 + 13*√70/3600 + 161/1800
                 -11*√(14*√70 + 245)/1440 - 13*√70/3600 + 161/1800 + 23*√(20*√70 + 350)/5760];
                [-59*√(-20*√70 + 350)/19440 - 13*√70/3600 + 4*√(14*√70 + 245)/2835 + 113*√(-14*√70 + 245)/34020 + 11*√(20*√70 + 350)/6480 + 161/1800
                 13*√70/3600 + √(-20*√70 + 350)/360 + 2*√(-14*√70 + 245)/315 + 161/1800
                 4*√(-20*√70 + 350)/1215 + 92*√(-14*√70 + 245)/8505 + 32/225
                 13*√70/3600 + 161/1800
                 -59*√(-20*√70 + 350)/19440 - 11*√(20*√70 + 350)/6480 - 13*√70/3600 - 4*√(14*√70 + 245)/2835 + 113*√(-14*√70 + 245)/34020 + 161/1800];
                [-√(20*√70 + 350)/360 - 13*√70/3600 + 161/1800 + 2*√(14*√70 + 245)/315
                 -11*√(-20*√70 + 350)/6480 + 4*√(-14*√70 + 245)/2835 + 13*√70/3600 + 113*√(14*√70 + 245)/34020 + 59*√(20*√70 + 350)/19440 + 161/1800
                 -4*√(20*√70 + 350)/1215 + 32/225 + 92*√(14*√70 + 245)/8505
                 -4*√(-14*√70 + 245)/2835 + 11*√(-20*√70 + 350)/6480 + 13*√70/3600 + 113*√(14*√70 + 245)/34020 + 59*√(20*√70 + 350)/19440 + 161/1800
                 -13*√70/3600 + 161/1800]
            ], (5,5))'

        b = [
                1/(4600/729+1300*√70/5103),
                1/(4600/729-1300*√70/5103),
                64/225,
                1/(4600/729-1300*√70/5103),
                1/(4600/729+1300*√70/5103)
            ]

        c = [
                1/2-√(5/9+2*√70/63)/2,
                1/2-√(5/9-2*√70/63)/2,
                1/2,
                1/2+√(5/9-2*√70/63)/2,
                1/2+√(5/9+2*√70/63)/2
            ]

        o = 10

        Tableau{T}(:Gauss5, o, a, b, c)
    end

    function _getTableauGLRK6(T=Float64)
        a = reshape([
                [+0.0428311230947925862600740355431832233817056253710109955996588599497364402,
                 -0.0147637259971974124753725910605206514420001878003408840631515394063317448,
                 +0.00932505070647775119143888450800314858828816322281640330116548464376180630,
                 -0.00566885804948351190092125641621650656214351346713436514621339322133564885,
                 +0.00285443331509933513092928583011602153367099528247654439080898011903375405,
                 -0.000812780171264762112299135651562540066903951464973606518543200348171357290];
                [+0.0926734914303788631865122917633203161433485534502333679517368248345204102,
                 +0.0901903932620346518924583784594290279153804731866863705724348100592850095,
                 -0.0203001022932395859524940805242724601067259796418104502230488160090068057,
                 +0.0103631562402464237307199458065599778725040265335003155034035217850387283,
                 -0.00488719292803767146341420376578964407137592740187386736001668777722125914,
                 +0.00135556105548506177551787075080010874364457167567922814822699418129216854];
                [+0.0822479226128438738077716511411289215554396236655129640997987178532804788,
                 +0.196032162333245006055759781563801382788760315019545751706035725357898777,
                 +0.116978483643172761847467585997387748702913901442302633827906329990978550,
                 -0.0204825277456560976298590118654006438219872235078237396303198827234514190,
                 +0.00798999189966233579720442148030827079362758194981571971201211474213850758,
                 -0.00207562578486633419359528915758164772805951363936400480613058750798236207];
                [+0.0877378719744515067137433602439480944914707643813859960054483074074552416,
                 +0.172390794624406967987712335438549785037133364423557021432857505376431512,
                 +0.254439495032001621324794183860176141227815026392429007286132542705408521,
                 +0.116978483643172761847467585997387748702913901442302633827906329990978551,
                 -0.0156513758091757022708430246449433269579993686461730105611661052393287579,
                 +0.00341432357674129871237641994523752520797162707650902709951900204619240137];
                [+0.0843066851341001107446302003355663380197666790663427630510907257181807085,
                 +0.185267979452106975248330960684647699902136873775246608504886307895791278,
                 +0.223593811046099099964215226188215519533323776351104952152409138196918370,
                 +0.254257069579585109647429252519047957512553782526415717878861475990963909,
                 +0.0901903932620346518924583784594290279153804731866863705724348100592850092,
                 -0.00701124524079369066636422067695386937993730270821137675241910493504753033];
                [+0.0864750263608499346324472067379289868303152022069955977178609202476442335,
                 +0.177526353208969968653987471088742034297089951090896196754060639999536273,
                 +0.239625825335829035595856428410992003967971316351739632802026053203292742,
                 +0.224631916579867772503496287486772348817539639661788864354647175338195299,
                 +0.195144512521266716260289347979378707272761134173713625208021159524901768,
                 +0.0428311230947925862600740355431832233817056253710109955996588599497364397]
            ], (6,6))'

        b = [
                0.0856622461895851725201480710863664467634112507420219911993177198994728803,
                0.180380786524069303784916756918858055830760946373372741144869620118570019,
                0.233956967286345523694935171994775497405827802884605267655812659981957101,
                0.233956967286345523694935171994775497405827802884605267655812659981957101,
                0.180380786524069303784916756918858055830760946373372741144869620118570019,
                0.0856622461895851725201480710863664467634112507420219911993177198994728804
            ]

        c = [
                0.0337652428984239860938492227530026954326171311438550875637251917366932495891,
                0.169395306766867743169300202490047326496775717802414964592736647073908251704,
                0.380690406958401545684749139159644032290694684929989324909302417712862532843,
                0.619309593041598454315250860840355967709305315070010675090697582287137467157,
                0.830604693233132256830699797509952673503224282197585035407263352926091748261,
                0.966234757101576013906150777246997304567382868856144912436274808263306750376
            ]

        o = 12

        Tableau{T}(:Gauss6, o, a, b, c)
    end

    @test TableauGauss(1) ≈ _getTableauGLRK1()
    @test TableauGauss(2) ≈ _getTableauGLRK2()
    @test TableauGauss(3) ≈ _getTableauGLRK3()
    @test TableauGauss(4) ≈ _getTableauGLRK4()
    @test TableauGauss(5) ≈ _getTableauGLRK5()
    @test TableauGauss(6) ≈ _getTableauGLRK6()

end


text:
function _transform(x::Array{nf_elem}, y::fmpz_mat)
  n = length(x)
  @assert n == nrows(y)
  m = ncols(y)
  z = Array{nf_elem}(m)
  for i in 1:m
    z[i] = x[1]^y[1, i]
    for j in 2:n
      z[i] = z[i]*x[j]^y[j, i]
    end
  end
  return z
end

function _make_row_primitive(x::fmpz_mat, j::Int)
  y = x[j, 1]
  for i in 1:ncols(x)
    y = gcd(y, x[j, i])
  end
  if y > 1
    for i in 1:ncols(x)
      x[j, i] = div(x[j, i], y)
    end
  end
end

function _make_row_primitive!(x::Array{fmpz, 1})
  y = x[1]
  for i in 2:length(x)
    y = gcd(y, x[i])
    if y == 1
      return x
    end
  end
  if y > 1
    for i in 1:ncols(x)
      x[i] = div(x[i], y)
    end
  end
end


text:
function compute_kinetic_limit(parameter_dictionary::Dict{String,Any})::Float64

    try

        # grab the parameter values from the dictionary -
        E1 = parameter_dictionary["PFK_concentration"]
        ATP = parameter_dictionary["ATP_concentration"]
        F6P = parameter_dictionary["F6P_concentration"]
        kcat = parameter_dictionary["kcat"]
        K_ATP = parameter_dictionary["K_ATP"]
        K_F6P = parameter_dictionary["K_F6P"]

        # compute the kinetic limit -
        kinetic_limit = kcat * E1 * (ATP / (K_ATP + ATP)) * (F6P / (K_F6P + F6P))

        # return -
        return kinetic_limit
    catch error
        throw(error)
    end
end

function compute_activity_variable(effector::Float64, parameter_dictionary::Dict{String,Any})::Float64

    try

        # grab the parameter values from the dictionary -
        W1 = parameter_dictionary["W1"]
        W2 = parameter_dictionary["W2"]
        KD = parameter_dictionary["KD"]
        n = parameter_dictionary["n"]

        # compute the f -
        f = ((effector / KD)^n) / (1 + (effector / KD)^n)

        # compute the v-variable -
        v_variable = (W1 + W2 * f) / (1 + W1 + W2 * f)

        # return -
        return v_variable
    catch error
        throw(error)
    end
end

function compute_activity_variable(effector_array::Array{Float64,1}, parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    v_variable_array = Array{Float64,1}()

    try

        # process each effector -
        for effector_value in effector_array

            # compute -
            value = compute_activity_variable(effector_value, parameter_dictionary)

            # cache -
            push!(v_variable_array, value)
        end

        # return -
        return v_variable_array
    catch error
        throw(error)
    end
end


text:
# SPDX-License-Identifier: X11
# 2020-08-30
# Crane and Turtle (200pt)

((X, Y) -> println((Y & 1 == 0 && Y / 4 <= X <= Y / 2) ? "Yes" : "No "))(
  parse.(Int, split(readline(), " "))...)

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# Code written by Christopher Coey

# Reads a SDPA file to create a MathProgBase conic model
# Detects an LP block but does not currently detect SOC or SOCRotated cones
# SDPA sparse format modified for integer constraints described at
#   http://www.opt.tu-darmstadt.de/scipsdp/downloads/data_format.txt

using GZip

function MPB.loadproblem!(m::MPB.AbstractConicModel, filename::AbstractString; all_bin=false)
    if endswith(filename,".dat-s") || endswith(filename,".dat-s.gz")
       loadsdpa!(m, filename)
    else
       error("unrecognized input format extension in $filename")
    end
end

function loadsdpa!(m::MPB.AbstractConicModel, filename::String; all_bin=false)
    if endswith(filename, ".gz")
        fd = gzopen(filename, "r")
    else
        fd = open(filename, "r")
    end
    (c, A, b, concones, varcones, vartypes) = sdpatoconicdata(fd, all_bin)
    close(fd)

    MathProgBase.loadproblem!(m, c, A, b, concones, varcones)
    MathProgBase.setvartype!(m, vartypes)
end

function sdpatoconicdata(io::IO, all_bin::Bool)
    numvars = 0
    numblocks = 0
    sizeblocks = Vector{Int}()
    c = Vector{Float64}()
    indline = 0

    while !eof(io)
        line = strip(readline(io))

        if !startswith(line, '*') && !startswith(line, '#') && !startswith(line, '\"') && length(line) > 0
            indline += 1
            spline = split(line)
            if indline == 1
                numvars = parse(Int, spline[1])
            elseif indline == 2
                numblocks = parse(Int, spline[1])
            elseif indline == 3
                sizeblocks = [parse(Int, spline[block]) for block in 1:numblocks]
            elseif indline == 4
                c = [parse(Float64, spline[col]) for col in 1:numvars]
                break
            end
        end
    end

    LPblocks = findall(size -> (sign(size) == -1), sizeblocks)
    LPblock = isempty(LPblocks) ? 0 : first(LPblocks)
    PSDblocks = findall(size -> (sign(size) == 1), sizeblocks)

    if LPblock != 0
        sizeblocks[LPblock] *= -1
    end

    numcons = sum(binomial((sizeblocks[block] + 1), 2) for block in PSDblocks)
    if LPblock != 0
        numcons += sizeblocks[LPblock]
    end
    b = zeros(numcons)
    A = spzeros(numcons, numvars)

    bPSDs = SparseMatrixCSC{Float64,Int64}[spzeros(size, size) for size in sizeblocks]
    APSDs = SparseMatrixCSC{Float64,Int64}[spzeros(size, size) for size in sizeblocks, col in 1:numvars]

    while !eof(io)
        line = strip(readline(io))

        if !startswith(line, '*') && !startswith(line, '#') && !startswith(line, '"') && length(line) > 0
            linevec = split(line)
            col = parse(Int, linevec[1])
            block = parse(Int, linevec[2])
            i = parse(Int, linevec[3])
            j = parse(Int, linevec[4])
            v = parse(Float64, linevec[5])

            if block == LPblock
                @assert i == j
                if col == 0
                    b[i] = -v
                else
                    A[i,col] = -v
                end
            else
                if col == 0
                    bPSDs[block][i,j] = -v
                else
                    APSDs[block,col][i,j] = -v
                end
            end
        elseif startswith(line, "*INTEGER")
            break
        end
    end

    dropzeros!(A)

    function ind_svec(n::Int, i::Int, j::Int)::Int
        indmin = min(i,j)
        indmax = max(i,j)

        accum = indmax + 1 - indmin
        for k in 1:(indmin - 1)
            accum += n + 1 - k
        end
        return accum
    end

    varcones = Tuple{Symbol,Vector{Int}}[(:Free, collect(1:numvars))]

    concones = Tuple{Symbol,Vector{Int}}[]
    if LPblock != 0
        rowend = sizeblocks[LPblock]
        push!(concones, (:NonNeg, collect(1:rowend)))
    else
        rowend = 0
    end

    sqrt2 = sqrt(2)

    for block in 1:numblocks
        if block == LPblock
            continue
        end

        rowendprev = rowend

        dropzeros!(bPSDs[block])
        (I,J,V) = findnz(bPSDs[block])
        for nzind in 1:length(V)
            row = rowendprev + ind_svec(sizeblocks[block], I[nzind], J[nzind])
            b[row] = (I[nzind] == J[nzind]) ? V[nzind] : V[nzind]*sqrt2
        end

        for col in 1:numvars
            dropzeros!(APSDs[block,col])
            (I,J,V) = findnz(APSDs[block,col])
            for nzind in 1:length(V)
                row = rowendprev + ind_svec(sizeblocks[block], I[nzind], J[nzind])
                A[row,col] = (I[nzind] == J[nzind]) ? V[nzind] : V[nzind]*sqrt2
            end
        end

        rowend = rowendprev + binomial((sizeblocks[block] + 1), 2)
        push!(concones, (:SDP, collect(rowendprev+1:rowend)))
    end

    vartypes = fill(:Cont, numvars)

    while !eof(io)
        line = strip(readline(io))

        if startswith(line, '*')
            col = parse(Int, split(strip(line[2:end]))[1])
            vartypes[col] = all_bin ? :Bin : :Int
        end
    end

    return (c, A, b, concones, varcones, vartypes)
end


text:
"""
    ColorChannels

Defines methods for conveniently swapping between arrays with color elements and channels.
"""
module ColorChannels

using AxisIndices
using AxisIndices.Interface: NamedDimsArray
using Reexport
using FFTW

@reexport using FixedPointNumbers
@reexport using Colors

if isdefined(ColorTypes, :XRGB) && isdefined(ColorTypes, :RGB1)
    Base.@deprecate_binding RGB1 XRGB
    Base.@deprecate_binding RGB4 RGBX
end
# backward compatibility for ColorTypes < v0.9
if !isdefined(ColorTypes, :XRGB)
    const XRGB = RGB1
    const RGBX = RGB4
end

using MappedArrays
using MappedArrays: AbstractMultiMappedArray, AbstractMappedArray

using AxisIndices.PaddedViews
using .ColorTypes: colorant_string
using Colors: Fractional

using Base: tail, @pure, Indices
import Base: float

# TODO: just use .+
# See https://github.com/JuliaLang/julia/pull/22932#issuecomment-330711997
plus(r::AbstractUnitRange, i::Integer) = broadcast(+, r, i)
plus(a::AbstractArray, i::Integer) = a .+ i

using .ColorTypes: AbstractGray, TransparentGray, Color3, Transparent3
Color1{T} = Colorant{T,1}
Color2{T} = Colorant{T,2}
Color4{T} = Colorant{T,4}
AColor{N,C,T} = AlphaColor{C,T,N}
ColorA{N,C,T} = ColorAlpha{C,T,N}
const NonparametricColors = Union{RGB24,ARGB32,Gray24,AGray32}
Color1Array{C<:Color1,N} = AbstractArray{C,N}
# Type that arises from reshape(reinterpret(To, A), sz):
const RRArray{To,From,N,M,P} = Base.ReshapedArray{To,N,Base.ReinterpretArray{To,M,From,P}}
const RGArray = Union{Base.ReinterpretArray{<:AbstractGray,M,<:Number,P}, Base.ReinterpretArray{<:Number,M,<:AbstractGray,P}} where {M,P}

# delibrately not export these constants to enable extensibility for downstream packages
const NumberLike = Union{Number,AbstractGray}
const Pixel = Union{Number,Colorant}
const GenericGrayImage{T<:NumberLike,N} = AbstractArray{T,N}
const GenericImage{T<:Pixel,N} = AbstractArray{T,N}

export
    StackedView,
    channel_view,
    color_view,
    raw_view,
    normed_view,
    reinterpretc,
    # conversions
#    float16,
    float32,
    float64,
    n0f8,
    n6f10,
    n4f12,
    n2f14,
    n0f16,
    # mapping values
    clamp01,
    clamp01!,
    clamp01nan,
    clamp01nan!,
    color_signed,
    scale_minmax,
    scale_signed,
    take_map,
    # spatial
    # channels
    channeldim,
    channel_axis,
    channel_axis_type,
    channel_keys,
    channel_indices,
    nchannel,
    select_channel,
    zeroarray

include("channels.jl")
include("stackedviews.jl")
include("convert_reinterpret.jl")
include("map.jl")
include("functions.jl")

is_channel(x::Symbol) = x === :channel || x === :Channel || x === :Color || x === :color
AxisIndices.@defdim channel is_channel

"""
    raw_view(img::AbstractArray{FixedPoint})

returns a "view" of `img` where the values are interpreted in terms of
their raw underlying storage. For example, if `img` is an `Array{N0f8}`,
the view will act like an `Array{UInt8}`.

See also: [`normed_view`](@ref)
"""
raw_view(a::AbstractArray{T}) where {T<:FixedPoint} = mappedarray(reinterpret, y->T(y,0), a)
raw_view(a::Array{T}) where {T<:FixedPoint} = reinterpret(FixedPointNumbers.rawtype(T), a)
raw_view(a::AbstractArray{T}) where {T<:Real} = a

"""
    normed_view([T], img::AbstractArray{Unsigned})

returns a "view" of `img` where the values are interpreted in terms of
`Normed` number types. For example, if `img` is an `Array{UInt8}`, the
view will act like an `Array{N0f8}`.  Supply `T` if the element
type of `img` is `UInt16`, to specify whether you want a `N6f10`,
`N4f12`, `N2f14`, or `N0f16` result.

See also: [`raw_view`](@ref)
"""
normed_view(::Type{T}, a::AbstractArray{S}) where {T<:FixedPoint,S<:Unsigned} = mappedarray(y->T(y,0),reinterpret, a)
normed_view(::Type{T}, a::Array{S}) where {T<:FixedPoint,S<:Unsigned} = reinterpret(T, a)
normed_view(::Type{T}, a::AbstractArray{T}) where {T<:Normed} = a
normed_view(a::AbstractArray{UInt8}) = normed_view(N0f8, a)
normed_view(a::AbstractArray{T}) where {T<:Normed} = a

# PaddedViews support
# This make sure Colorants as `fillvalue` are correctly filled, for example, let
# `PaddedView(ARGB(0, 0, 0, 0), img)` correctly filled with transparent color even when
# `img` is of eltype `RGB`
function PaddedViews.filltype(::Type{FC}, ::Type{C}) where {FC<:Colorant, C<:Colorant}
    # rand(RGB, 4, 4) has eltype RGB{Any} but it isn't a concrete type
    # although the consensus[1] is to not make a concrete eltype, this op is needed to make a
    # type-stable colorant construction in _filltype without error; there's no RGB{Any} thing
    # [1]: https://github.com/JuliaLang/julia/pull/34948
    T = eltype(C) === Any ? eltype(FC) : eltype(C)
    return _filltype(FC, base_colorant_type(C){T})
end
_filltype(::Type{<:Colorant}, ::Type{C}) where {C<:Colorant} = C
_filltype(::Type{FC}, ::Type{C}) where {FC<:Color3, C<:AbstractGray} =
    base_colorant_type(FC){promote_type(eltype(FC), eltype(C))}
_filltype(::Type{FC}, ::Type{C}) where {FC<:TransparentColor, C<:AbstractGray} =
    alphacolor(FC){promote_type(eltype(FC), eltype(C))}
_filltype(::Type{FC}, ::Type{C}) where {FC<:TransparentColor, C<:Color3} =
    alphacolor(C){promote_type(eltype(FC), eltype(C))}

# Support transpose
Base.transpose(a::AbstractMatrix{C}) where {C<:Colorant} = permutedims(a, (2,1))
function Base.transpose(a::AbstractVector{C}) where C<:Colorant
    ind = axes(a, 1)
    out = similar(Array{C}, (oftype(ind, Base.OneTo(1)), ind))
    outr = reshape(out, ind)
    copy!(outr, a)
    out
end

end


text:
function restrict(x̅::AbstractVector, n̅, n)
    s = @. (n̅ - n) ÷ 2
    cart = CartesianIndices(@. UnitRange(1+s, n+s))
    lin = LinearIndices(n̅)
    ind = reshape(lin[cart], prod(n))
    view(x̅, ind)
end

function restrict(A̅::AbstractMatrix, n̅, n)
    s = @. (n̅ - n) ÷ 2
    cart = CartesianIndices(@. UnitRange(1+s, n+s))
    lin = LinearIndices(n̅)
    ind = reshape(lin[cart], prod(n))
    view(A̅, ind, ind)
end


text:
module ChargedParticle3d

    include("charged_particle_3d/singular_field.jl")
    include("charged_particle_3d/symmetric_field.jl")
    include("charged_particle_3d/solovev_iter.jl")
    include("charged_particle_3d/solovev_iter_xpoint.jl")
    include("charged_particle_3d/theta_pinch_canonical.jl")
    include("charged_particle_3d/theta_pinch_noncanonical.jl")
    include("charged_particle_3d/tokamak_iter_cylindrical.jl")
    include("charged_particle_3d/tokamak_small_noncanonical.jl")
    include("charged_particle_3d/tokamak_small_cartesian.jl")
    include("charged_particle_3d/tokamak_small_cylindrical.jl")
    include("charged_particle_3d/tokamak_small_toroidal.jl")

end


text:
using Test

Test.@testset "Package Bayesian_EMR_MSM" begin
    include("estimations/test_estimations.jl")
end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
cli.jl defines the `jlcoverage.cli` module containing functions for the `jlcoverage` CLI.

Notes
-----
* CLI functions are defined in a .jl file so that testing and code quality tools can by
  applied to the CLI source code.
"""
module cli

# --- Exports

export parse_args, run

# --- Imports

# Standard library
using Logging

# External packages
using ArgParse
using Coverage
using CoverageTools

# Local modules
using ..jlcoverage

# --- Functions/Methods

"""
    parse_args(; raw_args::Vector{<:AbstractString}=ARGS)::Dict

Parse and return CLI arguments contained in `raw_args`. By default, `raw_args` is set to
`ARGS`, the command-line arguments provided to the executable that called `parse_args()`.
"""
function parse_args(; raw_args::Vector{<:AbstractString}=ARGS)::Dict

    # Define command-line arguments
    description = "Generate coverage analysis report."
    arg_table = ArgParse.ArgParseSettings(; prog="jlcoverage", description=description)
    ArgParse.@add_arg_table! arg_table begin
        "--verbose", "-v"
        help = "enable verbose mode"
        action = :store_true

        "--version", "-V"
        help = "show version and exit"
        action = :store_true

        "paths"
        help =
            "files and directories to include coverage analysis. If `paths` is empty, " *
            "a coverage report is generated for all Julia source files contained in " *
            "(1) the `src` directory if the current directory is a Julia package " *
            "(i.e., the current directory contains `Project.toml` and `src`) or " *
            "(2) the current directory if it is not a Julia package."
        nargs = '*'
    end

    # Parse command-line arguments
    args::Dict = ArgParse.parse_args(raw_args, arg_table)
    args["paths"] = convert(Vector{String}, args["paths"])

    return args
end

"""
    run(paths::Vector; <keyword arguments>)

Run code coverage analysis for files and directories in `paths`.

# Keyword Arguments

* `verbose::Bool`: print more output to the console. Default: `false`
"""
function run(paths::Vector; verbose::Bool=false)
    # --- Check arguments

    # Ensure that `paths` contains strings
    paths = convert(Vector{String}, paths)

    # --- Preparations

    # Handle edge case
    if isempty(paths)
        if isfile("Project.toml") && isdir("src")
            @info "Detected Julia package. Generating report for files in `src` directory."
            paths = ["src"]
        else
            message =
                "Julia package not detected. Generating report for files in current " *
                "directory."
            @info message
            paths = ["."]
        end
    end

    # Set log level
    if !verbose
        disable_logging(Logging.Info)
    end

    # --- Process coverage data and display results

    # Process `*.cov` files
    coverage = Vector{CoverageTools.FileCoverage}()
    for path in paths
        # Ensure paths are absolute paths
        if !isabspath(path)
            path = abspath(path)
        end

        # Process coverage files
        if isdir(path)
            coverage = CoverageTools.merge_coverage_counts(
                coverage, Coverage.process_folder(path::String)
            )
        elseif isfile(path)
            coverage = CoverageTools.merge_coverage_counts(
                coverage, [Coverage.process_file(path)]
            )
        else
            @warn "$path not found. Skipping..."
            continue
        end
    end

    # Display coverage results
    display_coverage(coverage)

    return nothing
end

end  # End of jlcoverage.cli module


text:
export UpdateColley

"""
    UpdateColley

 Colley update rule, e.g., see "Whos's #1", Langville and Meyer, p.21

## Parameters (none)
```
"""
struct UpdateColley <: UpdateRule
# this is a batch calculation that takes no account of past ratings, and
# has no parameters
end

function update_info( rule::UpdateColley )
    info = Dict(
                :name => "Colley",
                :reference => "\"Whos's #1\", Langville and Meyer, p.21",
                :computation => "simultaneous",
                :state_model => "none",
                :input => "outcome",
                :output => "deterministic",
                :model => "single",
                :ties => false,
                :factors => false,
                :parameters => [],
                :record => false
                )
    return info
end
update_info( ::Type{UpdateColley} ) =  update_info( UpdateColley() )

function update_ratings( rule::UpdateColley,
                         input_ratings::RatingsList,
                         input_competitions::DataFrame)
    n = size(input_competitions,1)
    m = length( input_ratings.players )
    I = player_indexes( input_ratings.players )

    # construct Colley matrices and vectors
    C = diagm(0 => 2*ones(Int,m))
    wins = zeros(Int, m)
    losses = zeros(Int, m)
    d = input_competitions
    for i=1:n
        if d[i,Outcome] == 1
            wins[I[d[i,PlayerA]]] += 1
            losses[I[d[i,PlayerB]]] += 1
        elseif d[i,Outcome] == -1
            wins[I[d[i,PlayerB]]] += 1
            losses[I[d[i,PlayerA]]] += 1
        end
        C[ I[d[i,PlayerA]], I[d[i,PlayerB]] ] -= 1
        C[ I[d[i,PlayerB]], I[d[i,PlayerA]] ] -= 1
        C[ I[d[i,PlayerA]], I[d[i,PlayerA]] ] += 1
        C[ I[d[i,PlayerB]], I[d[i,PlayerB]] ] += 1
    end

    # solve Colley's equation
    b = 1 .+ 0.5*(wins .- losses)
    r = C \ b
    ratings = Dict{String, Float64}()
    for player in input_ratings.players
        ratings[player] = r[ I[player] ]
    end

    # output ratings list
    output_ratings = RatingsList(input_ratings.players, ratings )
    return output_ratings
end

# Colley's predictions are deterministic, unless you apply some model after the fact
function predict_outcome(rule::UpdateColley,
                         ratingA::Real, ratingB::Real,
                         factorA::Union{Missing,Real}, factorB::Union{Missing,Real})
    rating_diff = ratingA - ratingB
    s = sign(rating_diff)
    if s==1
        return (1.0, 0.0, 0.0)
    elseif s==-1
        return (0.0, 1.0, 0.0)
    else
        return (0.0, 0.0, 1.0)
    end
end


text:
using LinearAlgebra

function readInput(filename::String)
	x::Vector{Int} = []
    y::Vector{Int} = []
    folds::Vector{Tuple{String, Int}} = []
    open(filename) do file
        for line in readlines(file)
            if isempty(line)
                continue
            end

            if occursin(",", line)
                x_line, y_line = parse.(Int, split(line, ","))
                push!(x, x_line)
                push!(y, y_line)
            else
                if occursin("x", line)
                    x_ = parse.(Int, line[length("fold along x=") + 1:end])
                    push!(folds, ("x", x_))
                else
                    y_ = parse.(Int, line[length("fold along y=") + 1:end])
                    push!(folds, ("y", y_))
                end
            end
        end
    end
    max_x = maximum(x) + 1
    max_y = maximum(y) + 1
    m = falses(max_x, max_y)

    for (xᵢ, yᵢ) in zip(x, y)
        m[xᵢ + 1, yᵢ + 1] = true
    end

    return BitMatrix(transpose(m)), folds
end

function print_matrix(mat)
    show(stdout, "text/plain", mat)
    println()
end

function fold_y(mat, y::Int)
    r = size(mat, 1)
    top_half = mat[1:y, :]
    bot_half = mat[y+2:end, :]
    if y + 1 > r / 2
        top_half[end-size(bot_half, 1)+1:end, :] .|= reverse(bot_half, dims=1)
        return top_half
    else
        bot_half[end-size(top_half, 1)+1:end, :] .|= reverse(top_half, dims=1)
        return reverse(bot_half, dims=1)
    end
end

function fold_x(mat, x::Int)
    c = size(mat, 2)
    left_half = mat[:, 1:x]
    rigt_half = mat[:, x+2:end]
    if x + 1 > c / 2
        left_half[:, end-size(rigt_half, 2)+1:end] .|= reverse(rigt_half, dims=2)
        return left_half
    else
        rigt_half[:, end-size(left_half, 2)+1:end] .|= reverse(left_half, dims=2)
        return reverse(rigt_half, dims=2)
    end
end

function part1(mat, folds)::Int
    dir, n = folds[1]
    if dir == "x"
        return sum(fold_x(mat, n))
    else
        return sum(fold_y(mat, n))
    end
end

function part2(mat, folds)
    for (dir, n) in folds
        if dir == "x"
            mat = fold_x(mat, n)
        else
            mat = fold_y(mat, n)
        end
    end

    return mat
end

mat, folds = readInput("input.txt")

println(part1(mat, folds))
print_matrix(part2(mat, folds))


text:
using BenchmarkTools
using CUDA
using Oceananigans
using Benchmarks

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]
Float_types = [Float32, Float64]
Ns = [32, 64, 128, 256]

# Define benchmarks

SUITE = BenchmarkGroup()

for Arch in Architectures, FT in Float_types, N in Ns
    @info "Setting up benchmark: ($Arch, $FT, $N)..."

    grid = RegularRectilinearGrid(FT, size=(N, N, N), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=Arch(), float_type=FT, grid=grid)

    time_step!(model, 1) # warmup

    benchmark = @benchmarkable begin
        @sync_gpu time_step!($model, 1)
    end samples=10

    SUITE[(Arch, FT, N)] = benchmark
end


text:
export CartesianTrajectory, cartesianTrajectory2dNodes, cartesian2dDensity

"""
    CartesianTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , kmin=(-0.5,-0.5)
                  , kmax=(0.5,0.5)
                  , kargs...)

returns a 2d cartesian trajectory.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* (`kmin=(-0.5,-0.5)`)            - minimum values of the covered k-space (for partial Fourier imaging)
* (`kmax=(-0.5,-0.5)`)            - maximum values of the covered k-space (for partial Fourier imaging)
"""
function CartesianTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , kmin=(-0.5,-0.5)
                  , kmax=(0.5,0.5)
                  , kargs...)
  nodes = cartesian2dNodes(numProfiles, numSamplingPerProfile; kmin=kmin,kmax=kmax)
  times = readoutTimes(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return  Trajectory("Cartesian", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, true, false)
end

function cartesian2dNodes(numProfiles, numSamplingPerProfile
                  ; kmin=(-0.5,-0.5)
                  , kmax=(0.5,0.5)
                  , kargs...)

  nodes = zeros(2,numSamplingPerProfile, numProfiles)
  posX = collect( -ceil(Int64, (numSamplingPerProfile-1)/2.):floor(Int64, (numSamplingPerProfile-1)/2.) ) / numSamplingPerProfile
  posY = collect( -ceil(Int64, (numProfiles-1)/2.):floor(Int64, (numProfiles-1)/2.) ) / numProfiles

  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      nodes[1,k,l] = posX[k]
      nodes[2,k,l] = posY[l]
    end
  end

  # rescale nodes according to the region covered by kmin & kmax
  nodes[1,:,:] .*= (kmax[1]-kmin[1])
  nodes[2,:,:] .*= (kmax[2]-kmin[2])
  # shift nodes to the region defined by kmin & kmax
  nodes[1,:,:] .+= 0.5*(kmin[1]+kmax[1])
  nodes[2,:,:] .+= 0.5*(kmin[2]+kmax[2])

  return reshape(nodes, 2, numSamplingPerProfile*numProfiles)
end

function cartesian2dDensity(numProfiles::Int64, numSamplingPerProfile::Int64)
  density = zeros(numSamplingPerProfile, numProfiles)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      density[k,l] = 1
    end
  end
  return vec(density)
end


text:
@target ptx function kernel_scalaradd(a::CuDeviceArray{Float32}, x)
    i = blockId_x() + (threadId_x()-1) * numBlocks_x()
    a[i] = a[i] + x

    return nothing
end

# TODO: get and compare dim tuple instead of xyz
@target ptx function kernel_lastvalue(a::CuDeviceArray{Float32},
                                      x::CuDeviceArray{Float32})
    i = blockId_x() + (threadId_x()-1) * numBlocks_x()
    max = numBlocks_x() * numThreads_x()
    if i == max
        x[1] = a[i]
    end

    return nothing
end


text:
using SSCF
using Test

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


text:
### INSTALLING PACKAGES

import Pkg
Pkg.update()
Pkg.add("Statistics")
Pkg.add("Distributions")
Pkg.add("Plots")
Pkg.add("StatsPlots")
Pkg.add("LaTeXStrings")

### GET A FEELING

# Create vectors
n = [1 2 3 4 5 6 7 8 9 10]
n = collect(1:10)
n = Vector(undef,10)
n = ones(10)

# Create matrices
A = [1 2 3 ; 4 5 6 ; 7 8 9]
A = Matrix(undef,5,3)
A = zeros(5,3)

# Ouput result
print(n)

# Loops
for i in 1:10
    print(i," ")
end

### APPLY EASY MATH

n2 = n.^2   # . operator is for piece-wise multiplication
n2 = n.*n

A = [1 0.5 0.3; 0.5 1 0.4; 0.3 0.4 1]
B = inv(A)
B = A'

# USING PACKAGES

m = mean(n)
using Statistics
m = mean(n)

### PLOTTING

using Plots

plot(A)

# plotarguments
labels = ["1" "2" "3"]

p1 = plot(A,xlabel="x",ylabel="y",label=labels)
p2 = plot(A,xlabel="x",ylabel="y",seriestype=:scatter,label="")
p3 = plot(A,xlabel="x",ylabel="y",yerror=A*0.1,label="")

# histogram, boxplot
p4 = histogram(A, xlabel="y",ylabel="Frequency",nbins=5, fillalpha=0.5, label="", title = "histogram")
using StatsPlots
p5 = boxplot(A, title = "boxplot",xticks=(1:length(labels),labels), label="")

# subplot
plot(p1,p2,p3,p4,p5)

# Plot sinus
t = collect(0:0.01:1);
T = [0.1 0.2 0.5]
y = sin.(2*pi./T.*t)
plot(t,y,size=(500,100))

using LaTeXStrings   # for latex in figures
plot(t,y,label=[L"y_0=0","y0=0","y_0 = 0"]) # by putting L in front, it's latex


text:
include("registration.jl")

# Classic Control
#=============================================================================#
register("CartPole-v0",
         :CartPoleEnv,
         "/classic_control/CartPole.jl",
	     max_episode_steps=200,
		 reward_threshold=195.0)

register("CartPole-v1",
         :CartPoleEnv,
         "/classic_control/CartPole.jl",
	     max_episode_steps=500,
		 reward_threshold=475.0)

register("Pendulum-v0",
         :PendulumEnv,
         "/classic_control/Pendulum.jl",
		 max_episode_steps=200)

register("MountainCarContinuous-v0",
         :Continuous_MountainCarEnv,
         "/classic_control/Continuous-MountainCar.jl",
		 max_episode_steps=999,
		 reward_threshold=90.0)

register("MountainCar-v0",
		 :MountainCarEnv,
		 "/classic_control/MountainCar.jl",
		 max_episode_steps=200,
		 reward_threshold=110.0)

register("Acrobot-v0",
		 :AcrobotEnv,
		 "/classic_control/Acrobot.jl",
		 max_episode_steps=500,
		 reward_threshold=100.0)

# Algorithmic environments
#=============================================================================#

register("Copy-v0",
		 :CopyEnv,
		 "/algorithmic/Copy.jl",
		 max_episode_steps=200,
		 reward_threshold=25.0)

register("RepeatCopy-v0",
		 :RepeatCopyEnv,
		 "/algorithmic/RepeatCopy.jl",
		 max_episode_steps=200,
		 reward_threshold=75.0)

register("ReversedAddition-v0",
		 :ReversedAdditionEnv,
		 "/algorithmic/ReversedAddition.jl",
		 max_episode_steps=200,
		 reward_threshold=25.0)

register("DuplicatedInput-v0",
		 :DuplicatedInputEnv,
		 "/algorithmic/DuplicatedInput.jl",
		 max_episode_steps=200,
		 reward_threshold=25.0)

register("Reverse-v0",
		 :ReverseEnv,
		 "/algorithmic/Reverse.jl",
		 max_episode_steps=200,
		 reward_threshold=25.0)

# Atari environments
#=============================================================================#

games_list = ["adventure", "air_raid", "alien", "amidar", "assault", "asterisk", "atlantis",
			  "bank_heist", "battle_zone", "beam_rider", "berzerk", "bowling", "boxing", "breakout", "carnival",
			  "centipede", "chopper_command", "crazy_climber", "defender", "demon_attack", "double_dunk",
			  "elevator_action", "enduro", "fishing_derby", "freeway", "frostbite", "gopher", "gravitar",
			  "hero", "ice_hockey", "jamesbond", "journey_escape", "kangaroo", "krull", "kung_fu_master",
			  "montezuma_revenge", "ms_pacman", "name_this_game", "pheonix", "pitfall", "pong", "pooyan",
			  "private_eye", "qbert", "riverraid", "road_runner", "robotank", "seaquest", "skiing",
			  "solaris", "space_invaders", "star_gunner", "tennis", "time_pilot", "tutankham",
			  "up_n_down", "venture", "video_pinball", "wizard_of_wor", "yars_revenge", "zaxxon"]

for game ∈ games_list
	for obs_type ∈ (:ram, :color, :grey)
		# space_invaders should yield SpaceInvaders-v0 and SpaceInvaders-ram-v0
		name = game |> (g) -> split(g, '_') .|> titlecase |> join
		obs_type == :ram && (name = "$(name)-ram")
		obs_type == :grey && (name = "$(name)-grey")

		# ElevatorAction-ram-v0 seems to yield slightly
		# non-deterministic observations about 10% of the time. We
		# should track this down eventually, but for now we just
		# mark it as nondeterministic.
		nondeterministic = game=="elevator_action" && obs_type==:ram

		register("$(name)-v0",
				 :AtariEnv,
				 "/atari/Atari.jl",
				 max_episode_steps = 10000,
				 nondeterministic = nondeterministic,
				 kwargs = Dict(
				 	:game_path => game,
				 	:obs_type  => obs_type,
				 	:repeat_action_probability => 0.25f0
				 	)
			)

		register("$(name)-v4",
				 :AtariEnv,
			 	 "/atari/Atari.jl",
				 max_episode_steps = 100000,
				 nondeterministic = nondeterministic,
				 kwargs = Dict(
				 	:game_path => game,
				 	:obs_type  => obs_type,
				 	)
			)

		frameskip = game=="space_invaders" ? 3 : 4

		register("$(name)Deterministic-v0",
				 :AtariEnv,
				 "/atari/Atari.jl",
				 max_episode_steps = 100000,
				 nondeterministic = nondeterministic,
				 kwargs = Dict(
				 	:game_path => game,
					:obs_type  => obs_type,
					:frameskip => frameskip
				 )
			)

		register("$(name)Deterministic-v4",
				 :AtariEnv,
				 "/atari/Atari.jl",
				 max_episode_steps = 100000,
				 nondeterministic = nondeterministic,
				 kwargs = Dict(
				 	:game_path => game,
					:obs_type  => obs_type,
					:frameskip => frameskip
				 )
			)

		register("$(name)NoFrameskip-v0",
				 :AtariEnv,
				 "/atari/Atari.jl",
				 max_episode_steps = frameskip * 100000,
				 nondeterministic = nondeterministic,
				 kwargs = Dict(
				 	:game_path => game,
					:obs_type  => obs_type,
					:frameskip => 1,
					:repeat_action_probability => 25f-2
				 )
			)

		# No frameskip. (Atari has no entropy source, so these are
	    # deterministic environments.)
		register("$(name)NoFrameskip-v4",
				 :AtariEnv,
				 "/atari/Atari.jl",
				 max_episode_steps = frameskip * 100000,
				 nondeterministic = nondeterministic,
				 kwargs = Dict(
				 	:game_path => game,
					:obs_type  => obs_type,
					:frameskip => 1
				 )
			)
	end
end


text:
module EnsembleSimulations

export ensemble_column_model_simulation

using DataDeps

using ..Observations: SyntheticObservations, batch, tupleit

using Oceananigans
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.Architectures: arch_array

function ensemble_column_model_simulation(observations;
                                          closure,
                                          Nensemble,
                                          Δt = 1.0,
                                          architecture = CPU(),
                                          tracers = :b,
                                          buoyancy = BuoyancyTracer(),
                                          kwargs...)

    observations = batch(observations)
    Nbatch = length(observations)

    # Assuming all observations are on similar grids...
    Nz = first(observations).grid.Nz
    Hz = first(observations).grid.Hz
    Lz = first(observations).grid.Lz

    column_ensemble_size = ColumnEnsembleSize(Nz=Nz, ensemble=(Nensemble, Nbatch))
    column_ensemble_halo_size = ColumnEnsembleSize(Nz=0, Hz=Hz)

    grid = RectilinearGrid(architecture,
                           size = column_ensemble_size,
                           halo = column_ensemble_halo_size,
                           topology = (Flat, Flat, Bounded),
                           z = (-Lz, 0))

    coriolis_ensemble = [FPlane(f=observations[j].metadata.coriolis.f) for i = 1:Nensemble, j=1:Nbatch]
    coriolis_ensemble = arch_array(architecture, coriolis_ensemble)

    closure_ensemble = [deepcopy(closure) for i = 1:Nensemble, j=1:Nbatch]
    closure_ensemble = arch_array(architecture, closure_ensemble)

    momentum_boundary_conditions =
        (; u = FieldBoundaryConditions(top = FluxBoundaryCondition(zeros(grid, Nensemble, Nbatch))))

    ensemble_tracer_bcs() = FieldBoundaryConditions(top = FluxBoundaryCondition(zeros(grid, Nensemble, Nbatch)),
                                                    bottom = GradientBoundaryCondition(zeros(grid, Nensemble, Nbatch)))

    tracers = tupleit(tracers)
    tracer_boundary_conditions = NamedTuple(name => ensemble_tracer_bcs() for name in tracers if name != :e)

    boundary_conditions = merge(momentum_boundary_conditions, tracer_boundary_conditions)

    ensemble_model = HydrostaticFreeSurfaceModel(; grid, tracers, buoyancy, boundary_conditions,
                                                 coriolis = coriolis_ensemble,
                                                 closure = closure_ensemble,
                                                 kwargs...)

    ensemble_simulation = Simulation(ensemble_model; Δt, stop_time=first(observations).times[end])

    return ensemble_simulation
end

end # module


text:
module Units

export SecYear, second, yr, kyr, myr, meter, km, km³

#####
##### Convenient definitions
#####

"""
    second

A `Float64` constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Δt = 1second`.
"""
const second = 1.0

"""
    SecYear

A `Float64` constant equal to 365.25 * 24 * 3600.
"""
const SecYear = 365.25*24*3600*second

"""
    year

A `Float64` constant equal to 365.25 * 24 * 3600 seconds. Useful for increasing the clarity of scripts, e.g. `Δt = 7yr`.
"""
const yr = 365.25*24*3600*second

"""
    kyr

A `Float64` constant equal to 1000`year`. Useful for increasing the clarity of scripts, e.g. `Δt = 10kyr`.
"""
const kyr = 1000yr

"""
    myr

A `Float64` constant equal to 10^6`year`. Useful for increasing the clarity of scripts, e.g. `Δt = 1myr`.
"""
const myr = 1000kyr

"""
    meter

A `Float64` constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `L = 1meter`.
"""
const meter = 1.0

"""
    km

A `Float64` constant equal to 1000`m`. Useful for increasing the clarity of scripts, e.g. `L = 50km`.
"""
const km = 1000meter

"""
    km³

A `Float64` constant equal to 1000`m³`. Useful for increasing the clarity of scripts, e.g. `Vol = 50km³`.
"""
const km³ = 1e9*meter

end


text:
module Batsrus
# BATSRUS data reader and analyzer.
#
# Hongyang Zhou, hyzhou@umich.edu

using Printf, Requires

export Data

"Type for the file information."
struct FileList
   "filename"
   name::String
   "file type"
   type::String
   "directory"
   dir::String
   "file size"
   bytes::Int
   "number of snapshots"
   npictinfiles::Int
end

"Primary data storage type"
struct Data{T<:AbstractFloat}
   "header information"
   head::NamedTuple
   "grid"
   x::Array{T}
   "variables"
   w::Array{T}
   "file information"
   list::FileList
end

include("io.jl")
include("select.jl")
include("vtk.jl")

include("unit/UnitfulBatsrus.jl")
using UnitfulRecipes, .UnitfulBatsrus

function __init__()
   @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
      include("plot/pyplot.jl")
   end
   @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
      include("plot/plots.jl")
   end
end

end


text:
# example 3: read all at once with readall (or readfasta)
#            default return type (ASCIIString)

module FastaIOExample3
using FastaIO

function read_fasta_file(filename::String)
    # note: FastaReader can take a readable IO descriptor as
    #       argument instead of a filename; it also supports
    #       do-notation:
    #
    #       FastaReader(filename) do fr
    #         # do something
    #       end
    fr = FastaReader(filename)

    out = readall(fr)

    # alternatively:
    # out = readfasta(filename)

    n = 0
    for (desc, seq) in out
        n += 1
        println("num=$n desc=$desc seq=$seq")
    end
    # here of course n == fr.num_parsed == length(out)
    println("read $n entries")
end

examples_dir = dirname(Base.source_path())
read_fasta_file(joinpath(examples_dir, "example.fasta.gz"))
end


text:
# Loading and saving utilities

function savemodel(model, path)
    @save path model = cpu(model)
end

function loadmodel(path)
    @load path model
    return model
end

function saveweights(model, path)
    @save path weights = params(cpu(model))
end

function loadweights(path)
    @save path weights
    return weights
end


text:
using RDatasets, DataFrames, Gadfly

set_default_plot_size(6inch, 3inch)

df = dataset("plm", "Cigar")
df = df[df[:State] .== 1, :]
ymin = df[:Sales] .- 20*0.23603334566204692
ymax = df[:Sales] .+ 20*0.34651701419196046

plot(df, x="Year", y="Sales", ymin=ymin, ymax=ymax,
     Geom.bar, Geom.errorbar)


text:
# Given a function f of z and a complex array z, evaluate the gradient wrt z at z
function gradient(f, z::AbstractArray{Complex{V}},
            cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, z)) where {T, V}
    #CHK && checktag(T, f, z)
    checktag(T, f, z)
    if chunksize(cfg) == length(z)
      vector_mode_gradient(f, z, cfg)
    else
      chunk_mode_gradient(f, z, cfg)
    end
end

# Given a complex dual, return dz and dzstar as Partial types. Ugly way of doing
# this
@inline function extract_gradient!(dz,dzstar,d::ComplexDual{T,V,N}) where {T,V,N}
    tmp, tmpstar = dz_partials(d)
    dz .= tmp
    dzstar .= tmpstar
    return dz,dzstar
end

# This is for making sure that the desired gradient has the same tag as the dual
@inline extract_gradient!(::Type{T},dz,dzstar,d::ComplexDual{T,V,N}) where {T,V,N} =
          extract_gradient!(dz,dzstar,d)

function vector_mode_gradient(f::F, z, cfg::ComplexGradientConfig{T},
                evalfcn=ForwardDiff.vector_mode_dual_eval) where {F,T}
    ydual = evalfcn(f, z, cfg)
    dz = similar(z, Complex{valtype(ydual)})
    dzstar = similar(z, Complex{valtype(ydual)})
    return extract_gradient!(T, dz, dzstar,ydual)
end

function extract_gradient_chunk!(::Type{T}, dz, dzstar, d::ComplexDual{T,V,N}, index, chunksize) where {T,V,N}
    offset = index - 1
    for i in 1:chunksize
        tmp, tmpstar = dz_partials(T, d, i)
        dz[i + offset] = tmp
        dzstar[i + offset] = tmpstar
    end
    return dz, dzstar
end

function chunk_mode_gradient_expr(work_array_definition::Expr,result_definition::Expr,
                                          result::Expr,compute_ydual::Expr,
                                          error_definition::Expr,
                                          seed_chunk_definition::Expr, unseed_definition::Expr,
                                          extract_definition::Expr,y_definition::Expr)
  return quote

    @assert length(z) >= N "chunk size cannot be greater than length(z) ($(N) > $(length(z)))"

    # precalculate loop bounds
    zlen = length(z)
    remainder = zlen % N
    lastchunksize = ifelse(remainder == 0, N, remainder)
    lastchunkindex = zlen - lastchunksize + 1
    middlechunks = 2:div(zlen - lastchunksize, N)

    # seed work vectors
    $(work_array_definition)
    rseeds = cfg.rseeds
    iseeds = cfg.iseeds

    index, chunksize = 1, N
    $(seed_chunk_definition)
    $(compute_ydual)
    $(error_definition)
    $(result_definition)

    $(extract_definition)
    $(unseed_definition)

    for c in middlechunks
      i = ((c - 1) * N + 1)
      index, chunksize = i, N
      $(seed_chunk_definition)
      $(compute_ydual)
      $(extract_definition)
      $(unseed_definition)
    end
    index, chunksize = lastchunkindex, lastchunksize
    $(seed_chunk_definition)
    $(compute_ydual)
    $(extract_definition)

    $(y_definition)

    return $(result)
  end
end

@eval function chunk_mode_gradient(f::F, z, cfg::ComplexGradientConfig{T,V,N}) where {F,T,V,N}
    $(chunk_mode_gradient_expr(quote
                                  zdual = cfg.duals
                                  seed!(zdual, z)
                                end,
                                quote
                                  dz = similar(z, Complex{valtype(ydual)})
                                  dzstar = similar(z, Complex{valtype(ydual)})
                                end,
                                :(dz, dzstar),
                                :(ydual = f(zdual)),
                                :(),
                                :(seed!(zdual, z, index, rseeds, iseeds, chunksize)),
                                :(seed!(zdual, z, index)),
                                :(extract_gradient_chunk!(T, dz, dzstar, ydual, index, chunksize)),
                                :()))
end


text:
##############################################################################################################################
# Some core functions needed to make thinks work smoothly :)
##############################################################################################################################

# Data member access functions (useful just in case member names change ...)
# Observation: These functions should be used whenever access to the internals of objects are needed
"""
Get a copy of the `x` field of types derived from `AbstractCell`
"""
getx(ac::T where T<:AbstractCell) = deepcopy(ac.x)

"""
Reference the `x` field of types derived from `AbstractCell`
"""
getx!(ac::T where T<:AbstractCell) = ac.x

"""
Get a copy of the `y` field of types derived from `AbstractCell`
"""
gety(ac::T where T<:AbstractCell) = deepcopy(ac.y)

"""
Reference the `y` field of types derived from `AbstractCell`
"""
gety!(ac::T where T<:AbstractCell) = ac.y

"""
Get a copy of the function from the `f` field of `AbstractCell`
"""
getf(c::T where T<:AbstractCell) = deepcopy(c.f)

"""
Reference the function from the `f` field of `AbstractCell`
"""
getf!(c::T where T<:AbstractCell) = c.f

# Number of variables for data vectors, matrices
nvars(x::Tuple{T} where T<:AbstractArray) = nvars(x[1])
nvars(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray) = nvars(x[1])
nvars(x::T where T<:AbstractVector) = 1
nvars(x::T where T<:AbstractMatrix) = size(x,1)

"""
Return the number of variables from a `DataCell` (e.g. the dimensionality of the dataset)
"""
nvars(x::T where T<:CellData) = nvars(getx!(x))

# Return tuple of data matrix, targets
"""
Return a Tuple containing the data and targets from a `DataCell`
"""
strip(x::T where T<:CellDataU) = getx!(x)
strip(x::T where T<:CellData) = (getx!(x), gety!(x))

# Size related functions for cell-type objects
size(::T where T<:Void) = 0
size(::T where T<:Void, i::Int) = 0
size(ac::T where T<:AbstractCell) = size(getx!(ac))							# Size of an abstract cells' data
size(ac::T where T<:AbstractCell, i::Int) = size(getx!(ac),i)						# Size along a dimension of an abstract cells' data
ndims(dc::T where T<:CellData) = ndims(getx!(dc))							# Number of dimensions

# Function needed to properly find NaNs in heterogeneous arrays
isnan(::T where T<:Void) = Bool(false)
isnan(x::T where T<:AbstractString) = Bool(false)
isnan(x::Array{T} where T<:AbstractString) = return BitArray(falses(size(x)))
isnan(x::DataArray{T} where T<:AbstractString) = return BitArray(falses(size(x)))
isnan(x::T where T<:Array) = isnan.(x)

# Function to check for voids
isvoid(::Void) = true
isvoid(x) = (x isa Void)

# Class related functions for DataCells (e.g. class sizes - normalized or not, class names)
"""
Returns a `SortedDict` with the label and associated sample count.
"""
classsizes(c::T where T<:CellDataU) = begin
	return SortedDict(Dict("unlabeled" => nobs(c)))
end

classsizes(c::T where T<:CellDataL) = begin							# Class sizes for DataCells; of the form Dict(clas sname => class size)
	labels = gety!(c)
	return SortedDict(countmap(labels))
end

classsizes(c::T where T<:CellDataLL) = begin 							# Class sizes for DataCells with multiple labels; of the form Vector{Dict(...)}
	labels = gety(c)
	return [SortedDict(countmap(labels[i,:])) for i = 1:nvars(gety!(c)) ]
end

"""
Normalized count version of `countmap`.
"""
countmapn(x::T where T<:AbstractVector) = begin
	cm = countmap(x)
	s = eltype(1.0)(sum(values(cm)))
	out = Dict{eltype(x),eltype(1.0)}()
	for k in keys(cm)
		push!(out,k=>eltype(1.0)(cm[k])/s)
	end
	return SortedDict(out)
end

countmapn(::T where T<:Void) = SortedDict(Dict())

"""
Returns a `SortedDict` with the label and associated normalized sample count.
"""
nclasssizes(c::T where T<:CellDataU) = SortedDict(Dict("unlabeled"=>1))				# Normalized class sizes for DataCells

nclasssizes(c::T where T<:CellDataL) = begin 							# Normalized class sizes for DataCells
	return countmapn(gety!(c))
end

nclasssizes(c::T where T<:CellDataLL) = begin 							# Normalized class sizes for DataCells with multiple labels
	normcmv = Vector{SortedDict}()
	for i in 1 : size(gety!(c),1)
		push!(normcmv, countmapn(gety!(c)[i,:]))
	end
	return normcmv
end

# Functions that return class names
"""
Return the class labels from a `DataCell`.
"""
classnames(c::T where T<:CellDataU) = []
classnames(c::T where T<:CellDataL) = begin
	try sort(unique(gety!(c)))
	catch
		unique(gety!(c))

	end
end
classnames(c::T where T<:CellDataLL) = begin
	try
		[sort(unique(gety!(c)[i,:])) for i in 1:nvars(gety!(c))]
	catch

		[unique(gety!(c)[i,:]) for i in 1:nvars(gety!(c))]
	end
end

# Functions that return class numbers
"""
Return the number of classes present in a `DataCell`.
"""
nclass(x::T where T<:CellDataU) = 0
nclass(x::T where T<:CellDataL) = length(classnames(x))
nclass(x::T where T<:CellDataLL) = [length(c) for c in classnames(x)]

# Other various deletion functions
deleteat!(::Void, args...)::Void = nothing

deleteat(x::T where T<:AbstractMatrix{N} where N, rows) = begin
	m::Int = size(x,1)
	n::Int = size(x,2)
	y::Vector{N} = squeeze(reshape(x,length(x),1),2)
	pos=Vector{Int}(length(rows)*n)

	k=1
	for i in rows
		@simd for j in 1:n
			@inbounds pos[k]=i+(j-1)*m
			k+=1
		end
	end
	pos = reshape(pos, length(pos))
	y = y[setdiff(1:length(y),pos)]

	return reshape(y, round(Integer,length(y)/n), n)
end

# Function that counts the layers in an AbstractCell
"""
	countlayers(x::Tuple)

Keeps track of the level of nesting in pipes.
"""
countlayers(x::T where T<:PTuple{AbstractCell}) = maximum(xi.layer for xi in x)

# Small function that builds the titles out of keywork argumens
kwtitle(root::String, kwargs)::String = begin
	title = root
	if isempty(kwargs)
		return title
	end
	title *= ": "
	n = length(kwargs)
	for k in 1:length(kwargs)
		title*=string(kwargs[k][1])*"="*string(kwargs[k][2])*(k==n ? "" : " ")
	end
	return title
end

# Positional argument selectors
idx(args...) = FunctionCell(getindex, args, "Index selector (args=$args)")
idx(f::Function, args...) = FunctionCell((x,args...)->f(getindex(x,args...)), args, "Index selector+function (args=$args)")


text:
include("Shutter.jl")
include("ErrorMeasures.jl")
include("CoilSensitivity.jl")
include("Regridding.jl")


text:
using DiffEqProblemLibrary, Sundials, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

prob = prob_ode_2Dlinear
integrator = init(prob,CVODE_BDF())
step!(integrator)
integrator(integrator.t,Val{1})

for i in integrator
    @show i.t
end

integrator.t == 1.0
step!(integrator)
integrator.t > 1.0


text:
using Cxx
include(Pkg.dir("Cxx","test","llvmincludes.jl"))
cxx"""
#include "llvm/Option/Option.h"
#include "clang/Tooling/JSONCompilationDatabase.h"
#include "clang/Driver/Compilation.h"
#include "clang/Driver/Driver.h"
#include "clang/Driver/Job.h"
#include "clang/Driver/Tool.h"
#include "clang/Basic/DiagnosticIDs.h"
#include "clang/Basic/DiagnosticOptions.h"
#include "clang/Frontend/TextDiagnosticPrinter.h"
#include "clang/Frontend/FrontendDiagnostic.h"
using namespace llvm;
using namespace llvm::opt;
using namespace clang;
using namespace clang::tooling;

/// \brief Retrieves the clang CC1 specific flags out of the compilation's jobs.
///
/// Returns NULL on error.
static const llvm::opt::ArgStringList *getCC1Arguments(
    clang::DiagnosticsEngine *Diagnostics,
    clang::driver::Compilation *Compilation) {
  // We expect to get back exactly one Command job, if we didn't something
  // failed. Extract that job from the Compilation.
  const clang::driver::JobList &Jobs = Compilation->getJobs();
  if (Jobs.size() != 1 || !isa<clang::driver::Command>(*Jobs.begin())) {
    SmallString<256> error_msg;
    llvm::raw_svector_ostream error_stream(error_msg);
    Jobs.Print(error_stream, "; ", true);
    Diagnostics->Report(clang::diag::err_fe_expected_compiler_job)
        << error_stream.str();
    return nullptr;
  }

  // The one job we find should be to invoke clang again.
  const clang::driver::Command &Cmd =
      cast<clang::driver::Command>(*Jobs.begin());
  if (StringRef(Cmd.getCreator().getName()) != "clang") {
    Diagnostics->Report(clang::diag::err_fe_expected_clang_command);
    return nullptr;
  }

  return &Cmd.getArguments();
}

/// \brief Returns a clang build invocation initialized from the CC1 flags.
clang::CompilerInvocation *newInvocation(
    clang::DiagnosticsEngine *Diagnostics,
    const llvm::opt::ArgStringList &CC1Args) {
  assert(!CC1Args.empty() && "Must at least contain the program name!");
  clang::CompilerInvocation *Invocation = new clang::CompilerInvocation;
  clang::CompilerInvocation::CreateFromArgs(
      *Invocation, CC1Args.data() + 1, CC1Args.data() + CC1Args.size(),
      *Diagnostics);
  Invocation->getFrontendOpts().DisableFree = false;
  Invocation->getCodeGenOpts().DisableFree = false;
  Invocation->getDependencyOutputOpts() = DependencyOutputOptions();
  return Invocation;
}
"""
invocation = icxx"""
std::string error;
auto cdb = clang::tooling::JSONCompilationDatabase::loadFromFile("/home/kfischer/julia/compile_commands.json", error);
CompilerInvocation Inv;
std::vector<const char*> Argv;
CompileCommand TheCommand = cdb->getCompileCommands("/home/kfischer/julia/src/gf.c")[0];
for (const std::string &Str : TheCommand.CommandLine)
  Argv.push_back(Str.c_str());
const char *const BinaryName = Argv[0];
IntrusiveRefCntPtr<DiagnosticOptions> DiagOpts = new DiagnosticOptions();
TextDiagnosticPrinter DiagnosticPrinter(
    llvm::errs(), &*DiagOpts);
DiagnosticsEngine Diagnostics(
    IntrusiveRefCntPtr<clang::DiagnosticIDs>(new DiagnosticIDs()), &*DiagOpts,
    &DiagnosticPrinter, false);

const std::unique_ptr<clang::driver::Driver> Driver(new clang::driver::Driver(
    BinaryName, llvm::sys::getDefaultTargetTriple(), Diagnostics));
// Since the input might only be virtual, don't check whether it exists.
Driver->setCheckInputsExist(false);
const std::unique_ptr<clang::driver::Compilation> Compilation(
    Driver->BuildCompilation(llvm::makeArrayRef(Argv)));
const llvm::opt::ArgStringList *const CC1Args = getCC1Arguments(
    &Diagnostics, Compilation.get());
if (!CC1Args) {
  assert(false);
}
return newInvocation(&Diagnostics, *CC1Args);
"""

RemoteClang = Cxx.setup_instance_from_inovcation(invocation);
Cxx.initialize_instance!(RemoteClang; register_boot = false)
push!(Cxx.active_instances, RemoteClang)
RemoteClang = Cxx.CxxInstance{length(Cxx.active_instances)}()

nothing


text:
include("poly_factorization_project.jl")
println("Welcome to the Great Polynomial Factoriser")
println("")
println("")
println("----- Example Polynomials -----")
x = x_poly()
p1 = 2x^3 + 4x^2 - 3x
p2 = -6x^4 - 4x^2 - 3x + 3
println(p1)
println(p2)
println("")
println("mod 7")
x = x_poly(7)
p1 = 2x^5 + 5x - 3
p2 = -6x^4 - 4x^2 - 3x + 3
println(p1)
println(p2)
println("")
println("mod 17")
x = x_poly(17)
p1 = 16x^3 - 3
p2 = -6x^8 - 4x^2 - 3x + 1
println(p1)
println(p2)
println("")
println("mod 101")
x = x_poly(101)
a = 99x^8+ 50x^2 - 3
b = 3x^7+ 5x^7 + 10x^2 + 60x
@show a
@show b

println("")

println("----- Basic Operations -----")
@show a+b
@show a-b
@show a*b
@show a÷b
@show a^3
println("")
println("d/dx(", "$a) = ", derivative(a))
println("")
println("-----Factorisation Example------")
p = (7x^3 + 2x^2 + 8x + 1)*(x^2+x+1)

println("The Factorisation of: ", p," is:")
println("")
factorization = factor(p,true)
println("Here it is:")
pretty_factor(factorization)
pr = expand_factorization(factorization)
println("")
println("Reconstructing: ", pr)

println(" ")
println("---- Extended Euclid Algorithm on x^2 + 50 and 10x^3+10 (mod 101) ----")
egcd = extended_euclid_alg(x^2 + 50,10x^3+10)
pretty_print_egcd((((x^2 + 50)*(10x^3+10)).polynomial,b.polynomial), egcd)

println("")
println("----- Bechmarking -----")
println("Time taken to raise $a to the 1000th power:")
@time l = a^1000
println("")
multbenchZ()
println("")
multbenchZp()
println("")
powbench()
println("")

println(":)")