text:
module trunc_cyl_shell_examples
using FinEtools
using FinEtools.AlgoDeforLinearModule
using LinearAlgebra
using Arpack

function trunc_cyl_shell()
    println("""
    Vibration modes of truncated cylindrical shell.
    """)

    # t0 = time()

    E = 205000*phun("MPa");# Young's modulus
    nu = 0.3;# Poisson ratio
    rho = 7850*phun("KG*M^-3");# mass density
    OmegaShift = (2*pi*100) ^ 2; # to resolve rigid body modes
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 5; nl  = 12; nc = 40;
    tolerance = h/nh/100;
    neigvs = 20;

    MR = DeforModelRed3D
    fens,fes  = H8block(h,l,2*pi,nh,nl,nc)
    # Shape into a cylinder
    R = zeros(3, 3)
    for i = 1:count(fens)
        x, y, z = fens.xyz[i,:];
        rotmat3!(R, [0, z, 0])
        Q = [cos(psi*pi/180) sin(psi*pi/180) 0;
        -sin(psi*pi/180) cos(psi*pi/180) 0;
        0 0 1]
        fens.xyz[i,:] = reshape([x+Rmed-h/2, y-l/2, 0], 1, 3)*Q*R;
    end
    candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], thickness = h/1000)
    fens,fes = mergenodes(fens, fes,  tolerance, candidates);

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material=MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)), material)
    femm = associategeometry!(femm, geom)
    K =stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
    M =mass(femm, geom, u)

    # eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors
    # v (only if ritzvec=true), the number of converged eigenvalues nconv, the number
    # of iterations niter and the number of matrix vector multiplications nmult, as
    # well as the final residual vector resid.

    if true
        d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
        d = d .- OmegaShift;
        fs = real(sqrt.(complex(d)))/(2*pi)
        println("Eigenvalues: $fs [Hz]")
        mode = 7
        scattersysvec!(u, v[:,mode])
        File =  "trunc_cyl_shell.vtk"
        vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
        @async run(`"paraview.exe" $File`)
    end

    if true
        solver = AlgoDeforLinearModule.ssit
        v0 = rand(size(K,1), 2*neigvs)
        tol = 1.0e-2
        maxiter = 20
        lamb, v, nconv, niter, nmult, lamberr = solver(K + OmegaShift.*M, M; nev=neigvs, v0=v0, tol=tol, maxiter=maxiter)
        if nconv < neigvs
            println("NOT converged")
        end
        lamb = lamb .- OmegaShift;
        fs = real(sqrt.(complex(lamb)))/(2*pi)
        println("Eigenvalues: $fs [Hz]")
        println("Eigenvalue errors: $lamberr [ND]")
        mode = 7
        scattersysvec!(u, v[:,mode])
        File =  "trunc_cyl_shell.vtk"
        vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
        @async run(`"paraview.exe" $File`)
    end

    true

end # trunc_cyl_shell

function trunc_cyl_shell_nas()
    println("""
    Vibration modes of truncated cylindrical shell. NASTRAN input file.
    """)

    # t0 = time()

    E = 205000*phun("MPa");# Young's modulus
    nu = 0.3;# Poisson ratio
    rho = 7850*phun("KG*M^-3");# mass density
    OmegaShift = (2*pi*100) ^ 2; # to resolve rigid body modes
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 5; nl  = 12; nc = 40;
    tolerance = h/nh/100;
    neigvs = 20;

    MR = DeforModelRed3D
    output = import_NASTRAN("trunc_cyl_shell_2.nas")
    fens, fes = output["fens"], output["fesets"][1]

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    femm = associategeometry!(femm, geom)
    K = stiffness(femm, geom, u)
    M = mass(femm, geom, u)

    # eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors
    # v (only if ritzvec=true), the number of converged eigenvalues nconv, the number
    # of iterations niter and the number of matrix vector multiplications nmult, as
    # well as the final residual vector resid.

    if true
        d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
        d = d .- OmegaShift;
        fs = real(sqrt.(complex(d)))/(2*pi)
        println("Eigenvalues: $fs [Hz]")
        mode = 7
        scattersysvec!(u, v[:,mode])
        File =  "trunc_cyl_shell_nas.vtk"
        vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])
        @async run(`"paraview.exe" $File`)
    end

    true

end # trunc_cyl_shell

function allrun()
    println("#####################################################")
    println("# trunc_cyl_shell ")
    trunc_cyl_shell()
    println("#####################################################")
    println("# trunc_cyl_shell_nas ")
    trunc_cyl_shell_nas()
    return true
end # function allrun

end # module trunc_cyl_shell_examples


text:
# For convenience
typealias TupleN{T,N} NTuple{N,T}

# Cast any Tuple to an TupleN{T}
convert_ntuple{T}(::Type{T},d::T) = T # For zero-dimensional arrays
convert_ntuple{N,T}(::Type{T},d::NTuple{N,T}) = d
@generated function convert_ntuple{N,T}(::Type{T},d::NTuple{N})
    exprs = ntuple(i -> :(convert(T, d[$i])), Val{N})
    return Expr(:tuple, exprs...)
end

# Base gives up on tuples for promote_eltype... (TODO can we improve Base?)
@generated function promote_tuple_eltype{T <: Tuple}(::Union{T,Type{T}})
    t = Union{}
    for i = 1:length(T.parameters)
        tmp = T.parameters[i]
        if tmp <: Vararg
            tmp = tmp.parameters[1]
        end
        t = promote_type(t, tmp)
    end
    return quote
        $(Expr(:meta,:pure))
        $t
    end
end

# some convenience functions for non-static arrays, generators, etc...
@inline convert{T}(::Type{Tuple}, a::AbstractArray{T}) = (a...)::Tuple{Vararg{T}}
@inline function convert{N,T}(::Type{NTuple{N}}, a::AbstractArray{T})
    @boundscheck if length(a) != N
        error("Array of length $(length(a)) cannot be converted to a $N-tuple")
    end

    @inbounds return ntuple(i -> a[i], Val{N})
end

@inline function convert{N,T1,T2}(::Type{NTuple{N,T1}}, a::AbstractArray{T2})
    @boundscheck if length(a) != N
        error("Array of length $(length(a)) cannot be converted to a $N-tuple")
    end

    @inbounds return ntuple(i -> convert(T1,a[i]), Val{N})
end

# TODO try and make this generate fast code
@inline convert(::Type{Tuple}, g::Base.Generator) = (g...)
@inline function convert{N}(::Type{NTuple{N}}, g::Base.Generator)
    @boundscheck if length(g.iter) != N
        error("Array of length $(length(a)) cannot be converted to a $N-tuple")
    end

    @inbounds return ntuple(i -> g.f(g.iter[i]), Val{N})
end

#=
@generated function convert{N}(::Type{NTuple{N}}, g::Base.Generator)
    exprs = [:(g.f(g.iter[$j])) for j=1:N]
    return quote
        @boundscheck if length(g.iter) != N
            error("Array of length $(length(a)) cannot be converted to a $N-tuple")
        end

        @inbounds return $(Expr(:tuple, exprs...))
    end
end=#

ext: jl


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

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.
=#

# __openBF__ is a computational library written
# in [Julia](julialang.org) for blood flow simulations
# in straight elastic arteries. Navier-Stokes equations are reduced
# to 1D form and solved along the longitudinal direction. Numerical
# solution is achieved via a first-order finite volume solver
# based on the Godunov' method. For an usage example refer to
# [main.jl](main.html) file and its documentation.
module openBF

# ### Import openBF' files

# OpenBF' own types are contained in [BTypes.jl](BTypes.html) where
# data structures for vessels, blood, and numerical scheme
# are specified.
using BTypes

# Data structures and output files are initialised at the beginning
# of each simulation. [initialise.jl](initialise.html) contains all
# the functions needed to read model description and the
# global settings.
include("initialise.jl")

# Inlet and outlets boundary conditions are applied at each time step,
# and vessels extremity nodes are updated with functions defined in
# [boundary_conditions.jl](boundary_conditions.html).
include("boundary_conditions.jl")

# Interface problems (junctions and bifurcations) are solved via the
# method of characteristics. The function in
# [junctions.jl](junctions.html) discriminates between
# [conjunctions](conjunctions.html) and
# [bifurcations](bifurcations.jl), and calls the appropriate solver.
include("junctions.jl")
include("conjunctions.jl")
include("bifurcations.jl")

# [godunov.jl](godunov.html) contains Godunov' method functions.
include("godunov.jl")

# [MUSCL.jl](MUSCL.html) is where the MUSCL method is implemented.
include("MUSCL.jl")

# Conversions from pressure to cross-sectional area (and vice versa),
# or from cross-sectional area to wave speed (and vice versa) are
# handled by functions in [converter.jl](converter.html); Riemann
# invariants are computed in the same location as well.
include("converter.jl")

# Input and output writing functions are stored defined in
# [IOutils.jl](IOutils.html).
include("IOutils.jl")

# Convergence check functions are all defined in [check_convergence.jl](check_convergence.html).
include("check_convergence.jl")

include("anastomosis.jl")

end


text:
"""
    WeightedParticleRepresentation{S} <: ParticleRepresentation{S} <: AbstractFilterRepresentation{S}

Abstract type for the representation of the conditional distribution over the hidden state by weighted particles (samples).
"""
abstract type WeightedParticleRepresentation{S} <: ParticleRepresentation{S} end

"""
    get_weight(ensemble::WeightedParticleRepresentation, i)

Return the weight of the particle with index `i` from `ensemble`.
"""
function get_weight(ensemble::WeightedParticleRepresentation, i) end

"""
    list_of_weights(ensemble::WeightedParticleRepresentation)

Return a list of particle weights from `ensemble`.
"""
list_of_weights(ens::WeightedParticleRepresentation) = [get_weight(ens, i) for i in 1:no_of_particles(ens)]

function eff_no_of_particles(ensemble::WeightedParticleRepresentation)
    sum = 0.
    @inbounds for i in 1:no_of_particles(ensemble)
        sum += get_weight(ensemble, i)^2
    end
    return 1/sum
end

dim(ensemle::WeightedParticleRepresentation) = (particle_dim(ensemble) + 1) * no_of_particles(ensemble)


text:
function plot_carrier_replica(data=DataFrame!(CSV.File("data/carrier_replica.csv")), targetdirtex="plots/carrier_replica.tex", targetdirpng="plots/carrier_replica.png")
    println("Plotting the carrier replica benhcmarks...")
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Samples",
            ylabel = "Zeit / us",
            ymode = "log",
            title = "Zeit für carrierreplica mit einer Antenne",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "outer north east"
            },
            PlotInc({
                blue,
                mark="|",
                thin
            }, Coordinates(data.Samples, 10^(-9)*data.sCPU_time)),
            PlotInc({
                red,
                mark="x",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.GPU_time)),
            PlotInc({
                green,
                mark="x",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.sGPU_time)),
            Legend(["CPU","GPU","sGPU"])
        )
    )
    pgfsave(targetdirtex, pgfplot)
    pgfsave(targetdirpng, pgfplot, dpi = 300)
    println("Saved the carrierreplica plot under ", targetdirtex)
end

function plot_carrier_replica_all(
    datamin=DataFrame!(CSV.File("data/carrierreplica_min.csv")),
    datamed=DataFrame!(CSV.File("data/carrierreplica_med.csv")),
    datamean=DataFrame!(CSV.File("data/carrierreplica_mean.csv")),
    targetdirtex="plots/carrierreplica_all.tex",
    targetdirpng="plots/carrierreplica_all.png"
)
    println("Plotting the carrier replica benhcmarks...")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Abtwastwerte",
            ylabel = "Zeit / s",
            ymode = "log",
            title = "Laufzeit der gencarrierreplica!",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "outer north east"
            },
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmin",
                style ="{solid}",
            }, Coordinates(datamin.Samples, 10^(-9)*datamin.sCPU_time)),
            PlotInc({
                blue,
                mark="x",
                thin,
                style ="{dashed}",
            }, Coordinates(datamed.Samples, 10^(-9)*datamed.sCPU_time)),
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmean",
                style ="{dotted}",
            }, Coordinates(datamean.Samples, 10^(-9)*datamean.sCPU_time)),
            PlotInc({
                thick,
                color = "blue",
                fill = "blue",
                opacity = 0.2 },
            raw"fill between [of=CPUmin and CPUmean]"),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmin",
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.GPU_time)),
            PlotInc({
                red,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(datamed.Samples, 10^(-9)*datamed.GPU_time)),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmean",
                style ="{dotted}",
            },Coordinates(datamean.Samples, 10^(-9)*datamean.GPU_time)),
            PlotInc({
                thick,
                color = "red",
                fill = "red",
                opacity = 0.2 },
            raw"fill between [of=GPUmin and GPUmean]"),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmin",
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.sGPU_time)),
            PlotInc({
                green,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(datamed.Samples, 10^(-9)*datamed.sGPU_time)),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmean",
                style ="{dotted}",
            },Coordinates(datamean.Samples, 10^(-9)*datamean.sGPU_time)),
            PlotInc({
                thick,
                color = "green",
                fill = "green",
                opacity = 0.2 },
            raw"fill between [of=sGPUmin and sGPUmean]"),
            Legend([
                "StructArray CPU Minimum", "StructArray CPU Median", "StructArray CPU Mean", "",
                "CuArray GPU Minimum", "CuArray GPU Median", "CuArray GPU Mean", "",
                "StructArray GPU Minimum", "StructArray GPU Median", "StructArray GPU Mean", ""])
        )
    )
    pgfsave(targetdirtex, pgfplot)
    pgfsave(targetdirpng, pgfplot, dpi = 300)
    println("Saved the carrierreplica plot")
end

function plot_code_replica_all(
    datamin=DataFrame!(CSV.File("data/codereplica_min_jetson.csv")),
    datamed=DataFrame!(CSV.File("data/codereplica_med_jetson.csv")),
    datamean=DataFrame!(CSV.File("data/codeprelica_mean_jetson.csv")),
    targetdirtex="plots/codereplica_jetson.tex",
    targetdirpng="plots/codereplica_jetson.png"
)
    println("Plotting the carrier replica benhcmarks...")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Abtwastwerte",
            ylabel = "Programmlaufzeit [s]",
            ymode = "log",
            title = "Jetson",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "outer north east"
            },
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmin",
                style ="{solid}",
            }, Coordinates(datamin.Samples, 10^(-9)*datamin.CPU_time)),
            PlotInc({
                blue,
                mark="x",
                thin,
                style ="{dashed}",
            }, Coordinates(datamed.Samples, 10^(-9)*datamed.CPU_time)),
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmean",
                style ="{dotted}",
            }, Coordinates(datamean.Samples, 10^(-9)*datamean.CPU_time)),
            PlotInc({
                thick,
                color = "blue",
                fill = "blue",
                opacity = 0.2 },
            raw"fill between [of=CPUmin and CPUmean]"),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmin",
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.GPU_time)),
            PlotInc({
                red,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(datamed.Samples, 10^(-9)*datamed.GPU_time)),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmean",
                style ="{dotted}",
            },Coordinates(datamean.Samples, 10^(-9)*datamean.GPU_time)),
            PlotInc({
                thick,
                color = "red",
                fill = "red",
                opacity = 0.2 },
            raw"fill between [of=GPUmin and GPUmean]"),
            Legend([
                "StructArray CPU Minimum", "StructArray CPU Median", "StructArray CPU Mean", "",
                "CuArray GPU Minimum", "CuArray GPU Median", "CuArray GPU Mean", "",])
        )
    )
    pgfsave(targetdirtex, pgfplot)
    pgfsave(targetdirpng, pgfplot, dpi = 300)
    println("Saved the codereplica plot")
end

function plot_downconvert(data=DataFrame!(CSV.File("data/downconvert.csv")), targetdirtex="plots/downconvert.tex", targetdirpng="plots/downconvert.png")
    println("Plotting the downconvert benhcmarks...")
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Samples",
            ylabel = "Zeit / s",
            ymode = "log",
            title = "Zeit für downconvert",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "north west",
            },
            PlotInc({
                blue,
                mark="|",
                thin
            }, Coordinates(data.Samples, 10^(-9)*data.sCPU_median_1ant)),
            PlotInc({
                blue,
                mark="|",
                thin
            }, Coordinates(data.Samples, 10^(-9)*data.sCPU_median_4ant)),
            PlotInc({
                blue,
                mark="|",
                thin
            }, Coordinates(data.Samples, 10^(-9)*data.sCPU_median_16ant)),
            PlotInc({
                red,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.GPU_median_1ant)),
            PlotInc({
                red,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.GPU_median_4ant)),
            PlotInc({
                red,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.GPU_median_16ant)),
            # PlotInc({
            #     gray,
            #     mark="|",
            #     thin
            # },Coordinates(data.Samples, 10^(-9)*data.sGPU_median_1ant)),
            # PlotInc({
            #     gray,
            #     mark="|",
            #     thin
            # },Coordinates(data.Samples, 10^(-9)*data.sGPU_median_4ant)),
            # PlotInc({
            #     gray,
            #     mark="|",
            #     thin
            # },Coordinates(data.Samples, 10^(-9)*data.sGPU_median_16ant)),
            Legend([
            "CPU median Ant=1","CPU median Ant=4","CPU median Ant=16",
            "GPU median Ant=1","GPU median Ant=4","GPU median Ant=16",
            #"sGPU median Ant=1","sGPU median Ant=4","sGPU median Ant=16"
            ])
        )
    )
    pgfsave(targetdirtex, pgfplot)
    pgfsave(targetdirpng, pgfplot, dpi = 300)
    println("Saved the downconvert plot under", targetdirtex)
end

function plot_downconvert_gpu_gpu(data=DataFrame!(CSV.File("data/downconvert.csv")), targetdirtex="plots/downconvert.tex", targetdirpng="plots/downconvert.png")
    println("Plotting the downconvert benhcmarks...")
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Samples",
            ylabel = "Zeit / s",
            ymode = "log",
            title = "Zeit für downconvert",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "outer north east",
            },
            PlotInc({
                red,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.GPU_median_1ant)),
            PlotInc({
                red,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.GPU_median_4ant)),
            PlotInc({
                red,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.GPU_median_16ant)),
            PlotInc({
                gray,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.sGPU_median_1ant)),
            PlotInc({
                gray,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.sGPU_median_4ant)),
            PlotInc({
                gray,
                mark="|",
                thin
            },Coordinates(data.Samples, 10^(-9)*data.sGPU_median_16ant)),
            Legend([
            # "CPU median Ant=1","CPU median Ant=4","CPU median Ant=16",
            "GPU median Ant=1","GPU median Ant=4","GPU median Ant=16",
            "sGPU median Ant=1","sGPU median Ant=4","sGPU median Ant=16"
            ])
        )
    )
    pgfsave(targetdirtex, pgfplot)
    pgfsave(targetdirpng, pgfplot, dpi = 300)
    println("Saved the downconvert plot under", targetdirtex)
end

function plot_downconvert_all(
    datamin=DataFrame!(CSV.File("data/donwconvert_min.csv")),
)
    println("Plotting the downconvert benhcmarks...")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Abtwastwerte",
            ylabel = "Programmlaufzeit [s]",
            ymode = "log",
            title = "Desktop PC",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "outer north east"
            },
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmin1ant",
                style ="{solid}",
            }, Coordinates(data.Samples, 10^(-9)*datamin.sCPU_time_1ant)),
            PlotInc({
                blue,
                mark="x",
                thin,
                style ="{dashed}",
            }, Coordinates(data.Samples, 10^(-9)*datamed.sCPU_time_1ant)),
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmean1ant",
                style ="{dotted}",
            }, Coordinates(data.Samples, 10^(-9)*datamean.sCPU_time_1ant)),
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmin4ant",
                style ="{solid}",
            }, Coordinates(data.Samples, 10^(-9)*datamin.sCPU_time_4ant)),
            PlotInc({
                thick,
                color = "blue",
                fill = "blue",
                opacity = 0.2 },
            raw"fill between [of=CPUmin4ant and CPUmean4ant]"),
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmin16ant",
                style ="{solid}",
            }, Coordinates(data.Samples, 10^(-9)*datamin.sCPU_time_16ant)),
            PlotInc({
                blue,
                mark="x",
                thin,
                style ="{dashed}",
            }, Coordinates(data.Samples, 10^(-9)*datamed.sCPU_time_16ant)),
            PlotInc({
                blue,
                mark="x",
                thin,
                "name path=CPUmean16ant",
                style ="{dotted}",
            }, Coordinates(data.Samples, 10^(-9)*datamean.sCPU_time_16ant)),
            PlotInc({
                thick,
                color = "blue",
                fill = "blue",
                opacity = 0.2 },
            raw"fill between [of=CPUmin16ant and CPUmean16ant]"),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmin1ant",
                style ="{solid}",
            },Coordinates(data.Samples, 10^(-9)*datamin.GPU_time_1ant)),
            PlotInc({
                red,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(data.Samples, 10^(-9)*datamed.GPU_time_1ant)),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmean1ant",
                style ="{dotted}",
            },Coordinates(data.Samples, 10^(-9)*datamean.GPU_time_1ant)),
            PlotInc({
                thick,
                color = "red",
                fill = "red",
                opacity = 0.2 },
            raw"fill between [of=GPUmin1ant and GPUmean1ant]"),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmin4ant",
                style ="{solid}",
            },Coordinates(data.Samples, 10^(-9)*datamin.GPU_time_4ant)),
            PlotInc({
                red,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(data.Samples, 10^(-9)*datamed.GPU_time_4ant)),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmean4ant",
                style ="{dotted}",
            },Coordinates(data.Samples, 10^(-9)*datamean.GPU_time_4ant)),
            PlotInc({
                thick,
                color = "red",
                fill = "red",
                opacity = 0.2 },
            raw"fill between [of=GPUmin4ant and GPUmean4ant]"),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmin16ant",
                style ="{solid}",
            },Coordinates(data.Samples, 10^(-9)*datamin.GPU_time_16ant)),
            PlotInc({
                red,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(data.Samples, 10^(-9)*datamed.GPU_time_16ant)),
            PlotInc({
                red,
                mark="x",
                thin,
                "name path=GPUmean16ant",
                style ="{dotted}",
            },Coordinates(data.Samples, 10^(-9)*datamean.GPU_time_16ant)),
            PlotInc({
                thick,
                color = "red",
                fill = "red",
                opacity = 0.2 },
            raw"fill between [of=GPUmin16ant and GPUmean16ant]"),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmin1ant",
                style ="{solid}",
            },Coordinates(data.Samples, 10^(-9)*datamin.sGPU_time_1ant)),
            PlotInc({
                green,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(data.Samples, 10^(-9)*datamed.sGPU_time_1ant)),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmean1ant",
                style ="{dotted}",
            },Coordinates(data.Samples, 10^(-9)*datamean.sGPU_time_1ant)),
            PlotInc({
                thick,
                color = "green",
                fill = "green",
                opacity = 0.2 },
            raw"fill between [of=sGPUmin1ant and sGPUmean1ant]"),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmin4ant",
                style ="{solid}",
            },Coordinates(data.Samples, 10^(-9)*datamin.sGPU_time_4ant)),
            PlotInc({
                green,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(data.Samples, 10^(-9)*datamed.sGPU_time_4ant)),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmean4ant",
                style ="{dotted}",
            },Coordinates(data.Samples, 10^(-9)*datamean.sGPU_time_4ant)),
            PlotInc({
                thick,
                color = "green",
                fill = "green",
                opacity = 0.2 },
            raw"fill between [of=sGPUmin4ant and sGPUmean4ant]"),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmin16ant",
                style ="{solid}",
            },Coordinates(data.Samples, 10^(-9)*datamin.sGPU_time_16ant)),
            PlotInc({
                green,
                mark="x",
                thin,
                style ="{dashed}",
            },Coordinates(data.Samples, 10^(-9)*datamed.sGPU_time_16ant)),
            PlotInc({
                green,
                mark="x",
                thin,
                "name path=sGPUmean16ant",
                style ="{dotted}",
            },Coordinates(data.Samples, 10^(-9)*datamean.sGPU_time_16ant)),
            PlotInc({
                thick,
                color = "green",
                fill = "green",
                opacity = 0.2 },
            raw"fill between [of=sGPUmin16ant and sGPUmean16ant]"),
            Legend([
                "StructArray CPU Minimum Ant=1", "StructArray CPU Median Ant=1", "StructArray CPU Mean Ant=1", "",
                "StructArray CPU Minimum Ant=4", "StructArray CPU Median Ant=4", "StructArray CPU Mean Ant=4", "",
                "StructArray CPU Minimum Ant=16", "StructArray CPU Median Ant=16", "StructArray CPU Mean Ant=16", "",
                "CuArray GPU Minimum Ant=1", "CuArray GPU Median Ant=1", "CuArray GPU Mean Ant=1", "",
                "CuArray GPU Minimum Ant=4", "CuArray GPU Median Ant=4", "CuArray GPU Mean Ant=4", "",
                "CuArray GPU Minimum Ant=16", "CuArray GPU Median Ant=16", "CuArray GPU Mean Ant=16", "",
                "StructArray GPU Minimum Ant=1", "StructArray GPU Median Ant=1", "StructArray GPU Mean Ant=1", "",
                "StructArray GPU Minimum Ant=4", "StructArray GPU Median Ant=4", "StructArray GPU Mean Ant=4", "",
                "StructArray GPU Minimum Ant=16", "StructArray GPU Median Ant=16", "StructArray GPU Mean Ant=16", ""])
        )
    )
    pgfsave(targetdirtex, pgfplot)
    pgfsave(targetdirpng, pgfplot, dpi = 300)
    println("Saved the downconvert plot")
end

function plot_correlate()
    println("Plotting the correlation benhcmarks...")
    data = DataFrame!(CSV.File("data/correlate.csv"))
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Samples",
            ylabel = "Zeit / us",
            title = "Zeit für correlate mit einer Antenne",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "north west"
            },
            PlotInc({
                blue,
                mark="|",
                thin
            },
            Coordinates(data.Samples, data.sCPU_median)),
            PlotInc({
                red,
                mark="x",
                thin
            },
            Coordinates(data.Samples, data.GPU_median)),
            Legend(["CPU median","GPU median"])
        )
    )
    pgfsave("plots/correlate.tex", pgfplot)
    pgfsave("plots/correlate.png", pgfplot, dpi = 300)
    println("Saved the correlations plot under /plots")
end

function plot_correlate_all(
    datamin=DataFrame!(CSV.File("data/correlate_min.csv"))
)
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Abtwastwerte",
            ylabel = "Programmlaufzeit [s]",
            ymode = "log",
            title = "correlate",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "outer north east"
            },
            PlotInc({
                cyan,
                mark="x",
                thin,
                style ="{solid}",
            }, Coordinates(datamin.Samples, 10^(-9)*datamin.sCPU_time_1ant)),
            PlotInc({
                teal,
                mark="x",
                thin,
                style ="{solid}",
            }, Coordinates(datamin.Samples, 10^(-9)*datamin.sCPU_time_4ant)),
            PlotInc({
                blue,
                mark="x",
                thin,
                style ="{solid}",
            }, Coordinates(datamin.Samples, 10^(-9)*datamin.sCPU_time_16ant)),
            PlotInc({
                color="magenta",
                mark="x",
                thin,
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.GPU_time_1ant)),
            PlotInc({
                color="red",
                mark="x",
                thin,
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.GPU_time_4ant)),
            PlotInc({
                color="violet",
                mark="x",
                thin,
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.GPU_time_16ant)),
            PlotInc({
                color="lime",
                mark="x",
                thin,
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.sGPU_time_1ant)),
            PlotInc({
                color="green",
                mark="x",
                thin,
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.sGPU_time_4ant)),
            PlotInc({
                color="green!50!black",
                mark="x",
                thin,
                style ="{solid}",
            },Coordinates(datamin.Samples, 10^(-9)*datamin.sGPU_time_16ant)),
            Legend([
                "StructArray CPU Minimum Ant=1",
                "StructArray CPU Minimum Ant=4",
                "StructArray CPU Minimum Ant=16",
                "CuArray GPU Minimum Ant=1",
                "CuArray GPU Minimum Ant=4",
                "CuArray GPU Minimum Ant=16",
                "StructArray GPU Minimum Ant=1",
                "StructArray GPU Minimum Ant=4",
                "StructArray GPU Minimum Ant=16"])
        )
    )
    pgfsave("plots/correlate.tex", pgfplot)
    pgfsave("plots/correlate.png", pgfplot, dpi = 300)
    println("Saved the correlations plot under /plots")

end

function plot_code_replica()
    println("Plotting the codereplica benhcmarks...")
    data = DataFrame!(CSV.File("data/codereplica.csv"))
    pgfplot = @pgf TikzPicture(
        Axis({
            xlabel = "Samples",
            ylabel = "Zeit / us",
            title = "Zeit für codereplica mit einer Antenne",
            xmajorgrids,
            ymajorgrids,
            scaled_ticks = "false",
            legend_pos = "north west"
            },
            PlotInc({
                blue,
                mark="|",
                thin
            },
            Coordinates(data.Samples, data.CPU_median)),
            PlotInc({
                red,
                mark="x",
                thin
            },
            Coordinates(data.Samples, data.GPU_median)),
            Legend(["CPU median","GPU median"])
        )
    )
    pgfsave("plots/codereplica.tex", pgfplot)
    pgfsave("plots/codereplica.png", pgfplot, dpi = 300)
    println("Saved the codereplica plot under /plots")
end


text:
module ReleaseCmd

using LibGit2
using OrderedCollections
using UUIDs
using GitHub
using Crayons.Box
using TOML
using Pkg
using Pkg.Types
using Pkg.Types: RegistrySpec
using Comonicon.Tools: prompt, cmd_error
using ..Options
using ..PkgCmd
using ..SearchCmd: find_max_version, collect_registries
using ..IonBase: IonBase, read_github_auth, gitcmd

"""
    PRN{name}

Package Registry Name
"""
struct PRN{name} end

"""
    PRN(name::String)

Create a `PRN` (Pacakge Registry Name) object.
"""
PRN(name::String) = PRN{Symbol(name)}()

macro PRN_str(name::String)
    return PRN{Symbol(name)}
end

Base.show(io::IO, ::PRN{registry}) where {registry} = print(io, "Pacakge Registry ", string(registry))

Base.@kwdef struct VersionTokens
    major::String = "major"
    minor::String = "minor"
    patch::String = "patch"
end

const VERSION_TOKENS = VersionTokens()
Base.show(io::IO, vt::VersionTokens) = print(io, "(", vt.major, ", ", vt.minor, ", ", vt.patch, ")")

Base.in(version::String, tokens::VersionTokens) = (version == tokens.major) ||
    (version == tokens.minor) || (version == tokens.patch)

function is_version_number(version)
    occursin(r"[0-9]+.[0-9]+.[0-9]+", version) ||
        occursin(r"v[0-9]+.[0-9]+.[0-9]+", version)
end

struct Project
    path::String
    toml::String
    pkg::Pkg.Types.Project
    ion::Options.Ion
    git::Cmd
    branch::String
    quiet::Bool
end

function Project(path::String=pwd(); branch="master", quiet=false, ion::Options.Ion=Options.read())
    toml = Base.current_project(path)
    toml === nothing && cmd_error("cannot find (Julia)Project.toml in $path")
    path = dirname(toml)
    pkg = Pkg.Types.read_project(toml)
    git = gitcmd(path)
    return Project(path, toml, pkg, ion, git, branch, quiet)
end

Base.show(io::IO, p::Project) = print(io, "Project(", p.path, ")")

function current_branch(p::Project)
    return readchomp(`$(p.git) rev-parse --abbrev-ref HEAD`)
end

function commit_toml(project::Project; push::Bool=false)
    git = project.git
    version_number = project.pkg.version
    run(`$git add $(project.toml)`)
    run(`$git commit -m"bump version to $version_number"`)
    push && gitpush(project)
    return
end

function gitpush(project::Project)
    run(`$(project.git) push origin $(project.branch)`)
    println(DARK_GRAY_FG("="^80))
end

function reset_last_commit(project::Project; push=false)
    git = project.git
    run(`$git revert --no-edit --no-commit HEAD`)
    run(`$git commit -m "revert version bump due to an error occured in IonCLI"`)
    push && gitpush(project)
    return
end

function checkout(f, p::Project)
    old = current_branch(p)

    if old != p.branch
        @info "checking out to $(p.branch)"
        run(`$(p.git) checkout $(p.branch)`)
    end

    f()

    if old != p.branch
        run(`$(p.git) checkout $old`)
    end
end

function query_project_registry(project::Project)
    registries = collect_registries(project.ion)

    @info "filtering registries"
    matches = filter(registries) do rs::RegistrySpec
        d = Pkg.Types.read_registry(joinpath(rs.path, "Registry.toml"))
        haskey(d["packages"], string(project.pkg.uuid))
    end

    isempty(matches) && return
    return matches
end

function register(registry::String, project::Project)
    if isempty(registry) # registered package
        @info "querying project registry"
        matches = query_project_registry(project)
        matches === nothing && cmd_error(
            "$(project.pkg.name) is not registered " *
            "in local registries, please specify " *
            "registry name using --registry=<name>"
        )

        length(matches) == 1 || cmd_error(
            "this package is registered in the following registries: " *
            join([isnothing(each.name) ? "unknown" : each.name for each in matches], ", ") *
            "please specify registry name using --registry=<name>"
        )

        path = first(matches).path
        return register(PRN(basename(path)), project)
    else
        return register(PRN(registry), project)
    end
end

function register(registry::PRN, project::Project)
    cmd_error("register workflow is not defined for $registry")
end

function registrator_msg(project)
    msg = "Released via [Ion CLI](https://github.com/Roger-luo/IonCLI.jl)\n"
    msg *= "@JuliaRegistrator register"
    # msg *= "testing"
    if project.branch == "master"
        return msg
    else
        return msg * " branch=$(project.branch)"
    end
end

function read_head(git, branch="master")
    return readchomp(`$git rev-parse --verify HEAD`)
end

function read_remote_push(git, remote="origin")
    return readchomp(`$git config --get remote.$remote.url`)
end

function github_repo(git, remote="origin")
    url = read_remote_push(git, remote)
    github_https = "https://github.com/"
    github_ssh = "git@github.com:"
    if startswith(url, github_https)
        if endswith(url, ".git")
            return url[length(github_https)+1:end-4]
        else
            return url[length(github_https)+1:end]
        end
    elseif startswith(url, github_ssh)
        return url[length(github_ssh)+1:end-4]
    else
        return
    end
end

function update_version!(project::Project, version)
    if is_version_number(version)
        version_number = VersionNumber(version)
    elseif version in VERSION_TOKENS
        version_number = bump_version(project, version)
    else
        cmd_error("invalid version $version")
    end

    if !project.quiet
        latest_version = find_max_version(project.pkg.name)

        if latest_version === nothing
            println("package not found in local registries")
        else
            println("latest registered version: ", LIGHT_CYAN_FG(string(latest_version)))
            if latest_version > version_number
                @warn "input version is smaller than registered version"
            end
        end

        println(" "^10, "current version: ", LIGHT_CYAN_FG(string(project.pkg.version)))
        println(" "^7, "version to release: ", LIGHT_CYAN_FG(string(version_number)))
        if !prompt("do you want to update Project.toml?")
            exit(0)
        end
    end

    write_version(project, version_number)
    println(" ", LIGHT_GREEN_FG("✔"), "  Project.toml has been updated to ", LIGHT_CYAN_FG(string(version_number)))
    return project
end

function bump_version(project::Project, token::String)
    if project.pkg.version === nothing
        return bump_version(v"0.0.0", token)
    else
        return bump_version(project.pkg.version, token)
    end
end

function bump_version(version::VersionNumber, token::String)
    if token == VERSION_TOKENS.major
        return VersionNumber(version.major+1, 0, 0)
    elseif token == VERSION_TOKENS.minor
        return VersionNumber(version.major, version.minor+1, 0)
    elseif token == VERSION_TOKENS.patch
        return VersionNumber(version.major, version.minor, version.patch+1)
    else
        cmd_error("invalid token $token")
    end
end

function write_version(project::Project, version::VersionNumber)
    project.pkg.version = version
    open(project.toml, "w+") do f # following whatever Pkg does
        TOML.print(f, to_dict(project); sorted=true, by=key -> (Pkg.Types.project_key_order(key), key)) do x
            if x isa UUID || x isa VersionNumber
                x = string(x)
            end
            x
        end
    end
end

to_dict(p::Project) = to_dict(p.pkg)

function to_dict(project::Pkg.Types.Project)
    project_keys = [:name, :uuid, :authors, :version, :deps, :compat, :extras, :targets]
    t = []
    for key in project_keys
        push_maybe!(t, project, key)
    end

    # write other part back
    for key in keys(project.other)
        if !(key in string.(project_keys))
            push!(t, project.other[key])
        end
    end
    return OrderedDict(t)
end

function push_maybe!(t::Vector, project::Pkg.Types.Project, name::Symbol)
    key = string(name)
    if hasfield(Pkg.Types.Project, name)
        member = getfield(project, name)
        if member !== nothing
            push!(t, key => member)
        end
    else
        member = get(project.other, key, nothing)
        if member !== nothing
            push!(t, key => member)
        end
    end
    return t
end

function release(version::String, path::String=pwd(), registry="", branch="master")
    project = Project(path; branch=branch)
    # new version needs to be pushed
    # so the JuliaRegistrator can find
    # it later
    checkout(project) do
        if LibGit2.isdirty(LibGit2.GitRepo(project.path))
            cmd_error("package repository is dirty, please commit or stash changes.")
        end

        committed_changes = false
        if version != "current"
            IonBase.Registry.update()
            update_version!(project, version)
            commit_toml(project; push=true)
            committed_changes = true
        elseif prompt("do you want to release current version $(LIGHT_CYAN_FG(string(project.pkg.version)))?")
            # make sure we sync with remote
            gitpush(project)
        else # do nothing
            return
        end

        try
            register(registry, project)
        catch e
            if committed_changes
                reset_last_commit(project; push=true)
                println(" ", LIGHT_GREEN_FG("✔"), "  revert version bump commit:")
            end
        end
    end
    return
end

function useless_animation(auth::Base.Event, summon::Base.Event, interrupted_or_done::Base.Event)
    anim_chars = ["◐","◓","◑","◒"]
    ansi_enablecursor = "\e[?25h"
    ansi_disablecursor = "\e[?25l"
    t = Timer(0; interval=1/20)
    print_lock = ReentrantLock()
    printloop_should_exit = interrupted_or_done.set
    return @async begin
        try
            count = 1
            while !printloop_should_exit
                lock(print_lock) do
                    print(ansi_disablecursor)
                    print("\e[1G  ", CYAN_FG(anim_chars[mod1(count, 4)]))
                    print("    ")
                    if !auth.set
                        print("authenticating...")
                    elseif !summon.set
                        print("summoning JuliaRegistrator...")
                    else
                        printloop_should_exit = true
                    end
                end
                printloop_should_exit = interrupted_or_done.set
                count += 1
                wait(t)
            end
        catch e
            notify(interrupted_or_done)
            lock(print_lock) do
                println("\e[1G  ", RED_FG("❌"), "  fail to register $(project.pkg.name), error msg:")
                println(e.msg)
            end
            rethrow(e)
        finally
            print(ansi_enablecursor)
        end
    end
end

function register(::PRN"General", project::Project)
    @info "register project in General"
    auth_done = Base.Event()
    summon_done = Base.Event()
    interrupted_or_done = Base.Event()

    print_task = useless_animation(auth_done, summon_done, interrupted_or_done)

    github_token = read_github_auth()
    auth = GitHub.authenticate(github_token)
    notify(auth_done)

    HEAD = read_head(project.git)
    comment_json = Dict{String, Any}(
        "body" => registrator_msg(project),
    )

    repo = github_repo(project.git)
    if repo === nothing
        cmd_error("not a GitHub repository")
    end

    comment = GitHub.create_comment(repo, HEAD, :commit; params=comment_json, auth=auth)
    notify(summon_done)
    notify(interrupted_or_done)
    wait(print_task)
    println("\e[1G ", LIGHT_GREEN_FG("✔"), "  JuliaRegistrator has been summoned, check it in the following URL:")
    println("  ", CYAN_FG(string(comment.html_url)))
    return comment
end

end

"""
release a package.

# Arguments

- `version`: version number you want to release. Can be a specific version, "current" or either of $(ReleaseCmd.VERSION_TOKENS)
- `path`: path to the project you want to release.

# Options

- `-r,--registry <registry name>`: registry you want to register the package.
    If the package has not been registered, ion will try to register
    the package in the General registry. Or the user needs to specify
    the registry to register using this option.

- `-b, --branch <branch name>`: branch you want to register.
"""
@cast function release(version::String="current", path::String=pwd(); registry="", branch="master")
    ReleaseCmd.release(version, path, registry, branch)
    return
end


text:
type MarketModelCashFlow
  timeIndex::Int
  amount::Float64
end

MarketModelCashFlow() = MarketModelCashFlow(1, 0.0)
clone(mmcf::MarketModelCashFlow) = MarketModelCashFlow(mmcf.timeIndex, mmcf.amount)

type MarketModelPathWiseCashFlow
  timeIndex::Int
  amount::Vector{Float64}
end

MarketModelPathWiseCashFlow(n::Int) = MarketModelPathWiseCashFlow(1, zeros(n))

## Discounter ##
type MarketModelDiscounter
  beforeSize::Int
  beforeWeight::Float64
end

function MarketModelDiscounter(paymentTime::Float64, rateTimes::Vector{Float64})
  check_increasing_times(rateTimes)
  beforeTime = findfirst(rateTimes, paymentTime)

  # handle the case where the payment is in the last period or after the last period
  if beforeTime > length(rateTimes) - 1
    beforeTime = length(rateTimes) - 1
  end

  if beforeTime == 0
    if paymentTime > rateTimes[end]
      beforeTime = length(rateTimes) - 1
    else
      beforeTime = 1
    end
  end

  beforeWeight = 1.0 - (paymentTime - rateTimes[beforeTime]) / (rateTimes[beforeTime + 1] - rateTimes[beforeTime])

  return MarketModelDiscounter(beforeTime, beforeWeight)
end

function numeraire_bonds(mmd::MarketModelDiscounter, curveState::CurveState, numeraire::Int)
  preDF = discount_ratio(curveState, mmd.beforeSize, numeraire)
  if mmd.beforeWeight == 1.0
    return preDF
  end

  postDF = discount_ratio(curveState, mmd.beforeSize+1, numeraire)
  if mmd.beforeWeight == 0.0
    return postDF
  end

  return ^(preDF, mmd.beforeWeight) * ^(postDF, 1.0 - mmd.beforeWeight)
end

type MarketModelPathwiseDiscounter
  beforeTimeIdx::Int
  numberRates::Int
  beforeWeight::Float64
  postWeight::Float64
  taus::Vector{Float64}
end

function MarketModelPathwiseDiscounter(paymentTime::Float64, rateTimes::Vector{Float64})
  check_increasing_times(rateTimes)

  numberRates = length(rateTimes) - 1
  beforeTimeIdx = findfirst(rateTimes, paymentTime)

  # handle the case where the payment is in the last period or after the last period
  if beforeTimeIdx > length(rateTimes) - 1
    beforeTimeIdx = length(rateTimes) - 1
  end

  beforeWeight = 1.0 - (paymentTime - rateTimes[beforeTimeIdx]) / (rateTimes[beforeTimeIdx + 1] - rateTimes[beforeTimeIdx])
  postWeight = 1.0 - beforeWeight
  taus = Vector{Float64}(numberRates)
  for i in eachindex(taus)
    taus[i] = rateTimes[i+1] - rateTimes[i]
  end

  return MarketModelPathwiseDiscounter(beforeTimeIdx, numberRates, beforeWeight, postWeight, taus)
end

function get_factors(mmdisc::MarketModelPathwiseDiscounter, ::Matrix{Float64}, Discounts::Matrix{Float64}, currentStep::Int)
  factors = Vector{Float64}(mmdisc.numberRates + 1)
  preDF = Discounts[currentStep, mmdisc.beforeTimeIdx]
  postDF = Discounts[currentStep, mmdisc.beforeTimeIdx + 1]

  for i = mmdisc.beforeTimeIdx+1:mmdisc.numberRates
    factors[i+1] = 0.0
  end

  if mmdisc.postWeight == 0.0
    factors[1] = preDF

    for i = 1:mmdisc.beforeTimeIdx
      factors[i+1] = -preDF * mmdisc.taus[i] * Discounts[currentStep, i+1] / Discounts[currentStep, i]
    end

    factors[mmdisc.beforeTimeIdx + 1] = 0.0

    return factors
  end

  df = preDF * ^(postDF / preDF, mmdisc.postWeight)

  factors[1] = df

  for i = 1:mmdisc.beforeTimeIdx
    factors[i+1] = -df * mmdisc.taus[i] * Discounts[currentStep, i+1] / Discounts[currentStep, i]
  end

  factors[mmdisc.beforeTimeIdx + 1] *= mmdisc.postWeight

  return factors
end

function check_increasing_times(times::Vector{Float64})
  length(times) > 0 || error("at least one time is required")
  times[1] > 0.0 || error("first time must be greater than zero")

  issorted(times) || error("non increasing rate times")
end

function check_increasing_times_and_calculate_taus(times::Vector{Float64})
  check_increasing_times(times)

  nTimes = length(times)
  nTimes > 1 || error("at least two times are required")

  taus = Vector{Float64}(nTimes - 1)

  for i in eachindex(taus)
    taus[i] = times[i+1] - times[i]
  end

  return taus
end

function merge_times(times::Vector{Vector{Float64}})
  allTimes = Vector{Float64}()
  for i in eachindex(times)
    allTimes = vcat(allTimes, times[i])
  end

  sort!(allTimes)
  allTimes = unique(allTimes)

  isPresent = similar(times, BitArray{1})

  for i in eachindex(times)
    isPresent[i] = BitArray{1}(length(times[i]))
    for j in eachindex(allTimes)
      isPresent[i][j] = findfirst(times[i], allTimes[j]) > 0 ? true : false
    end
  end

  return allTimes, isPresent
end

function is_in_subset(mainSet::Vector{Float64}, subSet::Vector{Float64})
  result = falses(length(mainSet))
  dimsubSet = length(subSet)
  if dimsubSet == 0
    return result
  end
  dimSet = length(mainSet)

  setElement = subSetElement = 0.0

  dimSet >= dimsubSet || error("set is required to be larger or equal than subset")

  for i in eachindex(mainSet) # loop in set
    j = 1
    setElement = mainSet[i]
    while true # loop in subset
      subSetElement = subSet[j]
      result[i] = false
      # if smaller no hope, leave false and go to next i
      if setElement < subSetElement
        break
      end
      # match!
      if setElement == subSetElement
        result[i] = true
        break
      end

      # if larger, leave false if at end or go to next j
      if j == dimsubSet
        break
      end

      j += 1
    end
  end

  return result
end


text:
abstract type Arr end
abstract type Op end
abstract type Fn end

struct Α <: Arr end
struct Ω <: Arr end
struct Apply{T} <: Arr # This could both be a funciton or an array....
    f::T
    r::Arr
end
struct Apply2{T} <: Arr
    f::T
    l::Arr
    r::Arr
end
struct JlVal <: Arr
    val
end

struct Op1Sym{c} <: Op end

struct Op1{c, L} <: Fn
    l::L
end
struct OpSymPair{c} <: Op
    r # Op1Sym
end

struct Op2Sym{c} <: Op end
struct Op2Partial{c, R} <: Op
    r::R
end
struct Op2{c, L, R} <: Fn
    l::L
    r::R
end

struct ConcArr <: Arr
    l
    r
end

struct PrimFn{c} <: Fn end

struct UDefFn{arity} <: Fn
    ast
end

# Library functions
digit(V,n)=(mod(n,V),(n-mod(n,V))/V)
function ⊤(V,state) # alan's encode function
   V, state
   z=copy(V)
   for i=length(V):-1:1
        z[i], state = digit(z[i],state)
   end
   z
end
function ⊤(V,n::AbstractArray)
    V′ = convert(Array, V)
    n′ = convert(Array, n)
    reduce(vcat, [⊤(V′,nn)' for nn in n′])'
end
ι(V,W)=[v∈W ? findfirst(isequal(v), W) : 1+length(W) for v in V]

function prefix_scan(f, x, v0) # optimizations are just dispatch on f!
    fst = f(v0, x[1])
    y = Array(typeof(fst), length(x))
    y[1] = fst
    for i=2:length(y)
      y[i] =f(y[i-1], x[i])
    end
    y
end

# More definitions
const prim_fns = [
  ( '+', (:.+, :.+)),
  ( '-', (:-, :-)),
  ( '÷', (:(1/ω), :/)),
  ( '×', (:sign, :(α.*ω))),
  ( '⌈', (:ceil, :max)),
  ( '⌊', (:floor, :min)),
  ( ',', (:(ω[:]), :vcat)),
  ( '*', (:exp, :(.^))),
  ( '!', (:factorial, :(factorial(ω) / (factorial(α) * factorial(ω-α))))), # Factorial on floats
  ( '|', (:abs, :(ω%α))),
  ( '⊛', (:log, :log)),
  ( '○', (:(π*ω), nothing)),
  ( '~', (:!, nothing)),
  ( '⌹', (:inv, :\)),
  ( 'ι', (:(1:ω[1]), :(ι(α, ω)))), # ω[1] is not really right
  ( '?', (:(rand(1:ω[1])), nothing)),
  ( '>', (nothing, :(.>))),
  ( '<', (nothing, :(.<))),
  ( '=', (nothing, :(.==))),
  ( '≠', (nothing, :(.!=))),
  ( '≤', (nothing, :≤)),
  ( '≥', (nothing, :≥)),
  ( '∧', (nothing, :&)),
  ( '∨', (nothing, :|)),
  ( ']', (nothing, :(getindex(ω, α)))),
  ( '⊤', (nothing, :⊤)),
  ( '⍋', (:sortperm, nothing)),
  ( '⍒', (:(sortperm(ω, rev=true)), nothing)),
  ( 'ρ', (:([size(ω)...]), :(reshape(ω, tuple(α...))))),
  ( '⍉', (:(x->permutedims(x, ntuple(i->ndims(x) - i + 1, ndims(x)))), :permutedims)),
  ( '⌽', (:(flipdim(ω, 1)), nothing)),
  ( '↓', (nothing, :(ω[α+1:length(ω)]))),
  ( '↑', (nothing, :(ω[1:α]))),
]
const function_names = map(x->x[1], prim_fns)
const numeric = ['⁻', ' ', '0':'9'..., '.']
const monadic_operators = ['\\', '/', '⌿', '⍀', '¨', '↔', '∘', '⍨']
const dyadic_operators = ['.', '⋅']
const argnames = ['α', 'ω']

# Identity elements
identity(::PrimFn{'+'}, T::Type) = zero(T)
identity(::PrimFn{'+'}, ::Type{Bool}) = 0
identity(::PrimFn{'-'}, ::Type{Bool}) = 0
identity(::PrimFn{'-'}, T::Type) = zero(T)
identity(::PrimFn{'×'}, T::Type) = one(T)
identity(::PrimFn{'∨'}, T::Type) = zero(T)
identity(::PrimFn{'∧'}, T::Type) = one(T)
identity(::PrimFn{'⌈'}, T::Type) = typemin(T)
identity(::PrimFn{'⌊'}, T::Type) = typemax(T)
identity(::PrimFn{'='}, T::Type) = true
identity(::PrimFn{','}, T::Type) = [] # ...
identity(op::Op1, T::Type) = identity(op.l, T)
identity(x,T::Type) = error("Could not determine the identity element for the operation $(x) for argtype $T")


text:
using Test
using IntervalSets

import ClimaCore
import ClimaCoreMakie: ClimaCoreMakie, Makie

if haskey(ENV, "CI")
    using GLMakie
end
OUTPUT_DIR = mkpath(get(ENV, "CI_OUTPUT_DIR", tempname()))

@testset "spectral element 2D cubed-sphere" begin
    R = 6.37122e6

    domain = ClimaCore.Domains.SphereDomain(R)
    mesh = ClimaCore.Meshes.EquiangularCubedSphere(domain, 6)
    grid_topology = ClimaCore.Topologies.Topology2D(mesh)
    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
    space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = ClimaCore.Fields.coordinate_field(space)

    u = map(coords) do coord
        u0 = 20.0
        α0 = 45.0
        ϕ = coord.lat
        λ = coord.long

        uu = u0 * (cosd(α0) * cosd(ϕ) + sind(α0) * cosd(λ) * sind(ϕ))
        uv = -u0 * sind(α0) * sind(λ)
        ClimaCore.Geometry.UVVector(uu, uv)
    end

    fig = ClimaCoreMakie.viz(u.components.data.:1)
    @test fig !== nothing

    fig_png = joinpath(OUTPUT_DIR, "2D_cubed_sphere.png")
    GLMakie.save(fig_png, fig)
    @test isfile(fig_png)
end

@testset "spectral element rectangle 2D" begin
    domain = ClimaCore.Domains.RectangleDomain(
        ClimaCore.Geometry.XPoint(0) .. ClimaCore.Geometry.XPoint(2π),
        ClimaCore.Geometry.YPoint(0) .. ClimaCore.Geometry.YPoint(2π),
        x1periodic = true,
        x2periodic = true,
    )

    n1, n2 = 2, 2
    Nq = 4
    mesh = ClimaCore.Meshes.RectilinearMesh(domain, n1, n2)
    grid_topology = ClimaCore.Topologies.Topology2D(mesh)
    #quad = ClimaCore.Spaces.Quadratures.GLL{Nq}()
    quad = ClimaCore.Spaces.Quadratures.ClosedUniform{Nq + 1}()
    space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = ClimaCore.Fields.coordinate_field(space)

    sinxy = map(coords) do coord
        cos(coord.x + coord.y)
    end

    fig = ClimaCoreMakie.viz(sinxy)
    @test fig !== nothing

    fig_png = joinpath(OUTPUT_DIR, "2D_rectangle.png")
    GLMakie.save(fig_png, fig)
    @test isfile(fig_png)
end


text:
function Cos()
    return Cos(())
end

function value(obj::Cos, arg0::Derivative)
    return jcall(obj, "value", Derivative, (Derivative,), arg0)
end

function value(obj::Cos, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::Cos, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq
using DomainSets

# using Plots

# local sol
begin #@testset "Test 01: Brusselator equation 2D" begin
       @parameters x y t
       @variables u(..) v(..)
       Dt = Differential(t)
       Dx = Differential(x)
       Dy = Differential(y)
       Dxx = Differential(x)^2
       Dyy = Differential(y)^2

       ∇²(u) = Dxx(u) + Dyy(u)

       brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.

       x_min = y_min = t_min = 0.0
       x_max = y_max = 1.0
       t_max = 11.5

       α = 10.

       u0(x,y,t) = 22(y*(1-y))^(3/2)
       v0(x,y,t) = 27(x*(1-x))^(3/2)

       eq = [Dt(u(x,y,t)) ~ 1. + v(x,y,t)*u(x,y,t)^2 - 4.4*u(x,y,t) + α*∇²(u(x,y,t)) + brusselator_f(x, y, t),
             Dt(v(x,y,t)) ~ 3.4*u(x,y,t) - v(x,y,t)*u(x,y,t)^2 + α*∇²(v(x,y,t))]

       domains = [x ∈ Interval(x_min, x_max),
                  y ∈ Interval(y_min, y_max),
                  t ∈ Interval(t_min, t_max)]

       bcs = [u(x,y,0) ~ u0(x,y,0),
              u(0,y,t) ~ u(1,y,t),
              u(x,0,t) ~ u(x,1,t),

              v(x,y,0) ~ v0(x,y,0),
              v(0,y,t) ~ v(1,y,t),
              v(x,0,t) ~ v(x,1,t)]

       @named pdesys = PDESystem(eq,bcs,domains,[x,y,t],[u(x,y,t),v(x,y,t)])

       # Method of lines discretization
       N = 32

       dx = 1/N
       dy = 1/N

       order = 2

       discretization = MOLFiniteDifference([x=>dx, y=>dy], t, approx_order=order)

       #MethodOfLines.generate_code(pdesys, discretization)
       # Convert the PDE problem into an ODE problem
       println("Discretization:")
       @time prob = discretize(pdesys,discretization)

       println("Solve:")
       @time sol = solve(prob, TRBDF2(),saveat=0.01)

       # Solve reference problem

       xyd_brusselator = range(0,stop=1,length=N)
       brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
       limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a
       function brusselator_2d_loop(du, u, p, t)
              A, B, alpha, dx = p
              alpha = alpha/dx^2
              @inbounds for I in CartesianIndices((N, N))
                     i, j = Tuple(I)
           x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
           ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
           du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
           B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
           du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
                       A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
           end
       end
       p = (3.4, 1., 10., step(xyd_brusselator))

       function init_brusselator_2d(xyd)
           N = length(xyd)
           u = zeros(N, N, 2)
           for I in CartesianIndices((N, N))
               x = xyd[I[1]]
               y = xyd[I[2]]
               u[I,1] = 22*(y*(1-y))^(3/2)
               u[I,2] = 27*(x*(1-x))^(3/2)
           end
           u
       end
       u0_manual = init_brusselator_2d(xyd_brusselator)
       prob = ODEProblem(brusselator_2d_loop,u0_manual,(0.,11.5),p)

       msol = solve(prob,TRBDF2(),saveat=0.01) # 2.771 s (5452 allocations: 65.73 MiB)

       # get variables for reshape
       Nx = floor(Int64, (x_max - x_min) / dx) + 1
       Ny = floor(Int64, (y_max - y_min) / dy) + 1

       @variables u[1:Nx,1:Ny](t)
       @variables v[1:Nx,1:Ny](t)

       t = sol[t]
       @testset "." begin
       for k in div(length(t), 2):length(t)
              solu = reshape([sol[u[(i-1)*Ny+j]][k] for i in 1:Nx for j in 1:Ny],(Nx,Ny))[2:end,2:end]
              msolu = msol.u[k][:,:,1]
              @test solu ≈ msolu rtol = 0.1

              solv = reshape([sol[v[(i-1)*Ny+j]][k] for i in 1:Nx for j in 1:Ny],(Nx,Ny))[2:end,2:end]
              msolv = msol.u[k][:,:,2]
              @test solv ≈ msolv rtol = 0.1
       end
       end

       # Nx = floor(Int64, (x_max - x_min) / dx) + 1
       # Ny = floor(Int64, (y_max - y_min) / dy) + 1

       #  @variables u[1:Nx,1:Ny](t)
       #  @variables v[1:Nx,1:Ny](t)
       #  t = sol[t]
       #   anim = @animate for k in 1:length(t)
       #          solu = reshape([sol[u[(i-1)*Ny+j]][k] for i in 1:Nx for j in 1:Ny],(Nx,Ny))
       #          #solv = reshape([sol[v[(i-1)*Ny+j]][k] for i in 1:Nx for j in 1:Ny],(Nx,Ny))
       #          heatmap(solu[2:end, 2:end], title="$(t[k])")
       #   end
       #   gif(anim, "plots/Brusselator2Dsol.gif", fps = 8)
end


text:
SUITE["transition_range"] = BenchmarkGroup()
SUITE["transition_range"]["local"] = BenchmarkGroup()
SUITE["transition_range"]["utc"] = BenchmarkGroup()

SUITE["transition_range"]["local"]["standard"] = begin
    @benchmarkable TimeZones.transition_range($(DateTime(2015,1,1)), $(tz"America/Winnipeg"), $Local)
end
SUITE["transition_range"]["local"]["non-existent"] = begin
    @benchmarkable TimeZones.transition_range($(DateTime(2015,3,8,2)), $(tz"America/Winnipeg"), $Local)
end
SUITE["transition_range"]["local"]["ambiguous"] = begin
    @benchmarkable TimeZones.transition_range($(DateTime(2015,11,1,1)), $(tz"America/Winnipeg"), $Local)
end
SUITE["transition_range"]["utc"] = begin
    @benchmarkable TimeZones.transition_range($(DateTime(2015,1,1)), $(tz"America/Winnipeg"), $UTC)
end

SUITE["interpret"] = BenchmarkGroup()
SUITE["interpret"]["local"] = BenchmarkGroup()
SUITE["interpret"]["utc"] = BenchmarkGroup()

SUITE["interpret"]["local"]["standard"] = begin
    @benchmarkable TimeZones.interpret($(DateTime(2015,1,1)), $(tz"America/Winnipeg"), $Local)
end
SUITE["interpret"]["local"]["non-existent"] = begin
    @benchmarkable TimeZones.interpret($(DateTime(2015,3,8,2)), $(tz"America/Winnipeg"), $Local)
end
SUITE["interpret"]["local"]["ambiguous"] = begin
    @benchmarkable TimeZones.interpret($(DateTime(2015,11,1,1)), $(tz"America/Winnipeg"), $Local)
end
SUITE["interpret"]["utc"] = begin
    @benchmarkable TimeZones.interpret($(DateTime(2015,1,1)), $(tz"America/Winnipeg"), $UTC)
end


text:
include("utils.jl")

@testset "ProjectiveSpace" begin
    @testset "Real" begin
        M = ProjectiveSpace(2)
        @testset "Basics" begin
            @test repr(M) == "ProjectiveSpace(2, ℝ)"
            @test representation_size(M) == (3,)
            @test manifold_dimension(M) == 2
            @test !is_point(M, [1.0, 0.0, 0.0, 0.0])
            @test !is_vector(M, [1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0])
            @test_throws DomainError is_point(M, [2.0, 0.0, 0.0], true)
            @test !is_point(M, [2.0, 0.0, 0.0])
            @test !is_vector(M, [1.0, 0.0, 0.0], [1.0, 0.0, 0.0])
            @test_throws DomainError is_vector(M, [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], true)
            @test injectivity_radius(M) == π / 2
            @test injectivity_radius(M, ExponentialRetraction()) == π / 2
            @test injectivity_radius(M, [1.0, 0.0, 0.0]) == π / 2
            @test injectivity_radius(M, [1.0, 0.0, 0.0], ExponentialRetraction()) == π / 2
        end
        types = [Vector{Float64}]
        TEST_STATIC_SIZED && push!(types, MVector{3,Float64})

        TEST_FLOAT32 && push!(types, Vector{Float32})
        basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
        @testset "Type $T" for T in types
            x = [1.0, 0.0, 0.0]
            v = [0.0, 1.0, 0.0]
            y = exp(M, x, v)
            w = [0.0, 1.0, -1.0]
            z = exp(M, x, w)
            pts = convert.(T, [x, y, z])
            test_manifold(
                M,
                pts,
                test_injectivity_radius=false,
                test_project_point=true,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    ProjectionTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                point_distributions=[Manifolds.uniform_distribution(M, pts[1])],
                tvector_distributions=[
                    Manifolds.normal_tvector_distribution(M, pts[1], 1.0),
                ],
                test_forward_diff=false,
                test_reverse_diff=false,
                basis_types_vecs=(
                    DiagonalizingOrthonormalBasis([0.0, 1.0, 2.0]),
                    basis_types...,
                ),
                basis_types_to_from=basis_types,
                test_vee_hat=false,
                retraction_methods=[
                    ProjectionRetraction(),
                    PolarRetraction(),
                    QRRetraction(),
                ],
                inverse_retraction_methods=[
                    ProjectionInverseRetraction(),
                    PolarInverseRetraction(),
                    QRInverseRetraction(),
                ],
                is_tangent_atol_multiplier=1,
                test_inplace=true,
            )
        end

        @testset "retract/inverse_retract" begin
            x = [1.0, 0.0, 0.0]
            v = [0.0, 1.0, 0.0]
            y = retract(M, x, v, ProjectionRetraction())
            v2 = inverse_retract(M, x, y, ProjectionInverseRetraction())
            @test v ≈ v2
        end

        @testset "equivalence" begin
            x = [1.0, 0.0, 0.0]
            v = [0.0, 1.0, 0.0]
            @test isapprox(M, x, -x)
            @test isapprox(M, x, exp(M, x, π * v))
            @test log(M, x, -x) ≈ zero(v)
            @test isapprox(M, -x, vector_transport_to(M, x, v, -x), -v)
        end

        @testset "Distribution MVector tests" begin
            upd_mvector = Manifolds.uniform_distribution(M, @MVector [1.0, 0.0, 0.0])
            @test isa(rand(upd_mvector), MVector)

            gtpd_mvector =
                Manifolds.normal_tvector_distribution(M, (@MVector [1.0, 0.0, 0.0]), 1.0)
            @test isa(rand(gtpd_mvector), MVector)
        end
    end

    @testset "Complex" begin
        M = ProjectiveSpace(2, ℂ)
        @testset "Basics" begin
            @test repr(M) == "ProjectiveSpace(2, ℂ)"
            @test representation_size(M) == (3,)
            @test manifold_dimension(M) == 4
            @test Manifolds.allocation_promotion_function(M, exp!, (1,)) == complex
            @test !is_point(M, [1.0 + 0im, 0.0, 0.0, 0.0])
            @test !is_vector(M, [1.0 + 0im, 0.0, 0.0, 0.0], [0.0 + 0im, 1.0, 0.0])
            @test_throws DomainError is_point(M, [1.0, im, 0.0], true)
            @test !is_point(M, [1.0, im, 0.0])
            @test !is_vector(M, [1.0 + 0im, 0.0, 0.0], [1.0 + 0im, 0.0, 0.0])
            @test !is_vector(M, [1.0 + 0im, 0.0, 0.0], [-0.5im, 0.0, 0.0])
            @test_throws DomainError is_vector(
                M,
                [1.0 + 0im, 0.0, 0.0],
                [1.0 + 0im, 0.0, 0.0],
                true,
            )
            @test_throws DomainError is_vector(
                M,
                [1.0 + 0im, 0.0, 0.0],
                [-0.5im, 0.0, 0.0],
                true,
            )
            @test injectivity_radius(M) == π / 2
            @test injectivity_radius(M, ExponentialRetraction()) == π / 2
            @test injectivity_radius(M, [1.0 + 0im, 0.0, 0.0]) == π / 2
            @test injectivity_radius(M, [1.0 + 0im, 0.0, 0.0], ExponentialRetraction()) ==
                  π / 2
        end
        types = [Vector{ComplexF64}]
        @testset "Type $T" for T in types
            x = [0.5 + 0.5im, 0.5 + 0.5im, 0]
            v = [0.0, 0.0, 1.0 - im]
            y = im * exp(M, x, v)
            w = [0.5, -0.5, 0.5im]
            z = (sqrt(0.5) - sqrt(0.5) * im) * exp(M, x, w)
            pts = convert.(T, [x, y, z])
            test_manifold(
                M,
                pts,
                test_injectivity_radius=false,
                test_project_point=true,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    ProjectionTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                test_forward_diff=false,
                test_reverse_diff=false,
                basis_types_to_from=(DefaultOrthonormalBasis(),),
                test_vee_hat=false,
                retraction_methods=[
                    ProjectionRetraction(),
                    PolarRetraction(),
                    QRRetraction(),
                ],
                inverse_retraction_methods=[
                    ProjectionInverseRetraction(),
                    PolarInverseRetraction(),
                    QRInverseRetraction(),
                ],
                is_tangent_atol_multiplier=1,
                exp_log_atol_multiplier=10.0^3,
                retraction_atol_multiplier=10.0,
            )
        end

        @testset "retract/inverse_retract" begin
            x = [0.5 + 0.5im, 0.5 + 0.5im, 0]
            v = [0.0, 0.0, 1.0 - im]
            y = retract(M, x, v, ProjectionRetraction())
            v2 = inverse_retract(M, x, y, ProjectionInverseRetraction())
            @test v ≈ v2
        end

        @testset "equivalence" begin
            x = [1.0 + 0im, 0.0, 0.0]
            v = [0.0, im, 0.0]
            s = sqrt(0.5) - sqrt(0.5) * im
            @test isapprox(M, x, s * x)
            @test isapprox(M, x, exp(M, x, π * v))
            @test log(M, x, s * x) ≈ zero(v)
            @test isapprox(M, s * x, vector_transport_to(M, x, v, s * x), s * v)
        end
    end

    @testset "Right Quaternion" begin
        M = ProjectiveSpace(2, ℍ)
        @testset "Basics" begin
            @test repr(M) == "ProjectiveSpace(2, ℍ)"
            @test representation_size(M) == (3,)
            @test manifold_dimension(M) == 8
            @test !is_point(M, Quaternion[1.0 + 0im, 0.0, 0.0, 0.0])
            @test !is_vector(
                M,
                Quaternion[1.0 + 0im, 0.0, 0.0, 0.0],
                Quaternion[0.0 + 0im, 1.0, 0.0],
            )
            @test_throws DomainError is_point(M, Quaternion[1.0, im, 0.0], true)
            @test !is_point(M, Quaternion[1.0, im, 0.0])
            @test !is_vector(
                M,
                Quaternion[1.0 + 0im, 0.0, 0.0],
                Quaternion[1.0 + 0im, 0.0, 0.0],
            )
            @test !is_vector(
                M,
                Quaternion[1.0 + 0im, 0.0, 0.0],
                Quaternion[-0.5im, 0.0, 0.0],
            )
            @test_throws DomainError is_vector(
                M,
                Quaternion[1.0 + 0im, 0.0, 0.0],
                Quaternion[1.0 + 0im, 0.0, 0.0],
                true,
            )
            @test_throws DomainError is_vector(
                M,
                Quaternion[1.0 + 0im, 0.0, 0.0],
                Quaternion[-0.5im, 0.0, 0.0],
                true,
            )
            @test injectivity_radius(M) == π / 2
            @test injectivity_radius(M, ExponentialRetraction()) == π / 2
            @test injectivity_radius(M, Quaternion[1.0 + 0im, 0.0, 0.0]) == π / 2
            @test injectivity_radius(
                M,
                Quaternion[1.0 + 0im, 0.0, 0.0],
                ExponentialRetraction(),
            ) == π / 2
        end
        types = [Vector{Quaternion{Float64}}]
        @testset "Type $T" for T in types
            x = [Quaternion(0.5, 0, 0, 0.5), Quaternion(0, 0, 0.5, 0.5), 0]
            v = [Quaternion(0), Quaternion(0), Quaternion(0.0, -0.5, -0.5, 0.0)]
            y = Quaternion(0, 0, 0, 1) * exp(M, x, v)
            w = [
                Quaternion(0.25, -0.25, 0.25, 0.25),
                Quaternion(0.25, 0.25, -0.25, -0.25),
                1,
            ]
            z = Quaternion(0.5, -0.5, 0.5, -0.5) * exp(M, x, w)
            pts = convert.(T, [x, y, z])
            test_manifold(
                M,
                pts,
                test_injectivity_radius=false,
                test_project_point=true,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    ProjectionTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                test_forward_diff=false,
                test_reverse_diff=false,
                basis_types_to_from=(DefaultOrthonormalBasis(),),
                test_vee_hat=false,
                retraction_methods=[
                    ProjectionRetraction(),
                    PolarRetraction(),
                    QRRetraction(),
                ],
                inverse_retraction_methods=[
                    ProjectionInverseRetraction(),
                    PolarInverseRetraction(),
                    QRInverseRetraction(),
                ],
                is_tangent_atol_multiplier=10,
                exp_log_atol_multiplier=10.0^3,
                retraction_atol_multiplier=10.0,
            )
        end

        @testset "retract/inverse_retract" begin
            x = [Quaternion(0.5, 0, 0, 0.5), Quaternion(0, 0, 0.5, 0.5), 0]
            v = [Quaternion(0), Quaternion(0), Quaternion(0.0, -0.5, -0.5, 0.0)]
            y = retract(M, x, v, ProjectionRetraction())
            v2 = inverse_retract(M, x, y, ProjectionInverseRetraction())
            @test v ≈ v2
        end

        @testset "equivalence" begin
            x = Quaternion[1.0 + 0im, 0.0, 0.0]
            v = Quaternion[0.0, im, 0.0]
            s = Quaternion(0.5, -0.5, 0.5, -0.5)
            @test isapprox(M, x, x * s)
            @test isapprox(M, x, exp(M, x, π * v))
            @test log(M, x, x * s) ≈ zero(v)
            @test isapprox(M, x * s, vector_transport_to(M, x, v, x * s), v * s)
        end
    end

    @testset "ArrayProjectiveSpace" begin
        M = ArrayProjectiveSpace(2, 2; field=ℝ)
        @test manifold_dimension(M) == 3
        @test repr(M) == "ArrayProjectiveSpace(2, 2; field = ℝ)"
        @test typeof(get_embedding(M)) === Euclidean{Tuple{2,2},ℝ}
        @test representation_size(M) == (2, 2)
        p = ones(2, 2)
        q = project(M, p)
        @test is_point(M, q)
        Y = [1.0 0.0; 0.0 1.1]
        X = project(M, q, Y)
        @test is_vector(M, q, X)

        M = ArrayProjectiveSpace(2, 2; field=ℂ)
        @test manifold_dimension(M) == 6
        @test repr(M) == "ArrayProjectiveSpace(2, 2; field = ℂ)"
        @test typeof(get_embedding(M)) === Euclidean{Tuple{2,2},ℂ}
        @test representation_size(M) == (2, 2)
    end
end


text:
using OptimalTransport
using Distances

import LinearAlgebra
import Random
using Test

Random.seed!(1234)

@testset "Initial map" begin
    # example from "Computational Optimal Transport"
    P = zeros(3, 3)
    OptimalTransport.initialmap!(P, [0.2, 0.5, 0.3], [0.5, 0.1, 0.4])
    @test P == [0.2 0 0; 0.3 0.1 0.1; 0 0 0.3]
end

# example taken from https://www.me.utexas.edu/~jensen/models/network/net8.html
@testset "Simple problem" begin
    C = [3 1 1000; 4 2 4; 1000 3 3]
    P, f, g = earthmover([5, 7, 3], [7, 3, 5], C)

    @test vec(sum(P; dims=1)) == [7, 3, 5]
    @test vec(sum(P; dims=2)) == [5, 7, 3]
    @test LinearAlgebra.dot(C, P) == 46
end

@testset "Total variation" begin
    histogram1 = LinearAlgebra.normalize!(rand(20), 1)
    histogram2 = LinearAlgebra.normalize!(rand(20), 1)

    C = (!==).(1:20, (1:20)')
    P, f, g = earthmover(histogram1, histogram2, C)

    @test vec(sum(P; dims=1)) ≈ histogram2
    @test vec(sum(P; dims=2)) ≈ histogram1
    @test LinearAlgebra.dot(C, P) ≈ totalvariation(histogram1, histogram2)
end


text:
import DataFrames
import PredictMD
import Test

Test.@test_throws(
    ErrorException,
    PredictMD.calculate_smote_pct_under(;pct_over=-111,minority_to_majority_ratio=0.5),
    )

Test.@test_throws(
    ErrorException,
    PredictMD.calculate_smote_pct_under(;pct_over=100,minority_to_majority_ratio=-1.5),
    )

features_df_0rows = DataFrames.DataFrame()
labels_df_0rows = DataFrames.DataFrame()
labels_df_3rows = DataFrames.DataFrame()
labels_df_3rows[:y] = [1,2,3]
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "",
                    minorityclass = "",
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "",
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_3rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    minority_to_majority_ratio = 1,
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    minority_to_majority_ratio = 1,
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_0rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    minority_to_majority_ratio = 1,
                    ),
    )
Test.@test_throws(
    ErrorException,
    PredictMD.smote(features_df_0rows,
                    labels_df_3rows,
                    Symbol[],
                    :y;
                    majorityclass = "majorityclass",
                    minorityclass = "minorityclass",
                    minority_to_majority_ratio = 1,
                    ),
    )


text:
##
### g(x, y) = |x²-y|
##
struct maxquadAL <: Regularizer
end

g(::maxquadAL, xy) = abs(xy[1]^2 - xy[2])

function prox_αg!(::maxquadAL, res, xy, γ)
    x̄, ȳ = xy[1], xy[2]

    ## Detect structure
    man = Euclidean(2)

    if ȳ < (x̄/(1+2γ))^2 - γ
        # printstyled("lower quad\n", color=:red)
        res .= [
            x̄ / (1+2γ),
            ȳ + γ
        ]
    elseif (x̄/(1-2γ))^2 + γ < ȳ
        # printstyled("upper quad\n", color=:red)
        res .= [
            x̄ / (1-2γ),
            ȳ - γ
        ]
    else
        # printstyled("Manifold!\n", color=:red)
        man = PlaneParabola()

        ## Compute prox output
        a = 16*γ^2
        b = 8*γ*(1-2γ)
        c = (1-2γ)^2
        ā = -2γ
        b̄ = γ - ȳ

        ts = roots([
            c*b̄ + x̄^2,
            b*b̄ + c*ā,
            a*b̄ + b*ā,
            a*ā
        ])
        ts = filter(c -> abs(imag(c)) < 1e-12, ts)
        length(ts) == 0 && @warn "Prox of maxquadAL troubled" ts
        ts = map(c -> real(c), ts)
        length(ts) == 0 && @warn "Prox of maxquadAL troubled" ts
        ts = filter(c -> 0 <= c <= 1, ts)
        length(ts) != 1 && @warn "Prox of maxquadAL troubled" ts
        t = first(ts)

        res .= [
            x̄ / (1+4γ*t-2γ),
            ȳ + 2γ*t - γ,
        ]

        !is_manifold_point(man, res) && @warn "prox_αg! returned point not on manifold"
    end

    return man
end

∇M_g!(::maxquadAL, M, ∇M_g, xy) = (∇M_g .= 0)
∇M_g(::maxquadAL, M, xy) = zeros(2)

∇²M_g_ξ(::maxquadAL, M, xy, ξ) = zeros(2)
∇²M_g_ξ!(::maxquadAL, M, res, xy, ξ) = (res .= 0)


text:
using Fishmetica
using Test
using DelimitedFiles

println("Test 1")
@time  include("test1.jl")

println("Test 2")
@time include("test2.jl")

println("Test 3")
@time include("test3.jl")

# write your own tests here
#@test 1 == 2

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
struct IB
    client::IBKRCore.IBClient
    host::String
    port::CxxInt
    client_id::CxxInt
    function IB(host, port, client_id)
        new(IBKRCore.IBClient(), host, port, client_id)
    end
end

function IB(; host=IB_DEFAULT_HOST, port=IB_DEFAULT_PORT, client_id=IB_DEFAULT_CLIENT_ID)
    IB(host, port, client_id)
end

function (ib::IB)(f::Function)
    IBKRCore.connect(ib.client, ib.host, ib.port, ib.client_id)
    ev = Threads.Event()
    t = Threads.@spawn process_messages_loop($(ib.client), $ev)
    wait(ev)
    try
        f()
    finally
        IBKRCore.isConnected(ib.client) && IBKRCore.disconnect(ib.client)
        # IBKRCore.issueSignal(c)
        wait(t)
    end
    return nothing
end

function process_messages_loop(client::IBKRCore.IBClient, ev::Event)
    try
        while !IBKRCore.isConnected(client) sleep(0.001) end
        notify(ev)
        while IBKRCore.isConnected(client)
            IBKRCore.processMessages(client)
        end
    finally
        IBKRCore.isConnected(client) && IBKRCore.disconnect(client)
    end
    @debug "DONE process_messages_loop"
end

# function start_request(ib::IB)
#     ev = Threads.Event()
#     push!
#
# function managed_accounts(ib::IB)
#     IBKRCore.reqManagedAccounts(c)


text:
module PipsSolver

using Plasmo
using JuMP
using Libdl
using MPI
using Distributed
using DataStructures

export pipsnlp_solve, distribute

include("PipsNlpInterface.jl")

using .PipsNlpInterface

include("distribute.jl")

end # module


text:
"""
    step!(solution::AbstractRungeKuttaSolution, problem::AbstractInitialValueProblem, solver::ExplicitRungeKuttaSolver, cache::ExplicitRungeKuttaCache)

computes a step of the [`AbstractRungeKuttaSolution`](@ref) of an [`AbstractInitialValueProblem`](@ref) using an [`ExplicitRungeKuttaSolver`](@ref).
"""
function step!(solution::AbstractRungeKuttaSolution, problem::AbstractInitialValueProblem, solver::ExplicitRungeKuttaSolver, cache::ExplicitRungeKuttaCache)
    @↓ n, v = cache
    @↓ u, t = solution
    k = solution.k isa Nothing ? cache.k : solution.k[n]
    @↓ f! = problem.rhs
    @↓ s, A, c, b = solver.tableau
    @↓ h = solver.stepsize
    v = u[n+1] # avoid allocs
    # Compute stages
    for i = 1:s
        zero!(v)
        for j = 1:i-1
            @. v += A[i,j] * k[j]
        end
        @. v = u[n] + h * v
        # @← k[i] = f(v, t[n] + h * c[i])
        f!(k[i], v, t[n] + h * c[i])
    end
    # Compute step
    zero!(v)
    for i = 1:s
        @. v += b[i] * k[i]
    end
    @. v = u[n] + h * v
    t[n+1] = t[n] + h
    return u[n+1], t[n+1]
end

"""
    step!(solution::AbstractRungeKuttaSolution, problem::AbstractInitialValueProblem, solver::ImplicitRungeKuttaSolver, cache::ImplicitRungeKuttaCache)

computes a step of the [`AbstractRungeKuttaSolution`](@ref) of an [`AbstractInitialValueProblem`](@ref) using an [`ImplicitRungeKuttaSolver`](@ref).
"""
function step!(solution::AbstractRungeKuttaSolution, problem::AbstractInitialValueProblem, solver::ImplicitRungeKuttaSolver, cache::ImplicitRungeKuttaCache)
    @↓ n, v, Δk, J = cache
    @↓ u, t = solution
    k = solution.k isa Nothing ? cache.k : solution.k[n]
    @↓ f!, Df! = problem.rhs
    @↓ s, A, c, b = solver.tableau
    @↓ ϵ, K = solver.newton
    @↓ h = solver.stepsize
    v = u[n+1] # avoid allocs
    # @← J = Df(v, u[n], t[n])
    Df!(J, v, u[n], t[n])
    Z = factorize(I - h * kron(A, J))
    # Compute stages
    zero!(k)
    for l = 1:K
        for i = 1:s
            zero!(v)
            for j = 1:s
                @. v += A[i,j] * k[j]
            end
            @. v = u[n] + h * v
            # @← Δk[i] = f(v, t[n] + h * c[i])
            f!(Δk[i], v, t[n] + h * c[i])
            @. Δk[i] -= k[i]
        end
        # temporary workaround for ldiv!(Z, Δk)
        Δk_ = vcat(Δk...)
        ldiv!(Z, Δk_)
        if norm(Δk_) < ϵ * norm(k)
            break
        end
        # temporary workaround for k .+= Δk
        for i in eachindex(k)
            L = length(k[i])
            k[i] .+= Δk_[(i-1)*L+1:i*L]
        end
    end
    # Compute step
    zero!(v)
    for i = 1:s
        @. v += b[i] * k[i]
    end
    @. v = u[n] + h * v
    t[n+1] = t[n] + h
    return u[n+1], t[n+1]
end

function step!(solution::AbstractRungeKuttaSolution, problem::AbstractInitialValueProblem, solver::ExplicitExponentialRungeKuttaSolver, cache::ExplicitExponentialRungeKuttaCache)
    @↓ n, Q, α, β, γ, E, E2 = cache
    @↓ u, t = solution
    @↓ f = problem.rhs.rhs
    @↓ h = solver.stepsize

    uₙ = u[n]
    tₙ = t[n+1] = t[n] + h

    Nu = f(uₙ, tₙ)
    aₙ = @. E2 * uₙ + Q * Nu
    Na = f(aₙ, tₙ)
    bₙ = @. E2 * uₙ + Q * Na
    Nb = f(bₙ, tₙ)
    cₙ = @. E2 * aₙ + Q * (2Nb - Nu)
    Nc = f(cₙ, tₙ)

    @. u[n+1] = E * uₙ + α * Nu + 2β * (Na + Nb) + γ * Nc
    return u[n+1], t[n+1]
end


text:
import Languages.name

"""
    title(doc)

Return the title metadata for `doc`.

See also: [`title!`](@ref), [`titles`](@ref), [`titles!`](@ref)
"""
title(d::AbstractDocument) = d.metadata.title

"""
    language(doc)

Return the language metadata for `doc`.

See also: [`language!`](@ref), [`languages`](@ref), [`languages!`](@ref)
"""
language(d::AbstractDocument) = d.metadata.language

"""
    author(doc)

Return the author metadata for `doc`.

See also: [`author!`](@ref), [`authors`](@ref), [`authors!`](@ref)
"""
author(d::AbstractDocument) = d.metadata.author

"""
    timestamp(doc)

Return the timestamp metadata for `doc`.

See also: [`timestamp!`](@ref), [`timestamps`](@ref), [`timestamps!`](@ref)
"""
timestamp(d::AbstractDocument) = d.metadata.timestamp

"""
    title!(doc, str)

Set the title of `doc` to `str`.

See also: [`title`](@ref), [`titles`](@ref), [`titles!`](@ref)
"""
function title!(d::AbstractDocument, nv::AbstractString)
    d.metadata.title = nv
end

"""
    language!(doc, lang::Language)

Set the language of `doc` to `lang`.

# Example
```julia-repl
julia> d = StringDocument("String Document 1")

julia> language!(d, Languages.Spanish())

julia> d.metadata.language
Languages.Spanish()
```

See also: [`language`](@ref), [`languages`](@ref), [`languages!`](@ref)
"""
function language!(d::AbstractDocument, nv::Language)
    d.metadata.language = nv
end

"""
    author!(doc, author)

Set the author metadata of doc to `author`.

See also: [`author`](@ref), [`authors`](@ref), [`authors!`](@ref)
"""
function author!(d::AbstractDocument, nv::AbstractString)
    d.metadata.author = nv
end

"""
    timestamp!(doc, timestamp::AbstractString)

Set the timestamp metadata of doc to `timestamp`.

See also: [`timestamp`](@ref), [`timestamps`](@ref), [`timestamps!`](@ref)
"""
function timestamp!(d::AbstractDocument, nv::AbstractString)
    d.metadata.timestamp = nv
end

"""
    titles(crps)

Return the titles for each document in `crps`.

See also: [`titles!`](@ref), [`title`](@ref), [`title!`](@ref)
"""
titles(c::Corpus) = map(d -> title(d), documents(c))

"""
    languages(crps)

Return the languages for each document in `crps`.

See also: [`languages!`](@ref), [`language`](@ref), [`language!`](@ref)
"""
languages(c::Corpus) = map(d -> language(d), documents(c))

"""
    authors(crps)

Return the authors for each document in `crps`.

See also: [`authors!`](@ref), [`author`](@ref), [`author!`](@ref)
"""
authors(c::Corpus) = map(d -> author(d), documents(c))

"""
    timestamps(crps)

Return the timestamps for each document in `crps`.

See also: [`timestamps!`](@ref), [`timestamp`](@ref), [`timestamp!`](@ref)
"""
timestamps(c::Corpus) = map(d -> timestamp(d), documents(c))

titles!(c::Corpus, nv::AbstractString) = title!.(documents(c), nv)
languages!(c::Corpus, nv::T) where {T <: Language} = language!.(documents(c), Ref(nv)) #Ref to force scalar broadcast
authors!(c::Corpus, nv::AbstractString) = author!.(documents(c), Ref(nv))
timestamps!(c::Corpus, nv::AbstractString) = timestamp!.(documents(c), Ref(nv))

"""
    titles!(crps, vec::Vector{String})
    titles!(crps, str)

Update titles of the documents in a Corpus.

If the input is a String, set the same title for all documents. If the input is a vector, set title of `i`th document to corresponding `i`th element in the vector `vec`. In the latter case, the number of documents must equal the length of vector.

See also: [`titles`](@ref), [`title!`](@ref), [`title`](@ref)
"""
function titles!(c::Corpus, nvs::Vector{String})
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        title!(d, nvs[i])
    end
end

"""
    languages!(crps, langs::Vector{Language})
    languages!(crps, lang::Language)

Update languages of documents in a Corpus.

If the input is a Vector, then language of the `i`th document is set to the `i`th element in the vector, respectively. However, the number of documents must equal the length of vector.

See also: [`languages`](@ref), [`language!`](@ref), [`language`](@ref)
"""
function languages!(c::Corpus, nvs::Vector{T}) where T <: Language
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        language!(d, nvs[i])
    end
end

"""
    authors!(crps, athrs)
    authors!(crps, athr)

Set the authors of the documents in `crps` to the `athrs`, respectively.

See also: [`authors`](@ref), [`author!`](@ref), [`author`](@ref)
"""
function authors!(c::Corpus, nvs::Vector{String})
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        author!(d, nvs[i])
    end
end

"""
    timestamps!(crps, times::Vector{String})
    timestamps!(crps, time::AbstractString)

Set the timestamps of the documents in `crps` to the timestamps in `times`, respectively.

See also: [`timestamps`](@ref), [`timestamp!`](@ref), [`timestamp`](@ref)
"""
function timestamps!(c::Corpus, nvs::Vector{String})
    length(c) == length(nvs) || throw(DimensionMismatch("dimensions must match"))
    for (i, d) in pairs(IndexLinear(), documents(c))
        timestamp!(d, nvs[i])
    end
end