text:
function labelposition(ranges, dim)
    a, b = extrema(ranges[dim])
    pos = Float32(((b - a) / 2) + a)
    axis_vec = unit(Point{N, Float32}, dim)
    normal = unit(Point{N, Float32}, mod1(dim + 1, 3))
    pos .* axis_vec .- (normal * 0.2f0)
end

function GeometryTypes.widths(x::Range)
    mini, maxi = Float32.(extrema(x))
    maxi - mini
end

function draw_axis(
        textbuffer::TextBuffer{N}, linebuffer, ranges,
        axisnames, visible, showaxis, showticks, showgrid,
        axiscolors, gridcolors, gridthickness, tickfont
    ) where N

    empty!(textbuffer); empty!(linebuffer)
    origin = Point{N, Float32}(map(minimum, ranges))
    is2d = N == 2 # sadly we need to special case on 2D case
    for i = 1:N
        axis_vec = unit(Point{N, Float32}, i)
        width = widths(ranges[i])
        stop = origin .+ (width .* axis_vec)
        if showaxis[i]
            append!(linebuffer, [origin, stop], axiscolors[i], 1.5f0)
        end
        if showticks[i]
            range = ranges[i]
            j = mod1(i + 1, N)
            tickdir = unit(Point{N, Float32}, j)
            tickdir, offset2 = if i != 2
                tickdir = unit(Vec{N, Float32}, j)
                tickdir, Float32(widths(ranges[j]) + 0.3f0) * tickdir
            else
                tickdir = unit(Vec{N, Float32}, 1)
                tickdir, Float32(widths(ranges[1]) + 0.3f0) * tickdir
            end
            if is2d
                offset2 = Vec2f0(0)
            end
            offset = -tickdir .* 0.1f0
            for tick in drop(range, 1)
                startpos = origin .+ ((Float32(tick - range[1]) * axis_vec) .+ offset) .+ offset2
                str = sprint(io-> print(io, round(tick, 2)))
                append!(textbuffer, startpos, str, tickfont[i]...)
            end
        end
        if showgrid[i]
            c = gridcolors[i]
            thickness = gridthickness[i]
            for _j = (i + 1):(i + N - 1)
                j = mod1(_j, N)
                dir = unit(Point{N, Float32}, j)
                range = ranges[j]
                for tick in drop(range, 1)
                    offset = Float32(tick - range[1]) * dir
                    append!(linebuffer, [origin .+ offset, stop .+ offset], c, thickness)
                end
            end
        end
        nametext = axisnames[i]
        # if !isempty(first(nametext))
        #     pos = labelposition(ranges, i)
        #     printat(textio, startpos, nametext...)
        # end
    end
    return
end

function axis(scene::Scene, x, y, attributes::Dict)
    axis(scene, to_node((x, y)), attributes)
end

function axis(scene::Scene, x, y, z, attributes::Dict)
    axis(scene, to_node((x, y, z)), attributes)
end

function axis(scene::Scene, ranges::Node{<: NTuple{N}}, attributes::Dict) where N
    textbuffer = TextBuffer(Point{N, Float32}(0))
    linebuffer = LinesegmentBuffer(Point{N, Float32}(0))
    scene = get_global_scene()
    attributes = axis_defaults(scene, attributes)
    tickfont = N == 2 ? :tickfont2d : :tickfont3d
    names = (
        :axisnames, :visible, :showaxis, :showticks,
        :showgrid, :axiscolors, :gridcolors, :gridthickness, tickfont
    )
    args = getindex.(attributes, names)
    lift_node(
        draw_axis,
        to_node(textbuffer), to_node(linebuffer), ranges, args...
    )
    bb = to_signal(lift_node(ranges) do ranges
        mini, maxi = Vec{N, Float32}(map(minimum, ranges)), Vec{N, Float32}(map(maximum, ranges))
        mini3d, w = Vec3f0(to_nd(mini, Val{3}, 0f0)), Vec3f0(to_nd(maxi .- mini, Val{3}, 0f0))
        HyperRectangle(mini3d, w)
    end)
    linebuffer.robj.boundingbox = bb
    viz = Context(linebuffer.robj, textbuffer.robj)
    insert_scene!(scene, :axis, viz, attributes)
end


text:
# Generates a batch of plots for a MonteCarlo run.

using TrajOptSOCPs

include("plotTrajectory.jl")
include("plotTrajectory3D.jl")
include("plotConstraintViolation_MonteCarlo.jl")

function batchPlotMonteCarlo(trajList, dynVio, constVio, nDim::Int64,
                            penalty::Float64 = 1.0, maxThrust::Float64 = 0.0,
                            maxAngleDeg::Float64 = 0.0)

    # First and foremost, we want to plot the trajectories
    ptList = [getParseTrajectory(traj, nDim) for traj in trajList]

    if nDim == 2
        pltTraj = plotTrajPos2D_Multiple(ptList)
        display(pltTraj)
    elseif nDim == 3
        pltTraj = plotTrajPos3D_Multiple(ptList)
        display(pltTraj)
    end

    pltVio = plotDynConstVio(dynVio, constVio, true)

    pltList = [pltTraj, pltVio]
    pltLabels = ["Trajectory", "Violation"]

    # Return the result as a dictionary
    return Dict(zip(pltLabels, pltList))
end

ext: jl

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
"""Generate quadtree root from .lax header."""
function quadtree(header::LaxQuadtreeHeader)
    root = Cell(SVector(header.min_x, header.min_y), SVector(header.max_x - header.min_x, header.max_y - header.min_y), Vector{UnitRange{Integer}}())
end

"""Add data to a quadtree cell with given index."""
function quadtree!(root::RegionTrees.Cell, index::Integer, data::Vector{UnitRange{T}}) where T <: Integer
    cell = quadtree!(root, index)
    cell.data = data
    return cell
end

"""Split a quadtree until the requested index is reached."""
function quadtree!(root::RegionTrees.Cell, index::Integer)
    places = zlevels(index)
    places = reverse(places, dims=1)  # bottom up

    for place in places
        (place != 0) && (isleaf(root)) && (split!(root))
        root = getindex(root, place)
    end
    return root
end

"""For a given Quadtree index, return one of [1,2,3,4] for each level in the quadtree.
Levels are defined as follows, where the index is the number.
0
1 2 3 4
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 - 84 etc
"""
function zlevels(index::Integer)
    places = Array{Integer,1}()

    # Determine places until root is reached
    while index > 4

        # index one level up in the quadtree
        parent_cell = (index - 1) >>> 2

        # place is position in Z order of four cells
        place = index - (parent_cell << 2)

        push!(places, place)
        index = parent_cell
    end

    push!(places, index)
    return places
end

"""Convert place in Z order to RegionTrees index."""
function Base.getindex(cell::Cell, place::Integer)
    # 3 4
    # 1 2
    place == 4 && return cell[2,2]
    place == 3 && return cell[1,2]
    place == 2 && return cell[2,1]
    place == 1 && return cell[1,1]
    place == 0 && return cell  # root cell
    error("Place can't be higher than 4.")
end


text:
export SteadyStateProblem

"""
    SteadyStateProblem([T=STD_REAL_PREC], ℒ, operators=true, variance=true)
    SteadyStateProblem([T=STD_REAL_PREC], H, J, operators=true, variance=true)

Returns the problem targeting the steady-state of the Liouvillian `ℒ` through
a minimisation of the cost function 𝒞 = Tr[ρ'ℒ'ℒ ρ]/Tr[ρ'ρ] .

If `variance=true` the cost function is sampled by evaluating ℒρ, which has
better convergence properies. The sampling is performed on ℒ'ℒ ρ otherwise.

See appendix of https://arxiv.org/abs/1902.10104 for more info.

If `operators=true` a memory efficient representation of the hamiltonian is used,
resulting in less memory consuption but higher CPU usage. This is needed for lattices
bigger than a certain threshold.
To use this feature, the hamiltonian must be provided as a `GraphOperator` object,
from `QuantumLattices` package.
"""
SteadyStateProblem(args...; kwargs...) = SteadyStateProblem(STD_REAL_PREC, args...; kwargs...)

# Dispatched if called with Hamiltonian and Jump operators
function SteadyStateProblem(T::Type{<:Number}, H::DataOperator, J::AbstractVector,
                            operators=false; kwargs...)
    if operators
        @warn """Converting QuantumOptics' DataOperators to memory-efficient Local Operators is not possible.
        Setting `operators=false`.
        To suppress this warning, don't request `operators=true`.
        """
        operators = false
    end
    return _build_steadystate_sparsemat_problem(T, H, J; kwargs...)
end

# Dispatched when called with the object for the whole liouvillian
function SteadyStateProblem(T::Type{<:Number}, ℒ; operators=true, variance=true, kwargs...)
    base = basis(ℒ)

    if operators
        if !variance
            throw("Can't use operators=true and variance=false. Operators are not
                   compatible with non-variance minimization.")
        end

        return LRhoKLocalOpProblem(T, ℒ)
    else # not operators
        H = SparseOperator(hamiltonian(ℒ))
        J = jump_operators(ℒ)
        _build_steadystate_sparsemat_problem(T, H, J; variance=variance, kwargs...)
    end
end

# Builds the problem if not using KLocal operators
function _build_steadystate_sparsemat_problem(T::Type{<:Number}, H, J; variance=false, fullmatrix=false)
    if variance
        if fullmatrix
            return LRhoSparseSuperopProblem(T, H, J)
        else
            return LRhoSparseOpProblem(T, H, J)
        end
    else # not variance
        if fullmatrix
            return LdagLSparseSuperopProblem(T, H, J)
        else
            return LdagLSparseOpProblem(T, H, J)
        end
    end
end


text:
# Promotion rules for "promiscuous" types like Intervals and SIUnits, which both
# simply wrap any Number, are often ambiguous. That is, which type should "win"
# -- is the promotion between an SIUnit and an ClosedInterval an SIQuantity{ClosedInterval}
# or is it an ClosedInterval{SIQuantity}? For our uses in AxisArrays, though, we can
# sidestep this problem by making Intervals *not* a subtype of Number. Then in
# order for them to plug into the promotion system, we *extend* the promoting
# operator behaviors to Union{Number, ClosedInterval}. This way other types can
# similarly define their own extensions to the promoting operators without fear
# of ambiguity -- there will simply be, e.g.,
#
# f(x::Number, y::Number) = f(promote(x,y)...) # in base
# f(x::Union{Number, ClosedInterval}, y::Union{Number, ClosedInterval}) = f(promote(x,y)...)
# f(x::Union{Number, T}, y::Union{Number, T}) = f(promote(x,y)...)
#
# In this way, these "promiscuous" types will never interact unless explicitly
# made subtypes of Number or otherwise defined with knowledge of eachother. The
# downside is that Intervals are not as useful as they could be; they really
# could be considered as <: Number themselves. We do this in general for any
# supported Scalar:
const Scalar = Union{Number, Dates.AbstractTime}
Base.promote_rule(::Type{ClosedInterval{T}}, ::Type{T}) where {T<:Scalar} = ClosedInterval{T}
Base.promote_rule(::Type{ClosedInterval{T}}, ::Type{S}) where {T,S<:Scalar} = ClosedInterval{promote_type(T,S)}

import Base: isless, <=, >=, ==, +, -, *, /, ^, //
# TODO: Is this a total ordering? (antisymmetric, transitive, total)?
isless(a::ClosedInterval, b::ClosedInterval) = isless(a.right, b.left)
# The default definition for <= assumes a strict total order (<=(x,y) = !(y < x))
<=(a::ClosedInterval, b::ClosedInterval) = a.left <= b.left && a.right <= b.right
+(a::ClosedInterval) = a
+(a::ClosedInterval, b::ClosedInterval) = ClosedInterval(a.left + b.left, a.right + b.right)
-(a::ClosedInterval) = ClosedInterval(-a.right, -a.left)
-(a::ClosedInterval, b::ClosedInterval) = a + (-b)
for f in (:(*), :(/), :(//))
    # For a general monotonic operator, we compute the operation over all
    # combinations of the endpoints and return the widest interval
    @eval function $(f)(a::ClosedInterval, b::ClosedInterval)
        w = $(f)(a.left, b.left)
        x = $(f)(a.left, b.right)
        y = $(f)(a.right, b.left)
        z = $(f)(a.right, b.right)
        ClosedInterval(min(w,x,y,z), max(w,x,y,z))
    end
end

# Extend the promoting operators to include Intervals. The comparison operators
# (<, <=, and ==) are a pain since they are non-promoting fallbacks that call
# isless, !(y < x) (which is wrong), and ===. So implementing promotion with
# Union{T, ClosedInterval} causes stack overflows for the base types. This is safer:
for f in (:isless, :(<=), :(>=), :(==), :(+), :(-), :(*), :(/), :(//))
    # We don't use promote here, though, because promotions can be lossy... and
    # that's particularly bad for comparisons. Just make an interval instead.
    @eval $(f)(x::ClosedInterval, y::Scalar) = $(f)(x, y..y)
    @eval $(f)(x::Scalar, y::ClosedInterval) = $(f)(x..x, y)
end

# And, finally, we have an Array-of-Structs to Struct-of-Arrays transform for
# the common case where the interval is constant over many offsets:
struct RepeatedInterval{T,S,A} <: AbstractVector{T}
    window::ClosedInterval{S}
    offsets::A # A <: AbstractVector
end
RepeatedInterval(window::ClosedInterval{S}, offsets::A) where {S,A<:AbstractVector} =
    RepeatedInterval{promote_type(ClosedInterval{S}, eltype(A)), S, A}(window, offsets)
Base.size(r::RepeatedInterval) = size(r.offsets)
Base.IndexStyle(::Type{<:RepeatedInterval}) = IndexLinear()
Base.getindex(r::RepeatedInterval, i::Int) = r.window + r.offsets[i]
+(window::ClosedInterval, offsets::AbstractVector) = RepeatedInterval(window, offsets)
+(offsets::AbstractVector, window::ClosedInterval) = RepeatedInterval(window, offsets)
-(window::ClosedInterval, offsets::AbstractVector) = RepeatedInterval(window, -offsets)
-(offsets::AbstractVector, window::ClosedInterval) = RepeatedInterval(-window, offsets)

# As a special extension to intervals, we allow specifying Intervals about a
# particular index, which is resolved by an axis upon indexing.
struct IntervalAtIndex{T}
    window::ClosedInterval{T}
    index::Int
end
atindex(window::ClosedInterval, index::Integer) = IntervalAtIndex(window, index)

# And similarly, an AoS -> SoA transform:
struct RepeatedIntervalAtIndexes{T,A<:AbstractVector{Int}} <: AbstractVector{IntervalAtIndex{T}}
    window::ClosedInterval{T}
    indexes::A # A <: AbstractVector{Int}
end
atindex(window::ClosedInterval, indexes::AbstractVector) = RepeatedIntervalAtIndexes(window, indexes)
Base.size(r::RepeatedIntervalAtIndexes) = size(r.indexes)
Base.IndexStyle(::Type{<:RepeatedIntervalAtIndexes}) = IndexLinear()
Base.getindex(r::RepeatedIntervalAtIndexes, i::Int) = IntervalAtIndex(r.window, r.indexes[i])


text:
#NOTE: We may want to convert this to monthly at some point

const ALL_COLS = [:name, :charitytype, :ein, :fisyr, :fileyear, :invrealized,
  :invrealizeddirect, :totrev, :totexp,
  :invrealizedindirect, :name, :netassets,
  :bookassets, :bookliabilities, :contributions,
  :pfexpbook, :pfexpcharity, :totinv, :category,:datefiscal,
  :returnamt, :returncheckamt,:netincome, :naics, :affiliationcode]

const CATEGORY_CONSOLIDATED_THRESHOLD = 50_000

#=const :Arts, :Education, :Environmental, :AnimalWelfare, :Health, :MentalHealth,
:MedicalDiseases, :MedicalResearch, :CrimeAndLegal, :Employment, :FoodAndAgriculture,
:Housing, :PublicSafety, :Recreation, :YouthDevelopment, :HumanServices, :ForeignAffairs,
:SocialAction, :CommunityImprovement, :GrantmakingFoundations, :ScienceandTechResearch,
:SocialScienceResearch, :SocietyBenefitAndMultipurpose, :ReligionRelated,
:MutualBenefitOrg, :Unknown, :Other=#

#processes the returns in a reasonable manner
function processReturns!(data::NCCSData; maxReturnCutoff::Float64 = MAX_RETURN_CUTOFF,
  minReturnCutoff::Float64 = MIN_RETURN_CUTOFF,
  winsorizeprop::Float64 = WINSORIZE_PROP)

  local N::Int = size(data.df,1) #length of dataframe
  local revenueCutoff::Float64 = REVENUE_CUTOFF
  local largeDifferenceCutoff::Float64 = LARGE_DIFFERENCE_CUTOFF #max discrepenacy between methods

  println("data.df N1: $N")
  println("data.df pf N1: $(sum(data.df[:charitytype] .== :pf))")
  #filter out organizations which are too small
  data.df[:tokeep] = trues(size(data.df,1))

  #dfByein::GroupedDataFrame = groupby(data.df, :ein)
  for subdf::SubDataFrame ∈ groupby(data.df, :ein)
    subN::Int = size(subdf,1)
    maxRevenue::Float64 = maximum(collect(skipmissing(subdf[:totrev])))
    if (maxRevenue < revenueCutoff)
      subdf[:tokeep] = falses(subN)
    end

    if subN ≠ length(unique(subdf[:fisyr])) # check if we have duplicate filings
      deleteOrg::Bool = false
      subdfByfisyr = groupby(subdf, :fisyr)
      for ssubdf::SubDataFrame ∈ subdfByfisyr
        nrows::Int = size(ssubdf, 1)
        if nrows > 1 # if we are at the org-year with the dup
          fileYears::Vector{Int} = unique(ssubdf[:fileyear])
          if length(fileYears) == nrows #if we can reconcile by fileYear
            maxFileYear::Int = maximum(fileYears) # keep the newest version
            ssubdf[ssubdf[:fileyear] .≠ maxFileYear, :tokeep] = false
          else
            deleteOrg = true #If we cannot correct, throw out the org
          end
        end
      end
      if deleteOrg #delete the org if we need to
        subdf[:tokeep] = false
      end
    end
  end

  #delete the invalid rows
  filter!(r->(r[:tokeep]),data.df)
  #data.df = data.df[data.df[:tokeep] .== true, :]
  N = size(data.df,1)

  println("data.df N2: $N")

  data.df[:return] = Vector{Union{Float64,Missing}}(missings(N))
  data.df[:lreturn] = Vector{Union{Float64,Missing}}(missings(N))

  #data.df[:returncheckamt] = Vector{Union{Float64,Missing}}(missings(N))
  data.df[:returncheck] = Vector{Union{Float64,Missing}}(missings(N))

  data.df[:largeDifference] = Vector{Union{Bool,Missing}}(missings(N))

  data.df[:lagnetassets] = Vector{MFloat64}(missing, N)
  data.df[:lag2netassets] = Vector{MFloat64}(missing, N)
  data.df[:lagcontributions] = Vector{MFloat64}(missing, N)

  sort!(data.df, [:ein, :fisyr])

  #double check we have no duplicate entries
  #@assert size(data.df,1) == size(unique(data.df[[:ein,:fisyr]]), 1)
  #ctr::Int = 0
  #dfByein = groupby(data.df, :ein)
  for subdf::SubDataFrame ∈ groupby(data.df, :ein)
    subN = size(subdf,1)

    if (subdf[1,:charitytype] == :co) || (subdf[1,:charitytype] == :pc)

      for i ∈ 2:subN

        if subdf[i-1, :fisyr] == subdf[i, :fisyr] - 1
          subdf[i, :lagnetassets] = subdf[i-1, :netassets]
          subdf[i, :lagcontributions] = subdf[i-1, :contributions]
          (i≥3) && (subdf[i, :lag2netassets] = subdf[i-2, :netassets])

          #NOTE: returns = unrealized + inv income s.t. unrealized = Δassets - netincome
          #write down 2 different methods
          returnCandidateAmt::MFloat64 = (subdf[i, :netassets] .- subdf[i, :lagnetassets]
            ) .+ subdf[i, :invrealizeddirect] .- subdf[i, :netincome]
          returnCandidate::MFloat64 = returnCandidateAmt / subdf[i, :lagnetassets]

          returnCheckCandidateAmt::MFloat64 = (subdf[i, :netassets] .- subdf[i, :lagnetassets]
            ) .+ subdf[i, :invrealizedindirect] .- subdf[i, :netincome]
          returnCheckCandidate::MFloat64 = returnCheckCandidateAmt / subdf[i, :lagnetassets]

          if ((!ismissing(returnCandidate)) && (returnCandidate < maxReturnCutoff) &&
              (returnCandidate > minReturnCutoff) &&
              (abs(returnCandidate - returnCheckCandidate) ≤ largeDifferenceCutoff)
               && abs(subdf[i, :totinv]) ≥ 1.) #make sure we have valid numbers here
            subdf[i, :returnamt] = returnCandidateAmt
            subdf[i, :return] = returnCandidate

            subdf[i, :returncheckamt] = returnCheckCandidateAmt
            subdf[i, :returncheck] = returnCheckCandidate
          end
        end
      end
    elseif subdf[1,:charitytype] == :pf
      for i ∈ 2:subN

        if subdf[i-1, :fisyr] == subdf[i, :fisyr] - 1
          subdf[i, :lagnetassets] = subdf[i-1, :netassets]
          subdf[i, :lagcontributions] = subdf[i-1, :contributions]
          (i≥3) && (subdf[i, :lag2netassets] = subdf[i-2, :netassets])

          #NOTE: returns = unrealized + inv income s.t. unrealized = Δassets - netincome
          returnCandidateAmt::MFloat64 = subdf[i, :returnamt]
          returnCandidate::MFloat64 = returnCandidateAmt / subdf[i, :lagnetassets]

          returnCheckCandidateAmt::MFloat64 = subdf[i, :returncheckamt]
          returnCheckCandidate::MFloat64 = returnCheckCandidateAmt / subdf[i, :lagnetassets]

          if ( (!ismissing(returnCandidate)) && (returnCandidate < maxReturnCutoff) &&
              (returnCandidate > minReturnCutoff) &&
              (!ismissing(returnCheckCandidate)) &&
              abs(returnCandidate - returnCheckCandidate) ≤ largeDifferenceCutoff)# && abs(subdf[i, :totinv]) ≥ 1. #make sure we have valid numbers here
            subdf[i, :returnamt] = returnCandidateAmt
            subdf[i, :return] = returnCandidate
            #subdf[i, :lreturn] = log(1.0+returnCandidate)

            subdf[i, :returncheckamt] = returnCheckCandidateAmt
            subdf[i, :returncheck] = returnCheckCandidate
          #=elseif ctr < 50
            ctr+=1
            println("###\nreturn: $(returnCandidate)")
            println("returnchk: $(returnCheckCandidate)")
            println("returnamt: $(returnCandidateAmt)")
            println("returnchk: $((returnCheckCandidateAmt))")=#
          end
        end
      end
    end
  end

  #winsorize!(data.df[:return], prop=winsorizeprop)
  data.df[:lreturn] = (log).(1.0 .+ data.df[:return])
  #NOTE: Already sorted, so don't need sort!(data.df, [:ein, :fisyr])

  return data
end

#cleans the wealth if possible
function processWealth!(data::NCCSData)

  data.df[:lnetassets] = (log).(data.df[:netassets])
  data.df[:ladjnetassets] = (log).(data.df[:adjnetassets])

  return data
end

#data processing related to the categorization
function processCategories!(data::NCCSData;
  minPointsPerCategory::Int = MIN_POINTS_PER_CATEGORY,
  categoryConsolidatedThreshold::Int = CATEGORY_CONSOLIDATED_THRESHOLD)::Nothing

  data.df[:categoryname] = (s::MSymbol->((ismissing(s) || (!haskey(category, s))) ?
    missing : category[s])).(Vector{MSymbol}(data.df[:category]))

  data.df[:tokeep] .= true
  data.df[:categoryfiltered] = deepcopy(data.df[:category])
  data.df[:categorynamefiltered] = deepcopy(data.df[:categoryname])

  dfBycategoryfisyr::GroupedDataFrame = groupby(data.df, [:category, :fisyr])
  for subdf::SubDataFrame ∈ dfBycategoryfisyr
    if size(subdf,1) < minPointsPerCategory
      subdf[:categoryfiltered] = :ZZ
      subdf[:categorynamefiltered] = category[:ZZ]
    end
  end

  data.df[:categoryconsolidated] = deepcopy(data.df[:categoryfiltered])
  data.df[:categorynameconsolidated] = deepcopy(data.df[:categorynamefiltered])

  dfBycategorynamefiltered::GroupedDataFrame = groupby(data.df, :categorynamefiltered)
  for subdf::SubDataFrame ∈ dfBycategorynamefiltered
    if size(subdf,1) < categoryConsolidatedThreshold
      subdf[:categoryconsolidated] = :ZZ
      subdf[:categorynameconsolidated] = category[:ZZ]
    end
  end

  data.df[:naics2] = (s::MString->(!ismissing(s) ? Symbol(s[1:2]) : missing)).(data.df[:naics])
  data.df[:naics2name] = (s::MSymbol->(naicsCodes[s])).(data.df[:naics2])
  for subdf::SubDataFrame ∈ groupby(data.df, :categoryname)
    cname::MSymbol = subdf[1,:categoryname]
    subdf[:naics2name] = (s::MSymbol->
      (ismissing(s) || s == :Unknown) ? ntee2naics2names[cname] : s).(subdf[:naics2name])
  end

  return nothing
end

function filternoncentral!(data::NCCSData)::Nothing

  #conform the common affiliation codes
  data.df[:affiliationcode2] = (s::MSymbol->affiliationcodes[s]).(
    Vector{MSymbol}(data.df[:affiliationcode]))

  #filter out non-central orgs
  filter!(r->conformedaffiliationcodes[r[:affiliationcode2]]==:central, data.df)

  return nothing
end

#adjusts for inflation
#currently just does wealth, eventually can use this for return fields
function adjustInflation!(data::NCCSData;
  fredPath::String = FRED_PATH,
  fredinflationfile::String = FRED_INFLATION_FILE,
  inflationseries::Symbol = :CPI,
  baseyear::Int = BASE_YEAR,
  freddategormat::String = FRED_DATE_FORMAT,
  adjustCutoff::Bool = true,
  revenuelowerbound::Float64 = REVENUE_LOWER_BOUND,
  wealthlowerbound::Float64 = WEALTH_LOWER_BOUND,
  lowerezrevenuelimit::Float64 = LOWER_EZ_REVENUE_LIMIT,
  prefilterlowerbound::Float64 = PREFILTER_REVENUE_LOWER_BOUND)::Nothing

  #read in the inflation info
  freddf::DataFrame = CSV.read("$fredPath\\$fredinflationfile.csv")

  #make a dictionary
  cpiindex::Dict =
    Dict(Dates.year(Date(freddf[i,:DATE], freddategormat))=>freddf[i, inflationseries]
      for i ∈ 1:size(freddf,1))

  #get vectors for performance purpuses
  nrows::Int = size(data.df, 1)
  basevalue::Float64 = cpiindex[baseyear]
  #=netassets::Vector{Float64} = data.df[:netassets]
  totrev::Vector{Float64} = data.df[:totrev]
  adjnetassets =
  adjtotrev =
  fisyr::Vector{Int} = data.df[:fisyr]
  @inbounds @simd for i ∈ 1:nrows
    if haskey(cpiindex, fisyr[i])
      adjnetassets[i] = netassets[i] * (basevalue / cpiindex[fisyr[i]])
      adjtotrev[i] = totrev[i] * (basevalue / cpiindex[fisyr[i]])
    end
  end

  data.df[:adjnetassets] = adjnetassets
  data.df[:adjtotrev] = adjtotrev=#

  data.df[:adjnetassets] = Vector{MFloat64}(undef, nrows)
  data.df[:adjtotrev] = Vector{MFloat64}(undef, nrows)
  data.df[:tokeep] = trues(nrows)

  for subdf ∈ groupby(data.df, :fisyr)
    currentfisyr::Int = subdf[1,:fisyr]

    #make sure its a valid fisyr and set the inflation factors
    current2basefactor::Float64 = (haskey(cpiindex, currentfisyr) ?
      basevalue / cpiindex[currentfisyr] : 0.0)
    base2currentfactor::Float64 = (haskey(cpiindex, currentfisyr) ?
      cpiindex[currentfisyr] / basevalue : 0.0)

    wealthcutoff::Float64 = wealthlowerbound * base2currentfactor
    revenuecutoff::Float64 = revenuelowerbound * base2currentfactor
    (lowerezrevenuelimit ≥ revenuecutoff) &&  @warn (
      "$currentfisyr revenue limit ($revenuecutoff) below IRS minimum" *
      "$(lowerezrevenuelimit). Expect selection bias.")

    (prefilterlowerbound ≥ revenuecutoff) && @warn (
      "$currentfisyr revenue limit ($revenuecutoff) below inflation adjusted revenue limit" *
      "$(prefilterlowerbound). Expect moving limit in real terms.")

    #adjust assets and revenue to the base year
    subdf[:adjnetassets] .= (f::MFloat64->f*current2basefactor).(subdf[:netassets])
    subdf[:adjtotrev] .= (f::MFloat64->f*current2basefactor).(subdf[:totrev])

    #delete rows below the limit
    subdf[:tokeep] .= ((fna::MFloat64, ftr::MFloat64) ->
      (!ismissing(fna)) && (!ismissing(ftr)) &&
      (wealthcutoff ≤ fna) && (revenuecutoff ≤ ftr)).(
      subdf[:netassets], subdf[:totrev])

#    println("$currentfisyr Minimum assets: $(minimum(subdf[subdf[:tokeep],:totrev])) " *
#      "Adj: $(minimum(subdf[subdf[:tokeep],:adjtotrev]))")
#    println("$currentfisyr Minimum assets: $(minimum(subdf[subdf[:tokeep],:netassets])) " *
#      "Adj: $(minimum(subdf[subdf[:tokeep],:adjnetassets]))")
  end

  data.df = data.df[data.df[:tokeep],:]
  #=if adjustCutoff #this is to account for the cutoff inflating forward
    inflatedrevenuecutoff::Float64 = REVENUE_LOWER_BOUND * basevalue/cpiindex[minimum(data.df[:fisyr])]
    println("Inflated revenue cutoff of $inflatedRevenueCutoff")
    data.df = data.df[
      (f::MFloat64-> (!ismissing(f)) && (f ≥ inflatedrevenuecutoff)).(data.df[:adjtotrev]), :]

    inflatedwealthcutoff::Float64 = WEALTH_LOWER_BOUND * basevalue/cpiindex[minimum(data.df[:fisyr])]
    data.df = data.df[
      (f::MFloat64-> (!ismissing(f)) && (f ≥ inflatedwealthcutoff)).(data.df[:adjnetassets]), :]
  end=#

  return nothing
end

#creates quantiles given a benchmark
function relativeperformance!(data::NCCSData, period::Int, lbenchmark::NSymbol;
  sorted::Bool = false, lretsym1y::NSymbol = :lreturn)::Nothing
  (!sorted) && sort!(data.df, [:ein, :fisyr])

  local N::Int = size(data.df,1)

  #first establish the field names
  lbenchstring::String = isnothing(lbenchmark) ? "" : string("net", lbenchmark)

  retstring1y::String = string(something(lretsym1y, "zero"))

  local retsym = Symbol("$retstring1y", lbenchstring, period, "yr")
  local Lretsym = Symbol("L$(period)$retstring1y", lbenchstring, period, "yr")
  local percentilesym::Symbol = Symbol("p$retstring1y", lbenchstring, period, "yr")
  local percentilebytypesym::Symbol = Symbol("p$retstring1y", lbenchstring, period, "yrtype")
  local Lpercentilesym::Symbol = Symbol("L$(period)p$retstring1y", lbenchstring, period, "yr")
  local Lpercentilebytypesym::Symbol = Symbol("L$(period)p$retstring1y", lbenchstring, period, "yrtype")

  #this will be useful
  local uniqueyrs::Vector{Int} = collect(
    minimum(data.df[:fisyr]):maximum(data.df[:fisyr]))#unique(data.df[:fisyr])

  #pre-allocate
  data.df[retsym] = Vector{MFloat64}(undef, N)
  data.df[Lretsym] = Vector{MFloat64}(undef, N)
  data.df[percentilesym] = Vector{MFloat64}(undef, N)
  data.df[percentilebytypesym] = Vector{MFloat64}(undef, N)
  data.df[Lpercentilesym] = Vector{MFloat64}(undef, N)
  data.df[Lpercentilebytypesym] = Vector{MFloat64}(undef, N)
  data.df[:ltempreturn] = Vector{MFloat64}(undef, N)

  #compute the excess returns, net of a benchmark if there is one
  data.df[:ltempreturn] .= isnothing(lretsym1y) ? zeros(N) : data.df[lretsym1y]
  data.df[:ltempreturn] .-= isnothing(lbenchmark) ? 0.0 : data.df[lbenchmark]

  #split the dataframe and compute the cumulative reutrns
  for subdf::SubDataFrame ∈ groupby(data.df, :ein)
    for i ∈ (1+period):size(subdf,1)
      #make sure we are covering the correct amount of time
      if subdf[i,:fisyr] == subdf[i-period,:fisyr] + period
        subdf[i, retsym] = sum(subdf[collect((i-(period-1)):i), :ltempreturn])
      end

    end
  end

  #split the dataframe and rank all of the returns
  for subdf::SubDataFrame ∈ groupby(data.df, :fisyr)

    F = ecdf(collect(skipmissing(subdf[retsym])))

    if isfinite(F(0.0)) #we can shortcut this loop if there are no valid ecdf values
      subdf[percentilesym] .= (f::MFloat64->ismissing(f) ? missing : F(f)).(subdf[retsym])

      #now rank within types
      for ssubdf::SubDataFrame ∈ groupby(subdf, :charitytype)
        F = ecdf(collect(skipmissing(ssubdf[retsym])))
        ssubdf[percentilebytypesym] .= (f::MFloat64->ismissing(f) ? missing : F(f)).(ssubdf[retsym])
      end
    end
  end

  #finally record the lag quantile of one period where available
  for subdf::SubDataFrame ∈ groupby(data.df, :ein)
    for i ∈ (1+2 * period):size(subdf,1)
      #make sure we are covering the correct amount of time
      if subdf[i,:fisyr] == subdf[i-period,:fisyr] + period

        #copy in the lag quantile
        subdf[i, Lretsym] = subdf[i - period,retsym]
        subdf[i, Lpercentilesym] = subdf[i - period,percentilesym]
        subdf[i, Lpercentilebytypesym] = subdf[i - period,percentilebytypesym]
      end
    end
  end

  return nothing
end

#ranks performance as a quantile
#NOTE: Returns can still be mis-aligned by year
function relativeperformance!(data::NCCSData; sorted::Bool = false,
    persistenceperiods::Vector{Int} = PERSISTENCE_PERIODS,
    lbenchmarks::Vector{NSymbol} = LBENCHMARKS,
    lbmxfields = LBM_XFIELDS)::Nothing

  for s ∈ lbenchmarks
    #need to build up column names so we don't overwrite anything

    #different procedure if no benchmark or we are not beta adjusting
    if isnothing(s)
      (i::Int->relativeperformance!(data, i, s)).(persistenceperiods)
    else
      regressbynonprofit!(data, bmfield=s, YField=:lreturn, XFields=lbmxfields[s])

      (i::Int->relativeperformance!(data, i, s)).(persistenceperiods)
      #(i::Int->relativeperformance(data, i, predlbenchmark,
      #  lretsym1y = nothing)).(persistenceperiods)
    end
  end

  return nothing
end

#lags a variable in the time series
function lagwithin!(data::NCCSData, target::Symbol;
  period::Symbol = :fisyr,
  lags::Int=1,
  laggedname = lags ≠ 1 ? Symbol(:L, lags, target) : Symbol(:L, target),
  group=:ein,
  sorted::Bool=false)::Nothing

  if !sorted
    sort!(data.df, [group, period])
  end

  #pre-allocate the space
  data.df[laggedname] = Vector{Union{eltype(data.df[target]), Missing}}(undef, size(data.df, 1))

  for subdf ∈ groupby(data.df, group)
    Nsub::Int = size(subdf, 1)

    #only run these routines if we need to
    if Nsub > lags
      for i ∈ (lags+1):Nsub #iterate for all values
        fisyr::Int = subdf[i, :fisyr]

        #first see if we can lag the easy way
        if subdf[i-lags, :fisyr] == fisyr - lags
          subdf[i,laggedname] = subdf[i-lags, target]
        else #now check the hard way (finding the lagged year)
          location::NInt = findfirst(isequal(fisyr-lags), subdf[1:(i-1),:fisyr])

          #record the lagged value if it is available
          if !isnothing(location)
            subdf[i,laggedname] = subdf[location, target]
          end
        end
      end
    end
  end

  return nothing
end

#creates a differenced column
function differencewithin!(data::NCCSData, target::Symbol;
  period::Symbol = :fisyr,
  differencedname = Symbol(:D, target),
  group=:ein,
  sorted::Bool=false,
  createlag::Bool=true,
  deletelag::Bool=true,
  laggedname::Symbol = deletelag ? Symbol(:L, target, :_temp) : Symbol(:L, target)
  )::Nothing

  createlag && lagwithin!(data, target,
    period=period,
    laggedname=laggedname,
    group=group,
    sorted=sorted)

  data.df[differencedname] = data.df[target] .- data.df[laggedname]

  deletelag && deletecols!(data.df, laggedname)

  return nothing

end

#does an initial filter to get down the file size
function processAll!(data::NCCSData; allcols::Vector{Symbol} = ALL_COLS,
  revenueCutoff::Float64 = REVENUE_LOWER_BOUND,
  wealthCutoff::Float64 = WEALTH_LOWER_BOUND,
  refreshoutsidedata::Bool = false,
  filternoncentral::Bool = true,
  refreshopen990data = true
  )

  data = removeFields!(data, setdiff(names(data.df), allcols))
  data.df[:rows] = collect(1:size(data.df,1))
  println("N0: $(size(data.df,1))")
  deleterows!(data.df, data.df[(!).(completecases(data.df[[:fisyr, :netassets, :totrev, :totexp]])),:rows])

  adjustInflation!(data)
  processReturns!(data)
  processWealth!(data)
  processCategories!(data)
  mergeoutsidedatasources!(data, refreshoutsidedata=refreshoutsidedata)
  relativeperformance!(data, sorted=false) #NOTE: Assumes already sorted
  mergeopen990data!(data, refreshopen990data=refreshopen990data)

  println("N3: $(size(data.df,1)) N3-Return: $(sum((!ismissing).(data.df[:return])))")
  filternoncentral && filternoncentral!(data)
  filternoncentral && println("Filtered out intermediate organziations.")
  filter!(r->(!ismissing(r[:adjtotrev])) && (!ismissing(r[:adjnetassets])) &&
    (r[:adjtotrev] ≥ revenueCutoff) && (r[:adjnetassets] ≥ wealthCutoff) , data.df)
  #data.df = data.df[((!ismissing).(data.df[:adjtotrev])) .&
  #  (data.df[:adjtotrev] .≥ revenueCutoff), :]
  #data.df = data.df[((!ismissing).(data.df[:adjnetassets])) .&
  #  (data.df[:adjnetassets] .≥ wealthCutoff), :]

  println("N4: $(size(data.df,1)) N4-Return: $(sum((!ismissing).(data.df[:return])))")

  return data
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function mdl_converter_ode!(
    device::PSY.DynamicInverter{PSY.AverageConverter, O, IC, DC, P, F},
) where {
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component

    #Obtain inner variables for component
    md = get_inner_vars(device)[md_var]
    mq = get_inner_vars(device)[mq_var]
    Vdc = get_inner_vars(device)[Vdc_var]
    θ_oc = get_inner_vars(device)[θ_oc_var]

    #Transform reference frame to grid reference frame
    m_ri = dq_ri(θ_oc + pi / 2) * [md; mq]

    #Update inner_vars
    get_inner_vars(device)[Vr_cnv_var] = m_ri[R] * Vdc
    get_inner_vars(device)[Vi_cnv_var] = m_ri[I] * Vdc
end


text:
using StaticArrays
using HCubature
using ForwardDiff

using Revise
using BoatDynamics
import BoatDynamics: BoatDynamics as bd
using Test

@testset "BoatDynamics.jl" begin
    include("kinetic_energy_test.jl")
    include("approximate_inertia_test.jl")
end


text:
module Yields

import Interpolations
import ForwardDiff
using LinearAlgebra
using UnicodePlots

# don't export type, as the API of Yields.Zero is nicer and
# less polluting than Zero and less/equally verbose as ZeroYieldCurve or ZeroCurve
export rate, discount, accumulation, forward, Yield, Rate, rate, spot

abstract type CompoundingFrequency end
Base.Broadcast.broadcastable(x::T) where {T<:CompoundingFrequency} = Ref(x)

"""
    Continuous()

A type representing continuous interest compounding frequency.

# Examples

```julia-repl
julia> Rate(0.01,Continuous())
Rate(0.01, Continuous())
```

See also: [`Periodic`](@ref)
"""
struct Continuous <: CompoundingFrequency end

"""

```julia-repl
julia> Continuous(0.01)
Rate(0.01, Continuous())
```

See also: [`Periodic`](@ref)
"""
Continuous(rate) = Rate(rate, Continuous())

"""
    Periodic(frequency)

A type representing periodic interest compounding with the given frequency

# Examples

Creating a semi-annual bond equivalent yield:

```julia-repl
julia> Rate(0.01,Periodic(2))
Rate(0.01, Periodic(2))
```

See also: [`Continuous`](@ref)
"""
struct Periodic <: CompoundingFrequency
    frequency::Int
end

"""
    Periodic(rate,frequency)

A convenience constructor for Rate(rate,Periodic(frequency)).

# Examples

Creating a semi-annual bond equivalent yield:

```julia-repl
julia> Periodic(0.01,2)
Rate(0.01, Periodic(2))
```

See also: [`Continuous`](@ref)
"""
Periodic(x, frequency) = Rate(x, Periodic(frequency))

struct Rate{N<:Real,T<:CompoundingFrequency}
    value::N
    compounding::T
end

# Base.:==(r1::Rate,r2::Rate) = (r1.value == r2.value) && (r1.compounding == r2.compounding)

"""
    Rate(rate[,frequency=1])
    Rate(rate,frequency::CompoundingFrequency)

Rate is a type that encapsulates an interest `rate` along with its compounding `frequency`.

Periodic rates can be constructed via `Rate(rate,frequency)` or `Rate(rate,Periodic(frequency))`.

Continuous rates can be constructed via `Rate(rate, Inf)` or `Rate(rate,Continuous())`.

# Examples

```julia-repl
julia> Rate(0.01,Continuous())
Rate(0.01, Continuous())

julia> Rate(0.01,Periodic(2))
Rate(0.01, Periodic(2))

julia> Rate(0.01)
Rate(0.01, Periodic(1))

julia> Rate(0.01,2)
Rate(0.01, Periodic(2))

julia> Rate(0.01,Periodic(4))
Rate(0.01, Periodic(4))

julia> Rate(0.01,Inf)
Rate(0.01, Continuous())

julia> Rate(0.01,Continuous())
Rate(0.01, Continuous())
```
"""
Rate(rate) = Rate(rate, Periodic(1))
Rate(x, frequency::T) where {T<:Real} = isinf(frequency) ? Rate(x, Continuous()) : Rate(x, Periodic(frequency))

"""
    convert(T::CompoundingFrequency,r::Rate)

Returns a `Rate` with an equivalent discount but represented with a different compounding frequency.

# Examples

```julia-repl
julia> r = Rate(Periodic(12),0.01)
Rate(0.01, Periodic(12))

julia> convert(Periodic(1),r)
Rate(0.010045960887181016, Periodic(1))

julia> convert(Continuous(),r)
Rate(0.009995835646701251, Continuous())
```
"""
function Base.convert(T::CompoundingFrequency, r::Rate{<:Real,<:CompoundingFrequency})
    convert(T, r, r.compounding)
end
function Base.convert(to::Continuous, r, from::Continuous)
    return r
end

function Base.convert(to::Periodic, r, from::Continuous)
    return Rate(to.frequency * (exp(r.value / to.frequency) - 1), to)
end

function Base.convert(to::Continuous, r, from::Periodic)
    return Rate(from.frequency * log(1 + r.value / from.frequency), to)
end

function Base.convert(to::Periodic, r, from::Periodic)
    c = convert(Continuous(), r, from)
    return convert(to, c, Continuous())
end

function rate(r::Rate{<:Real,<:CompoundingFrequency})
    r.value
end

function Base.isapprox(a::Rate{N,T}, b::Rate{N,T}; atol::Real = 0, rtol::Real = atol > 0 ? 0 : √eps()) where {T<:Periodic,N<:Real}
    return (a.compounding.frequency == b.compounding.frequency) && isapprox(rate(a), rate(b); atol, rtol)
end

function Base.isapprox(a::Rate{N,T}, b::Rate{N,T}; atol::Real = 0, rtol::Real = atol > 0 ? 0 : √eps()) where {T<:Continuous,N<:Real}
    return isapprox(rate(a), rate(b); atol, rtol)
end

# the fallback for rates not of the same type
function Base.isapprox(a::T, b::N; atol::Real = 0, rtol::Real = atol > 0 ? 0 : √eps()) where {T<:Rate,N<:Rate}
    return isapprox(convert(b.compounding, a), b; atol, rtol)
end

"""
An AbstractYield is an object which can be used as an argument to:

- zero-coupon spot rates via [`zero`](@ref)
- discount factor via [`discount`](@ref)
- accumulation factor via [`accumulation`](@ref)

It can be be constructed via:

- zero rate curve with [`Zero`](@ref)
- forward rate curve with [`Forward`](@ref)
- par rate curve with [`Par`](@ref)
- typical OIS curve with [`OIS`](@ref)
- typical constant maturity treasury (CMT) curve with [`CMT`](@ref)
"""
abstract type AbstractYield end

# make interest curve broadcastable so that you can broadcast over multiple`time`s in `interest_rate`
Base.Broadcast.broadcastable(ic::T) where {T<:AbstractYield} = Ref(ic)

struct YieldCurve{T,U,V} <: AbstractYield
    rates::T
    maturities::U
    zero::V # function time -> continuous zero rate
end

# Forward curves

"""
    ForwardStarting(curve,forwardstart)

Rebase a `curve` so that `discount`/`accumulation`/etc. are re-based so that time zero from the new curves perspective is the given `forwardstart` time.

# Examples

```julia-repl
julia> zero = [5.0, 5.8, 6.4, 6.8] ./ 100
julia> maturity = [0.5, 1.0, 1.5, 2.0]
julia> curve = Yields.Zero(zero, maturity)
julia> fwd = Yields.ForwardStarting(curve, 1.0)

julia> discount(curve,1,2)
0.9275624570410582

julia> discount(fwd,1) # `curve` has effectively been reindexed to `1.0`
0.9275624570410582
```

# Extended Help

While `ForwardStarting` could be nested so that, e.g. the third period's curve is the one-period forward of the second period's curve, it will be more efficient to reuse the initial curve from a runtime and compiler perspective.

`ForwardStarting` is not used to construct a curve based on forward rates. See  [`Forward`](@ref) instead.
"""
struct ForwardStarting{T,U} <: AbstractYield
    curve::U
    forwardstart::T
end

function discount(c::ForwardStarting, to)
    discount(c.curve, c.forwardstart, to + c.forwardstart)
end

"""
    zero(curve,time)
    zero(curve,time,CompoundingFrequency)

Return the zero rate for the curve at the given time. If not specified, will use `Periodic(1)` compounding.
"""
Base.zero(c::YieldCurve, time) = convert(Periodic(1), Continuous(c.zero(time)))
function Base.zero(c::YieldCurve, time, cf::Periodic)
    z = Continuous(c.zero(time))
    i = convert(cf, z)
    return i
end

function Base.zero(c::YieldCurve, time, cf::Continuous)
    return Continuous(c.zero(time))
end

"""
    Constant(rate::Real, cf::CompoundingFrequency=Periodic(1))
    Constant(r::Rate)

Construct a yield object where the spot rate is constant for all maturities. If `rate` is not a `Rate` type, will assume `Periodic(1)` for the compounding frequency

# Examples

```julia-repl
julia> y = Yields.Constant(0.05)
julia> discount(y,2)
0.9070294784580498     # 1 / (1.05) ^ 2
```
"""
struct Constant{T} <: AbstractYield
    rate::T
end

function Constant(rate::T, cf::C = Periodic(1)) where {T<:Real,C<:CompoundingFrequency}
    return Constant(Rate(rate, cf))
end

Base.zero(c::Constant, time) = c.rate
Base.zero(c::Constant, time, cf::C) where {C<:CompoundingFrequency} = convert(cf, c.rate)
rate(c::Constant) = c.rate
rate(c::Constant, time) = c.rate
discount(c::T, time) where {T<:Real} = discount(Constant(c), time)
discount(r::Constant, time) = 1 / accumulation(r, time)

accumulation(r::Constant, time) = accumulation(r.rate.compounding, r, time)
accumulation(c::T, time) where {T>:Real} = accumulation(Constant(c), time)
accumulation(::Continuous, r::Constant, time) = exp(rate(r.rate) * time)
accumulation(::Periodic, r::Constant, time) = (1 + rate(r.rate) / r.rate.compounding.frequency)^(r.rate.compounding.frequency * time)

"""
    Step(rates,times)

Create a yield curve object where the applicable rate is the effective rate of interest applicable until corresponding time. If `rates` is not a `Vector{Rate}`, will assume `Periodic(1)` type.

# Examples

```julia-repl
julia>y = Yields.Step([0.02,0.05], [1,2])

julia>rate(y,0.5)
0.02

julia>rate(y,1.5)
0.05

julia>rate(y,2.5)
0.05
```
"""
struct Step <: AbstractYield
    rates
    times
end

Step(rates) = Step(rates, collect(1:length(rates)))

function rate(y::Step, time)
    i = findfirst(t -> time <= t, y.times)
    if isnothing(i)
        return y.rates[end]
    else
        return y.rates[i]
    end
end

function discount(y::Step, time)
    v = 1 / (1 + y.rates[1])^min(y.times[1], time)

    if y.times[1] >= time
        return v
    end

    for i = 2:length(y.times)

        if y.times[i] >= time
            # take partial discount and break
            v /= (1 + y.rates[i])^(time - y.times[i-1])
            break
        else
            # take full discount and continue
            v /= (1 + y.rates[i])^(y.times[i] - y.times[i-1])
        end

    end

    return v
end

"""
    Zero(rates,maturities)

Construct a yield curve with given zero-coupon spot `rates` at the given `maturities`. If `rates` is not a `Vector{Rate}`, will assume `Periodic(1)` type.
"""
function Zero(rates::Vector{<:Rate}, maturities)
    # bump to a constant yield if only given one rate
    length(rates) == 1 && return Constant(first(rates))

    continuous_zeros = rate.(convert.(Continuous(), rates))
    return YieldCurve(
        rates,
        maturities,
        linear_interp([0.0; maturities], [first(continuous_zeros); continuous_zeros])
    )
end

#fallback if `rates` aren't `Rate`s. Assume `Periodic(1)` per Zero docstring
function Zero(rates, maturities)
    Zero(Periodic.(rates, 1), maturities)
end

function Zero(rates)
    # bump to a constant yield if only given one rate
    maturities = collect(1:length(rates))
    return Zero(rates, maturities)
end

"""
    Par(rate, maturity)

Construct a curve given a set of bond equivalent yields and the corresponding maturities. Assumes that maturities <= 1 year do not pay coupons and that after one year, pays coupons with frequency equal to the CompoundingFrequency of the corresponding rate.

# Examples

```julia-repl

julia> par = [6.,8.,9.5,10.5,11.0,11.25,11.38,11.44,11.48,11.5] ./ 100
julia> maturities = [t for t in 1:10]
julia> curve = Par(par,maturities);
julia> zero(curve,1)
Rate(0.06000000000000005, Periodic(1))

```
"""
function Par(rates::Vector{<:Rate}, maturities)
    # bump to a constant yield if only given one rate
    if length(rates) == 1
        return Constant(rate[1])
    end
    return YieldCurve(
        rates,
        maturities,
        # assume that maturities less than or equal to 12 months are settled once, otherwise semi-annual
        # per Hull 4.7
        bootstrap(rates, maturities, [m <= 1 ? nothing : 1 / r.compounding.frequency for (r, m) in zip(rates, maturities)])
    )
end

function Par(rates::Vector{T}, maturities) where {T<:Real}
    return Par(Rate.(rates), maturities)
end

"""
    Forward(rate_vector,maturities)

Takes a vector of 1-period forward rates and constructs a discount curve.
"""
function Forward(rates, maturities)
    # convert to zeros and pass to Zero
    disc_v = similar(rates)
    v = 1.0

    for i = 1:length(rates)
        Δt = maturities[i] - (i == 1 ? 0 : maturities[i-1])
        v *= discount(Constant(rates[i]), Δt)
        disc_v[i] = v
    end

    z = (1.0 ./ disc_v) .^ (1 ./ maturities) .- 1 # convert disc_v to zero
    return Zero(z, maturities)
end

Forward(rates) = Forward(rates, collect(1:length(rates)))

"""
Takes CMT yields (bond equivalent), and assumes that instruments <= one year maturity pay no coupons and that the rest pay semi-annual.
"""
function CMT(rates::Vector{T}, maturities) where {T<:Real}
    rs = map(zip(rates, maturities)) do (r, m)
        if m <= 1
            Rate(r, Periodic(1 / m))
        else
            Rate(r, Periodic(2))
        end
    end

    CMT(rs, maturities)
end

function CMT(rates::Vector{<:Rate}, maturities)
    return YieldCurve(
        rates,
        maturities,
        # assume that maturities less than or equal to 12 months are settled once, otherwise semi-annual
        # per Hull 4.7
        bootstrap(rates, maturities, [m <= 1 ? nothing : 0.5 for m in maturities])
    )
end

"""
    OIS(rates,maturities)
Takes Overnight Index Swap rates, and assumes that instruments <= one year maturity are settled once and other agreements are settled quarterly with a corresponding CompoundingFrequency

"""
function OIS(rates::Vector{T}, maturities) where {T<:Real}
    rs = map(zip(rates, maturities)) do (r, m)
        if m <= 1
            Rate(r, Periodic(1 / m))
        else
            Rate(r, Periodic(4))
        end
    end

    return OIS(rs, maturities)
end
function OIS(rates::Vector{<:Rate}, maturities)
    return YieldCurve(
        rates,
        maturities,
        # assume that maturities less than or equal to 12 months are settled once, otherwise quarterly
        # per Hull 4.7
        bootstrap(rates, maturities, [m <= 1 ? nothing : 1 / 4 for m in maturities])
    )
end

abstract type ObservableQuote end

"""
    ZeroCouponQuote(price, maturity)

Quote for a set of zero coupon bonds with given `price` and `maturity`.

# Examples

```julia-repl
julia> prices = [1.3, 0.1, 4.5]
julia> maturities = [1.2, 2.5, 3.6]
julia> swq = Yields.ZeroCouponQuote.(prices, maturities)
```
"""
struct ZeroCouponQuote <: ObservableQuote
    price
    maturity
end

"""
    SwapQuote(yield, maturity, frequency)

Quote for a set of interest rate swaps with the given `yield` and `maturity` and a given payment `frequency`.

# Examples

```julia-repl
julia> maturities = [1.2, 2.5, 3.6]
julia> interests = [-0.02, 0.3, 0.04]
julia> prices = [1.3, 0.1, 4.5]
julia> frequencies = [2,1,2]
julia> swq = Yields.SwapQuote.(interests, maturities, frequencies)
```
"""
struct SwapQuote <: ObservableQuote
    yield
    maturity
    frequency
    function SwapQuote(yield, maturity, frequency)
        frequency <= 0 && throw(DomainError("Payment frequency must be positive"))
        return new(yield, maturity, frequency)
    end
end

"""
    BulletBondQuote(yield, price, maturity, frequency)

Quote for a set of fixed interest bullet bonds with given `yield`, `price`, `maturity` and a given payment frequency `frequency`.

Construct a vector of quotes for use with SmithWilson methods, e.g. by broadcasting over an array of inputs.

# Examples

```julia-repl
julia> maturities = [1.2, 2.5, 3.6]
julia> interests = [-0.02, 0.3, 0.04]
julia> prices = [1.3, 0.1, 4.5]
julia> frequencies = [2,1,2]
julia> bbq = Yields.BulletBondQuote.(interests, maturities, prices, frequencies)
```
"""
struct BulletBondQuote <: ObservableQuote
    yield
    price
    maturity
    frequency

    function BulletBondQuote(yield, maturity, price, frequency)
        frequency <= 0 && throw(DomainError("Payment frequency must be positive"))
        return new(yield, maturity, price, frequency)
    end
end

"""
    SmithWilson(zcq::Vector{ZeroCouponQuote}; ufr, α)
    SmithWilson(swq::Vector{SwapQuote}; ufr, α)
    SmithWilson(bbq::Vector{BulletBondQuote}; ufr, α)
    SmithWilson(times<:AbstractVector, cashflows<:AbstractMatrix, prices<:AbstractVector; ufr, α)
    SmithWilson(u, qb; ufr, α)

Create a yield curve object that implements the Smith-Wilson interpolation/extrapolation scheme.

Positional arguments to construct a curve:

- Quoted instrument as the first argument: either a `Vector` of `ZeroCouponQuote`s, `SwapQuote`s, or `BulletBondQuote`s, or
- A set of `times`, `cashflows`, and `prices`, or
- A curve can be with `u` is the timepoints coming from the calibration, and `qb` is the internal parameterization of the curve that ensures that the calibration is correct. Users may prefer the other constructors but this mathematical constructor is also available.

Required keyword arguments:

- `ufr` is the Ultimate Forward Rate, the forward interest rate to which the yield curve tends, in continuous compounding convention.
- `α` is the parameter that governs the speed of convergence towards the Ultimate Forward Rate. It can be typed with `\\alpha[TAB]`
"""
struct SmithWilson{TU<:AbstractVector,TQb<:AbstractVector} <: AbstractYield
    u::TU
    qb::TQb
    ufr
    α

    # Inner constructor ensures that vector lengths match
    function SmithWilson{TU,TQb}(u, qb; ufr, α) where {TU<:AbstractVector,TQb<:AbstractVector}
        if length(u) != length(qb)
            throw(DomainError("Vectors u and qb in SmithWilson must have equal length"))
        end
        return new(u, qb, ufr, α)
    end
end

SmithWilson(u::TU, qb::TQb; ufr, α) where {TU<:AbstractVector,TQb<:AbstractVector} = SmithWilson{TU,TQb}(u, qb; ufr = ufr, α = α)

"""
    H_ordered(α, t_min, t_max)

The Smith-Wilson H function with ordered arguments (for better performance than using min and max).
"""
function H_ordered(α, t_min, t_max)
    return α * t_min + exp(-α * t_max) * sinh(-α * t_min)
end

"""
    H(α, t1, t2)

The Smith-Wilson H function implemented in a faster way.
"""
function H(α, t1::T, t2::T) where {T}
    return t1 < t2 ? H_ordered(α, t1, t2) : H_ordered(α, t2, t1)
end

H(α, t1, t2) = H(α, promote(t1, t2)...)

H(α, t1vec::AbstractVector, t2) = [H(α, t1, t2) for t1 in t1vec]
H(α, t1vec::AbstractVector, t2vec::AbstractVector) = [H(α, t1, t2) for t1 in t1vec, t2 in t2vec]
# This can be optimized by going to H_ordered directly, but it might be a bit cumbersome
H(α, tvec::AbstractVector) = H(α, tvec, tvec)

discount(sw::SmithWilson, t) = exp(-sw.ufr * t) * (1.0 + H(sw.α, sw.u, t) ⋅ sw.qb)
Base.zero(sw::SmithWilson, t) = Continuous(sw.ufr - log(1.0 + H(sw.α, sw.u, t) ⋅ sw.qb) / t)
Base.zero(sw::SmithWilson, t, cf::CompoundingFrequency) = convert(cf, zero(sw, t))

function SmithWilson(times::AbstractVector, cashflows::AbstractMatrix, prices::AbstractVector; ufr, α)
    Q = Diagonal(exp.(-ufr * times)) * cashflows
    q = vec(sum(Q, dims = 1))  # We want q to be a column vector
    QHQ = Q' * H(α, times) * Q
    b = QHQ \ (prices - q)
    Qb = Q * b
    return SmithWilson(times, Qb; ufr = ufr, α = α)
end

"""
    timepoints(zcq::Vector{ZeroCouponQuote})
    timepoints(bbq::Vector{BulletBondQuote})

Return the times associated with the `cashflows` of the instruments.
"""
function timepoints(qs::Vector{Q}) where {Q<:ObservableQuote}
    frequency = maximum(q.frequency for q in qs)
    timestep = 1 / frequency
    maturity = maximum(q.maturity for q in qs)
    return [timestep:timestep:maturity...]
end

"""
    cashflows(interests, maturities, frequency)
    timepoints(zcq::Vector{ZeroCouponQuote})
    timepoints(bbq::Vector{BulletBondQuote})

Produce a cash flow matrix for a set of instruments with given `interests` and `maturities`
and a given payment frequency `frequency`. All instruments are assumed to have their first payment at time 1/`frequency`
and have their last payment at the largest multiple of 1/`frequency` less than or equal to the input maturity.
"""
function cashflows(interests, maturities, frequencies)
    frequency = lcm(frequencies)
    fq = inv.(frequencies)
    timestep = 1 / frequency
    floored_mats = floor.(maturities ./ timestep) .* timestep
    times = timestep:timestep:maximum(floored_mats)
    # we need to determine the coupons in relation to the payment date, not time zero
    time_adj = floored_mats .% fq

    cashflows = [
        # if on a coupon date and less than maturity, pay coupon
        ((((t + time_adj[instrument]) % fq[instrument] ≈ 0) && t <= floored_mats[instrument]) ? interests[instrument] / frequencies[instrument] : 0.0) +
        (t ≈ floored_mats[instrument] ? 1.0 : 0.0) # add maturity payment
        for t in times, instrument = 1:length(interests)
    ]

    return cashflows
end

function cashflows(qs::Vector{Q}) where {Q<:ObservableQuote}
    yield = [q.yield for q in qs]
    maturity = [q.maturity for q in qs]
    frequency = [q.frequency for q in qs]
    return cashflows(yield, maturity, frequency)
end

# Utility methods for calibrating Smith-Wilson directly from quotes
function SmithWilson(zcq::Vector{ZeroCouponQuote}; ufr, α)
    n = length(zcq)
    maturities = [q.maturity for q in zcq]
    prices = [q.price for q in zcq]
    return SmithWilson(maturities, Matrix{Float64}(I, n, n), prices; ufr = ufr, α = α)
end

function SmithWilson(swq::Vector{SwapQuote}; ufr, α)
    times = timepoints(swq)
    cfs = cashflows(swq)
    ones(length(swq))
    return SmithWilson(times, cfs, ones(length(swq)), ufr = ufr, α = α)
end

function SmithWilson(bbq::Vector{BulletBondQuote}; ufr, α)
    times = timepoints(bbq)
    cfs = cashflows(bbq)
    prices = [q.price for q in bbq]
    return SmithWilson(times, cfs, prices, ufr = ufr, α = α)
end

# https://github.com/dpsanders/hands_on_julia/blob/master/during_sessions/Fractale%20de%20Newton.ipynb
newton(f, f′, x) = x - f(x) / f′(x)
function solve(g, g′, x0, max_iterations = 100)
    x = x0

    tolerance = 2 * eps(x0)
    iteration = 0

    while (abs(g(x) - 0) > tolerance && iteration < max_iterations)
        x = newton(g, g′, x)
        iteration += 1
    end

    return x
end

function bootstrap(rates, maturities, settlement_frequency; interp_function = linear_interp)
    discount_vec = zeros(length(rates)) # construct a placeholder discount vector matching maturities
    # we have to take the first rate as the starting point
    discount_vec[1] = discount(Constant(rates[1]), maturities[1])

    for t = 2:length(maturities)
        if isnothing(settlement_frequency[t])
            # no settlement before maturity
            discount_vec[t] = discount(Constant(rates[t]), maturities[t])
        else
            # need to account for the interim cashflows settled
            times = settlement_frequency[t]:settlement_frequency[t]:maturities[t]
            cfs = [rate(rates[t]) * settlement_frequency[t] for s in times]
            cfs[end] += 1

            function pv(v_guess)
                v = interp_function([[0.0]; maturities[1:t]], vcat(1.0, discount_vec[1:t-1], v_guess...))
                return sum(v.(times) .* cfs)
            end
            target_pv = sum(map(t2 -> discount(Constant(rates[t]), t2), times) .* cfs)
            root_func(v_guess) = pv(v_guess) - target_pv
            root_func′(v_guess) = ForwardDiff.derivative(root_func, v_guess)
            discount_vec[t] = solve(root_func, root_func′, rate(rates[t]))
        end

    end
    zero_vec = -log.(discount_vec) ./ maturities
    return linear_interp([0.0; maturities], [first(zero_vec); zero_vec])
end

## Generic and Fallbacks
"""
    discount(rate,to)
    discount(rate,from,to)

The discount factor for the `rate` for times `from` through `to`. If rate is a `Real` number, will assume a `Constant` interest rate.
"""
discount(yc, time) = exp(-yc.zero(time) * time)
discount(rate::Rate{<:Real,<:CompoundingFrequency}, from, to) = discount(Constant(rate), from, to)
discount(rate::Rate{<:Real,<:CompoundingFrequency}, to) = discount(Constant(rate), to)

discount(yc, from, to) = discount(yc, to) / discount(yc, from)

"""
    forward(curve,from,to,CompoundingFrequency=Periodic(1))

The forward `Rate` implied by the curve between times `from` and `to`.
"""
function forward(yc, from, to)
    return forward(yc, from, to, Periodic(1))
end

function forward(yc, from, to, cf::T) where {T<:CompoundingFrequency}

    r = Periodic((accumulation(yc, to) / accumulation(yc, from))^(1 / (to - from)) - 1, 1)
    return convert(cf, r)
end

function forward(yc, from)
    to = from + 1
    return forward(yc, from, to)
end

"""
    accumulation(rate,from,to)

The accumulation factor for the `rate` for times `from` through `to`. If rate is a `Real` number, will assume a `Constant` interest rate.
"""
function accumulation(y::T, time) where {T<:AbstractYield}
    return 1 ./ discount(y, time)
end
accumulation(rate::Rate{<:Real,<:CompoundingFrequency}, to) = accumulation(Constant(rate), to)

function accumulation(y::T, from, to) where {T<:AbstractYield}
    return 1 ./ discount(y, from, to)
end
accumulation(rate::Rate{<:Real,<:CompoundingFrequency}, from, to) = accumulation(Constant(rate), from, to)

## Curve Manipulations
struct RateCombination{T,U,V} <: AbstractYield
    r1::T
    r2::U
    op::V
end

rate(rc::RateCombination, time) = rc.op(rate(rc.r1, time), rate(rc.r2, time))
function discount(rc::RateCombination, time)
    a1 = discount(rc.r1, time)^(-1 / time) - 1
    a2 = discount(rc.r2, time)^(-1 / time) - 1
    return 1 / (1 + rc.op(a1, a2))^time
end

"""
    Yields.AbstractYield + Yields.AbstractYield

The addition of two yields will create a `RateCombination`. For `rate`, `discount`, and `accumulation` purposes the spot rates of the two curves will be added together.
"""
function Base.:+(a::AbstractYield, b::AbstractYield)
    return RateCombination(a, b, +)
end

function Base.:+(a::Constant, b::Constant)
    a_kind = rate(a).compounding
    rate_new_basis = rate(convert(a_kind, rate(b)))
    return Constant(
        Rate(
            rate(a.rate) + rate_new_basis,
            a_kind
        )
    )
end

function Base.:+(a::T, b) where {T<:AbstractYield}
    return a + Constant(b)
end

function Base.:+(a, b::T) where {T<:AbstractYield}
    return Constant(a) + b
end

"""
    Yields.AbstractYield - Yields.AbstractYield

The subtraction of two yields will create a `RateCombination`. For `rate`, `discount`, and `accumulation` purposes the spot rates of the second curves will be subtracted from the first.
"""
function Base.:-(a::AbstractYield, b::AbstractYield)
    return RateCombination(a, b, -)
end

function Base.:-(a::Constant, b::Constant)
    a_kind = rate(a).compounding
    rate_new_basis = rate(convert(a_kind, rate(b)))
    return Constant(
        Rate(
            rate(a.rate) - rate_new_basis,
            a_kind
        )
    )
end

function Base.:-(a::T, b) where {T<:AbstractYield}
    return a - Constant(b)
end

function Base.:-(a, b::T) where {T<:AbstractYield}
    return Constant(a) - b
end

linear_interp(xs, ys) = Interpolations.extrapolate(
    Interpolations.interpolate((xs,), ys, Interpolations.Gridded(Interpolations.Linear())),
    Interpolations.Line()
)

# used to display simple type name in show method
# https://stackoverflow.com/questions/70043313/get-simple-name-of-type-in-julia?noredirect=1#comment123823820_70043313
name(::Type{T}) where {T} = (isempty(T.parameters) ? T : T.name.wrapper)

function Base.show(io::IO, curve::T) where {T<:AbstractYield}
    println() # blank line for padding
    r = zero(curve, 1)
    ylabel = isa(r.compounding, Continuous) ? "Continuous" : "Periodic($(r.compounding.frequency))"
    kind = name(typeof(curve))
    l = lineplot(
        t -> rate(zero(curve, t)),
        0.0, #from
        30.0,  # to
        xlabel = "time",
        ylabel = ylabel,
        compact = true,
        name = "Zero rates",
        width = 60,
        title = "Yield Curve ($kind)"
    )
    show(io, l)
end

end


text:
include("../src/autodiff_module.jl")
include("../../__lib__/math/common/numder/src/numder_module.jl")

module dtest
using PyPlot
PyPlot.pygui(true)
using ..autodiff
using ..numder

x1 = autodiff.DualNumber(1.0)
x2 = autodiff.DualNumber(2.0)
x3 = autodiff.DualNumber(pi/2)

x4 = x1 * x2
x5 = sin(x3)
x6 = exp(x4)
x7 = x4 * x5
x8 = x6 + x7
x9 = x8 / x3

fcn(x) = begin
    x1 = x[1]
    x2 = x[2]
    x3 = x[3]
    (x1 * x2 * sin(x3) + exp(x1 * x2)) / x3
end

x_vec = [1.0, 2, pi/2]

L        = length(x_vec)
grad_vec = zeros(L)
e_i      = zeros(L)
for kk = 1:L
    e_i[kk]         = 1.0
    y_vec           = [autodiff.DualNumber(x_vec[kk], e_i[kk]) for kk = 1:L]
    grad_vec[kk]    = fcn(y_vec).v
    e_i[kk]         = 0.0
end

function test_elem()

    fcn(x)      = 3 * x^2 - 10 * x + 10 # log(x) # cos(x) # sin(x) # x^2
    theodiff(x) = 6 * x - 10 # 1 / x  # -sin(x) # cos(x)

    xs = collect( 0.0:0.1:3 )
    ds = autodiff.DualNumber.(xs)

    fs = fcn.(xs)
    theograd = theodiff.(xs)
    gs = fcn.(ds)

    f_from_d    = autodiff.value(gs)
    dfdx_from_d = autodiff.derivative(gs)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(xs, fs, marker = :s)
    PyPlot.plot(xs, f_from_d, marker = :*)
    PyPlot.plot(xs, dfdx_from_d, marker = :o, markersize = 10)
    PyPlot.plot(xs, theograd, marker = :^, markersize = 10)

    # return (gs, f_from_d, dfdx_from_d)
end

function test_dual()

    add(x1, x2) = x1 + x2
    mult(x1, x2) = x1 * x2
    div(x1, x2) = x1 / x2

    addv(x) = x[1] + x[2]
    multv(x) = x[1] * x[2]
    divv(x) = x[1] / x[2]

    f1 = 3.0
    f2 = 2.0

    d1 = autodiff.DualNumber(2.0, 1.0)
    d2 = autodiff.DualNumber(3.0, 0.0)
    dp = autodiff.DualNumber(5.0, 1.0)
    d3 = autodiff.DualNumber(-2.0, -1.5)
    d4 = autodiff.DualNumber(-2.0, 1.5)

    pd = add(d1, d2)
    p1 = add(d1, f1)
    p2 = add(f2, d2)

    md = mult(d1, d2)
    m1 = mult(d1, f1)
    m2 = mult(f2, d2)

    n_md = numder.gradient_fw(multv, [d1.x, d2.x])

    did = div(d1, d2)
    did2 = div(d2, d1)

    di1 = div(d1, f1)
    di2 = div(f2, d2)

    n_dd = numder.gradient_fw(divv, [d1.x, d2.x])

    sq = d1^d2

    a1 = abs(d3)
    a2 = abs(d4)

    s = sin(d1)
    c = cos(d1)
    t = tan(d1)

    sh = sinh(d1)
    ch = cosh(d1)
    th = tanh(d1)

    e = exp(d1)
    l = log(d1)

    r = sqrt(d1)

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
function getSDPUpperBound_gd(Sigma::Array{Float64, 2}, k::Int64, useSOCS::Bool=false, usePSDs::Bool=true)
# Note: greedy in terms of the rounding mechanism.

    n = size(Sigma, 1)

    sdp_gd = Model(Mosek.Optimizer)
    if usePSDs
        @variable(sdp_gd, X[1:n, 1:n], PSD)
    else
        @variable(sdp_gd, X[1:n, 1:n], Symmetric)
        @constraint(sdp_gd, diagnonneg[i=1:n], X[i,i]>=0.0)
        @constraint(sdp_gd, twobytwominor[i=1:n, j=1:n], [X[i,i]+X[j,j]; X[i,i]-X[j,j]; 2.0*X[i,j]] in SecondOrderCone())
    end
    @variable(sdp_gd, z[1:n]>=0.0)
    @constraint(sdp_gd, z.<=1.0)

    @constraint(sdp_gd, sum(diag(X))==1.0)
    @constraint(sdp_gd, sum(z)<=k)
    @constraint(sdp_gd, diagX[i=1:n], X[i,i]<=z[i])
    @constraint(sdp_gd, offDiagX[i=1:n, j=1:n], (i!=j)*(X[i,j]-0.5*z[i])<=0.0)
    @constraint(sdp_gd, offDiagX2[i=1:n, j=1:n], (i!=j)*(-X[i,j]-0.5*z[i])<=0.0)
    if useSOCS
        @variable(sdp_gd, U[1:n, 1:n])
        @constraint(sdp_gd, sum(U)<=k)
        @constraint(sdp_gd, X.-U.<=0.0)
        @constraint(sdp_gd, X.+U.>=0.0)
        @constraint(sdp_gd, perspectiveRelaxation[i=1:n], [z[i]+X[i,i]; 2*X[:,i];z[i]-X[i,i]] in SecondOrderCone())
    end

    @objective(sdp_gd, Max, LinearAlgebra.dot(Sigma, X))

    @suppress optimize!(sdp_gd)

    #Get greedy solution and evaluate objective
    indices_gd=sortperm(value.(z), rev=true)[1:k]
    Sigma_reduced=Sigma[indices_gd, indices_gd]
    λs,xs=eigs(Sigma_reduced, nev=1)
    λ=λs[1]

    return (JuMP.objective_value(sdp_gd)), λ
end

function getSDPUpperBound_d2008_gd(Sigma::Array{Float64, 2}, k::Int64, useSOCS::Bool=false)
# Note: greedy in terms of the rounding mechanism.

    A=cholesky(Sigma) #Need to use A.U
    n = size(Sigma, 1)

    sdp_gd = Model(SCS.Optimizer)
    if n<=13
        sdp_gd=Model(Mosek.Optimizer)
    end

    @variable(sdp_gd, X[1:n, 1:n], PSD)
    @variable(sdp_gd, P[1:n, 1:n, 1:n]) #let third dimension index each matrix
    @constraint(sdp_gd, imposePSDConstraint[i=1:n], P[:,:,i] in PSDCone())
    @constraint(sdp_gd, imposePSDConstraint2[i=1:n], X.-P[:,:,i] in PSDCone())

    @variable(sdp_gd, z[1:n]>=0.0)
    @constraint(sdp_gd, z.<=1.0)

    @constraint(sdp_gd, sum(diag(X))==1.0)
    @constraint(sdp_gd, imposebinaries[i=1:n], sum(P[j,j,i] for j=1:n)==z[i])
    if useSOCS
        @variable(sdp_gd, U[1:n, 1:n])
        @constraint(sdp_gd, sum(U)<=k)
        @constraint(sdp_gd, X.-U.<=0.0)
        @constraint(sdp_gd, X.+U.>=0.0)
        @constraint(sdp_gd, perspectiveRelaxation[i=1:n], [z[i]+X[i,i]; 2*X[:,i];z[i]-X[i,i]] in SecondOrderCone())
    end

    @objective(sdp_gd, Max, sum(LinearAlgebra.dot(A.U[i,:]*A.U[i,:]', P[:,:,i]) for i=1:n))
    #@show "About to optimize!"
    @suppress optimize!(sdp_gd)
    #@show getobjectivevalue(sdp_gd)

    #Get greedy solution and evaluate objective
    indices_gd=sortperm(value.(z), rev=true)[1:k]
    Sigma_reduced=Sigma[indices_gd, indices_gd]
    λs,xs=eigs(Sigma_reduced, nev=1)
    λ=λs[1]

    return (JuMP.objective_value(sdp_gd)), λ
end

function getSDPUpperBound_gd_highdim(Sigma::Array{Float64, 2}, k::Int64, useSOCS::Bool=false, useCuts=false)
# Note: greedy in terms of the rounding mechanism.

    n = size(Sigma, 1)

    sdp_gd = Model(Mosek.Optimizer)

    @variable(sdp_gd, X[1:n, 1:n], Symmetric)
    @variable(sdp_gd, z[1:n]>=0.0)
    @constraint(sdp_gd, z.<=1.0)

    @constraint(sdp_gd, sum(diag(X))==1.0)
    @constraint(sdp_gd, sum(z)<=k)
    @constraint(sdp_gd, diagX[i=1:n], X[i,i]<=z[i])
    @constraint(sdp_gd, offDiagX[i=1:n, j=1:n], (i!=j)*(X[i,j]-0.5*z[i])<=0.0)
    @constraint(sdp_gd, offDiagX2[i=1:n, j=1:n], (i!=j)*(-X[i,j]-0.5*z[i])<=0.0)
    if useSOCS
        @variable(sdp_gd, U[1:n, 1:n])
        @constraint(sdp_gd, sum(U)<=k)
        @constraint(sdp_gd, X.-U.<=0.0)
        @constraint(sdp_gd, X.+U.>=0.0)
        @constraint(sdp_gd, perspectiveRelaxation[i=1:n], [z[i]+X[i,i]; 2*X[:,i];z[i]-X[i,i]] in SecondOrderCone())
    end

    @objective(sdp_gd, Max, Compat.dot(Sigma, X))

    t_start=time()
    @suppress optimize!(sdp_gd)
    getobjectivevalue(sdp_gd)

    if useCuts
     j=0
     maxCuts=20
         while j<maxCuts
             λ, ϕ=eigs(value.(X),  nev=1, which=:SR)
                if real(λ[1])<=-1e-4
                     cutV=real(ϕ[:, 1])
                     @constraint(sdp_gd, Compat.dot(cutV*cutV', X)>=0.0)
                else
                    break
                end
             @suppress optimize!(sdp_gd)
             j+=1
        end
    end

    #Get greedy solution and evaluate objective
    indices_gd=sortperm(value.(z), rev=true)[1:k]
    Sigma_reduced=Sigma[indices_gd, indices_gd]
    λs,xs=eigs(Sigma_reduced, nev=1)
    λ=λs[1]

    return (getobjectivevalue(sdp_gd)), λ
end


text:
let
    @test include("../docs/make.jl")
end


text:
module BenchMapi

using BenchmarkTools
using ThreadsX

collatz(x) =
    if iseven(x)
        x ÷ 2
    else
        3x + 1
    end

function collatz_stopping_time(x)
    n = 0
    while true
        x == 1 && return n
        n += 1
        x = collatz(x)
    end
end

function consume(ns)
    t0 = time_ns()
    d = 0
    while d < ns
        d = Int(time_ns() - t0)
    end
    return d
end

constant(x) = _ -> x

const SUITE = BenchmarkGroup()

let s0 = SUITE["collatz"] = BenchmarkGroup()
    s0["base"] = @benchmarkable map(collatz_stopping_time, 1:1000_000)
    s0["tx"] =
        @benchmarkable ThreadsX.mapi(collatz_stopping_time, 1:1000_000; basesize = 12500)
end

let s0 = SUITE["consume-1ms"] = BenchmarkGroup()
    s0["base"] = @benchmarkable map(consume ∘ constant(1_000_000), 1:20)
    s0["tx"] = @benchmarkable ThreadsX.mapi(consume ∘ constant(1_000_000), 1:20)
end

let s0 = SUITE["consume-100us"] = BenchmarkGroup()
    s0["base"] = @benchmarkable map(consume ∘ constant(100_000), 1:20)
    s0["tx"] = @benchmarkable ThreadsX.mapi(consume ∘ constant(100_000), 1:20)
end

end  # module
BenchMapi.SUITE